thalfshear-darcy-rate-dependence.py - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
(HTM) git clone git://src.adamsgaard.dk/sphere
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
thalfshear-darcy-rate-dependence.py (16287B)
---
1 #!/usr/bin/env python
2 import matplotlib
3 matplotlib.use('Agg')
4 #matplotlib.rcParams.update({'font.size': 18, 'font.family': 'serif'})
5 import shutil
6
7 import os
8 import sys
9 import numpy
10 import sphere
11 import matplotlib.pyplot as plt
12 import scipy.optimize
13
14 matplotlib.rcParams.update({'font.size': 7, 'font.family': 'sans-serif'})
15 matplotlib.rc('text', usetex=True)
16 matplotlib.rcParams['text.latex.preamble']=[r"\usepackage{amsmath}"]
17
18 sids =\
19 ['halfshear-darcy-sigma0=80000.0-k_c=3.5e-13-mu=1.04e-07-ss=10000.0-A=70000.0-f=0.2',
20 'halfshear-darcy-sigma0=80000.0-k_c=3.5e-15-mu=1.04e-07-ss=10000.0-A=70000.0-f=0.2',
21 'halfshear-darcy-sigma0=10000.0-k_c=2e-16-mu=2.08e-07-ss=2000.0-A=4000.0-f=0.2',
22 'halfshear-darcy-sigma0=10000.0-k_c=2e-16-mu=2.08e-07-ss=2000.0-A=4125.0-f=0.2',
23 'halfshear-darcy-sigma0=10000.0-k_c=2e-16-mu=2.08e-07-ss=2000.0-A=4250.0-f=0.2',
24 'halfshear-darcy-sigma0=10000.0-k_c=2e-16-mu=2.08e-07-ss=2000.0-A=4375.0-f=0.2',
25 'halfshear-darcy-sigma0=10000.0-k_c=2e-16-mu=2.08e-07-ss=2000.0-A=4500.0-f=0.2',
26 'halfshear-darcy-sigma0=10000.0-k_c=2e-16-mu=2.08e-07-ss=2000.0-A=4625.0-f=0.2',
27 'halfshear-darcy-sigma0=10000.0-k_c=2e-16-mu=2.08e-07-ss=2000.0-A=4750.0-f=0.2',
28 'halfshear-darcy-sigma0=10000.0-k_c=2e-16-mu=2.08e-07-ss=2000.0-A=4875.0-f=0.2',
29 'halfshear-darcy-sigma0=10000.0-k_c=2e-16-mu=2.08e-07-ss=2000.0-A=4000.0-f=0.1',
30 'halfshear-darcy-sigma0=10000.0-k_c=2e-16-mu=2.08e-07-ss=2000.0-A=4500.0-f=0.1',
31 'halfshear-darcy-sigma0=10000.0-k_c=2e-16-mu=2.08e-07-ss=3000.0-A=4000.0-f=0.2',
32 'halfshear-darcy-sigma0=10000.0-k_c=2e-16-mu=2.08e-07-ss=3000.0-A=4250.0-f=0.2',
33 'halfshear-darcy-sigma0=10000.0-k_c=2e-16-mu=2.08e-07-ss=3000.0-A=4500.0-f=0.2',
34 ]
35 #['halfshear-darcy-sigma0=80000.0-k_c=3.5e-13-mu=1.797e-06-ss=10000.0-A=70000.0-f=0.2']
36 sids = \
37 ['halfshear-darcy-sigma0=80000.0-k_c=3.5e-13-mu=1.04e-07-ss=10000.0-A=70000.0-f=0.2']
38 outformat = 'pdf'
39 fluid = True
40 #threshold = 100.0 # [N]
41
42
43 def creep_rheology1(friction, n):
44 ''' return strain rate from friction (tau/N) value '''
45 return friction**n
46
47 def creep_rheology2(friction, n, A):
48 ''' return strain rate from friction (tau/N) value '''
49 return A*friction**n
50
51
52 for sid in sids:
53
54 print sid
55 sim = sphere.sim(sid, fluid=fluid)
56
57 #nsteps = 2
58 #nsteps = 10
59 nsteps = sim.status()
60
61 t = numpy.empty(nsteps)
62
63 tau = numpy.empty(sim.status())
64 N = numpy.empty(sim.status())
65 #v = numpy.empty(sim.status())
66 shearstrainrate = numpy.empty(sim.status())
67 shearstrain = numpy.empty(sim.status())
68 dilation = numpy.empty(sim.status())
69 for i in numpy.arange(sim.status()):
70 sim.readstep(i+1, verbose=False)
71 t_DEM_to_t_real = sim.mu[0]/1.787e-3
72 #tau = sim.shearStress()
73 tau[i] = sim.w_tau_x # defined shear stress
74 #tau[i] = sim.shearStress()[0] # measured shear stress along x
75 N[i] = sim.currentNormalStress() # defined normal stress
76 #v[i] = sim.shearVel()
77 shearstrainrate[i] = sim.shearStrainRate()*t_DEM_to_t_real
78 shearstrain[i] = sim.shearStrain()
79 t[i] = sim.currentTime()
80
81 if i == 0:
82 initial_height = sim.w_x[0]
83
84 dilation[i] = sim.w_x[0] - initial_height
85
86 # remove nonzero sliding velocities and their associated values
87 #idx = numpy.nonzero(v)
88 idx = numpy.nonzero(shearstrainrate)
89 #v_nonzero = v[idx]
90 shearstrainrate_nonzero = shearstrainrate[idx]
91 tau_nonzero = tau[idx]
92 N_nonzero = N[idx]
93 shearstrain_nonzero = shearstrain[idx]
94
95 t_nonzero = t[idx]
96
97 ### "Critical" state fit
98 # The algorithm uses the Levenberg-Marquardt algorithm through leastsq
99 idxfit = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
100 (shearstrainrate_nonzero < 0.1*t_DEM_to_t_real) &
101 (((t_nonzero > 6.0) & (t_nonzero < 8.0)) |
102 ((t_nonzero > 11.0) & (t_nonzero < 13.0)) |
103 ((t_nonzero > 16.0) & (t_nonzero < 18.0)) |
104 ((t_nonzero > 21.0) & (t_nonzero < 23.0)) |
105 ((t_nonzero > 26.0) & (t_nonzero < 28.0)) |
106 ((t_nonzero > 31.0) & (t_nonzero < 33.0)) |
107 ((t_nonzero > 36.0) & (t_nonzero < 38.0)) |
108 ((t_nonzero > 41.0) & (t_nonzero < 43.0)) |
109 ((t_nonzero > 46.0) & (t_nonzero < 48.0)) |
110 ((t_nonzero > 51.0) & (t_nonzero < 53.0)) |
111 ((t_nonzero > 56.0) & (t_nonzero < 58.0))))
112
113 #popt, pvoc = scipy.optimize.curve_fit(
114 #creep_rheology, tau/N, shearstrainrate)
115 popt, pvoc = scipy.optimize.curve_fit(
116 creep_rheology1, tau_nonzero[idxfit]/N_nonzero[idxfit],
117 shearstrainrate_nonzero[idxfit])
118 popt2, pvoc2 = scipy.optimize.curve_fit(
119 creep_rheology2, tau_nonzero[idxfit]/N_nonzero[idxfit],
120 shearstrainrate_nonzero[idxfit])
121 print '# Creep rheology 1'
122 print popt
123 print pvoc
124 n = popt[0] # stress exponent
125 #A = popt[1] # stress exponent
126 A = 1.
127
128 print '# Creep rheology 2'
129 print popt2
130 print pvoc2
131 n = popt2[0] # stress exponent
132 A = popt2[1] # stress exponent
133 #A = 1.
134
135
136 '''
137 # Investigate n and A evolution from cycle to cycle
138 idxfit_ = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
139 (shearstrainrate_nonzero < 0.1) &
140 (((t_nonzero > 1.0) & (t_nonzero < 3.0))))
141 popt_0, pvoc_0 = scipy.optimize.curve_fit(
142 creep_rheology2, tau_nonzero[idxfit_]/N_nonzero[idxfit_],
143 shearstrainrate_nonzero[idxfit_])
144
145 idxfit_ = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
146 (shearstrainrate_nonzero < 0.1) &
147 (((t_nonzero > 6.0) & (t_nonzero < 8.0))))
148 popt_1, pvoc_1 = scipy.optimize.curve_fit(
149 creep_rheology2, tau_nonzero[idxfit_]/N_nonzero[idxfit_],
150 shearstrainrate_nonzero[idxfit_])
151
152 idxfit_ = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
153 (shearstrainrate_nonzero < 0.1) &
154 (((t_nonzero > 11.0) & (t_nonzero < 13.0))))
155 popt_2, pvoc_2 = scipy.optimize.curve_fit(
156 creep_rheology2, tau_nonzero[idxfit_]/N_nonzero[idxfit_],
157 shearstrainrate_nonzero[idxfit_])
158
159 idxfit_ = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
160 (shearstrainrate_nonzero < 0.1) &
161 (((t_nonzero > 16.0) & (t_nonzero < 18.0))))
162 popt_3, pvoc_3 = scipy.optimize.curve_fit(
163 creep_rheology2, tau_nonzero[idxfit_]/N_nonzero[idxfit_],
164 shearstrainrate_nonzero[idxfit_])
165
166 idxfit_ = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
167 (shearstrainrate_nonzero < 0.1) &
168 (((t_nonzero > 21.0) & (t_nonzero < 23.0))))
169 popt_4, pvoc_4 = scipy.optimize.curve_fit(
170 creep_rheology2, tau_nonzero[idxfit_]/N_nonzero[idxfit_],
171 shearstrainrate_nonzero[idxfit_])
172
173 idxfit_ = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
174 (shearstrainrate_nonzero < 0.1) &
175 (((t_nonzero > 26.0) & (t_nonzero < 28.0))))
176 popt_5, pvoc_5 = scipy.optimize.curve_fit(
177 creep_rheology2, tau_nonzero[idxfit_]/N_nonzero[idxfit_],
178 shearstrainrate_nonzero[idxfit_])
179
180 idxfit_ = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
181 (shearstrainrate_nonzero < 0.1) &
182 (((t_nonzero > 31.0) & (t_nonzero < 33.0))))
183 popt_6, pvoc_6 = scipy.optimize.curve_fit(
184 creep_rheology2, tau_nonzero[idxfit_]/N_nonzero[idxfit_],
185 shearstrainrate_nonzero[idxfit_])
186
187 idxfit_ = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
188 (shearstrainrate_nonzero < 0.1) &
189 (((t_nonzero > 36.0) & (t_nonzero < 38.0))))
190 popt_7, pvoc_7 = scipy.optimize.curve_fit(
191 creep_rheology2, tau_nonzero[idxfit_]/N_nonzero[idxfit_],
192 shearstrainrate_nonzero[idxfit_])
193
194 idxfit_ = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
195 (shearstrainrate_nonzero < 0.1) &
196 (((t_nonzero > 41.0) & (t_nonzero < 43.0))))
197 popt_8, pvoc_8 = scipy.optimize.curve_fit(
198 creep_rheology2, tau_nonzero[idxfit_]/N_nonzero[idxfit_],
199 shearstrainrate_nonzero[idxfit_])
200
201 idxfit_ = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
202 (shearstrainrate_nonzero < 0.1) &
203 (((t_nonzero > 46.0) & (t_nonzero < 48.0))))
204 popt_9, pvoc_9 = scipy.optimize.curve_fit(
205 creep_rheology2, tau_nonzero[idxfit_]/N_nonzero[idxfit_],
206 shearstrainrate_nonzero[idxfit_])
207
208 idxfit_ = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
209 (shearstrainrate_nonzero < 0.1) &
210 (((t_nonzero > 51.0) & (t_nonzero < 53.0))))
211 popt_10, pvoc_10 = scipy.optimize.curve_fit(
212 creep_rheology2, tau_nonzero[idxfit_]/N_nonzero[idxfit_],
213 shearstrainrate_nonzero[idxfit_])
214
215 idxfit_ = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
216 (shearstrainrate_nonzero < 0.1) &
217 (((t_nonzero > 56.0) & (t_nonzero < 58.0))))
218 popt_11, pvoc_11 = scipy.optimize.curve_fit(
219 creep_rheology2, tau_nonzero[idxfit_]/N_nonzero[idxfit_],
220 shearstrainrate_nonzero[idxfit_])
221
222 n_list = [
223 popt_0[0],
224 popt_1[0],
225 popt_2[0],
226 popt_3[0],
227 popt_4[0],
228 popt_5[0],
229 popt_6[0],
230 popt_7[0],
231 popt_8[0],
232 popt_9[0],
233 popt_10[0],
234 popt_11[0]]
235
236 A_list = [
237 popt_0[1],
238 popt_1[1],
239 popt_2[1],
240 popt_3[1],
241 popt_4[1],
242 popt_5[1],
243 popt_6[1],
244 popt_7[1],
245 popt_8[1],
246 popt_9[1],
247 popt_10[1],
248 popt_11[1]]
249 cycle = numpy.arange(1, len(n_list)+1)
250
251 print(n_list)
252 print(A_list)
253 '''
254
255
256 friction = tau_nonzero/N_nonzero
257 x_min = numpy.floor(numpy.min(friction))
258 x_max = numpy.ceil(numpy.max(friction)) + 0.05
259 friction_fit =\
260 numpy.linspace(numpy.min(tau_nonzero[idxfit]/N_nonzero[idxfit]),
261 numpy.max(tau_nonzero[idxfit]/N_nonzero[idxfit]),
262 100)
263 strainrate_fit = A*friction_fit**n
264 '''
265 ### Consolidated state fit
266 #idxfit2 = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
267 #(shearstrainrate_nonzero < 0.1) &
268 #((t_nonzero > 0.0) & (t_nonzero < 2.5)))
269 idxfit2 = numpy.nonzero((shearstrain_nonzero < 0.1) &
270 (tau_nonzero/N_nonzero < 0.38))
271
272 popt2, pvoc2 = scipy.optimize.curve_fit(
273 creep_rheology2, tau_nonzero[idxfit2]/N_nonzero[idxfit2],
274 shearstrainrate_nonzero[idxfit2])
275 print '# Consolidated state'
276 print popt2
277 print pvoc2
278 n2 = popt2[0] # stress exponent
279 A2 = popt2[1] # prefactor
280
281 friction_fit2 =\
282 numpy.linspace(numpy.min(tau_nonzero[idxfit2]/N_nonzero[idxfit2]),
283 numpy.max(tau_nonzero[idxfit2]/N_nonzero[idxfit2]),
284 100)
285 strainrate_fit2 = A2*friction_fit2**n2
286 '''
287
288
289 ### Plotting
290 fig = plt.figure(figsize=(3.5,2.5))
291 ax1 = plt.subplot(111)
292 #ax1.semilogy(N/1000., v)
293 #ax1.semilogy(tau_nonzero/N_nonzero, v_nonzero, '+k')
294 #ax1.plot(tau/N, v, '.')
295 #CS = ax1.scatter(friction, v_nonzero, c=shearstrain_nonzero,
296 #linewidth=0)
297 CS = ax1.scatter(friction, shearstrainrate_nonzero,
298 #c=t_nonzero, linewidth=0.2,
299 c=shearstrain_nonzero, linewidth=0.05,
300 cmap=matplotlib.cm.get_cmap('afmhot'))
301 #CS = ax1.scatter(friction[idxfit], shearstrainrate_nonzero[idxfit],
302 #c=shearstrain_nonzero[idxfit], linewidth=0.2,
303 #cmap=matplotlib.cm.get_cmap('afmhot'))
304
305 ## plastic limit
306 x = [0.26, 0.26]
307 y = ax1.get_ylim()
308 limitcolor = '#333333'
309 ax1.plot(x, y, '--', linewidth=2, color=limitcolor)
310 ax1.text(x[0]+0.03, 2.0e-4*t_DEM_to_t_real,
311 'Yield strength', rotation=90, color=limitcolor,
312 bbox={'fc':'#ffffff', 'pad':3, 'alpha':0.7})
313
314
315 ## Fit
316 ax1.plot(friction_fit, strainrate_fit)
317 #ax1.plot(friction_fit2, strainrate_fit2)
318 ax1.annotate('$\\dot{{\\gamma}}$ = ' + '{:.2e} s$^{{-1}}$'.format(A) +
319 '$(\\tau/N)^{{{:.2f}}}$'.format(n),
320 xy = (friction_fit[40], strainrate_fit[40]),
321 xytext = (0.31+0.05, 2.0e-9*t_DEM_to_t_real),
322 arrowprops=dict(facecolor='blue', edgecolor='blue', shrink=0.1,
323 width=1, headwidth=4, frac=0.2),)
324 #xytext = (friction_fit[50]+0.15, strainrate_fit[50]-1.0e-5))#,
325 #arrowprops=dict(facecolor='black', shrink=0.05),)
326
327 ax1.set_yscale('log')
328 ax1.set_xlim([x_min, x_max])
329 #y_min = numpy.min(v_nonzero)*0.5
330 #y_max = numpy.max(v_nonzero)*2.0
331 y_min = numpy.min(shearstrainrate_nonzero)*0.5
332 y_max = numpy.max(shearstrainrate_nonzero)*2.0
333 ax1.set_ylim([y_min, y_max])
334 #ax1.set_xlim([friction_fit[0], friction_fit[-1]])
335 #ax1.set_ylim([strainrate_fit[0], strainrate_fit[-1]])
336
337 cb = plt.colorbar(CS)
338 cb.set_label('Shear strain $\\gamma$ [-]')
339
340 #ax1.set_xlabel('Effective normal stress [kPa]')
341 ax1.set_xlabel('Friction $\\tau/N$ [-]')
342 #ax1.set_ylabel('Shear velocity [m/s]')
343 ax1.set_ylabel('Shear strain rate $\\dot{\\gamma}$ [s$^{-1}$]')
344
345 plt.tight_layout()
346 filename = sid + '-rate-dependence.' + outformat
347 plt.savefig(filename)
348 print(filename)
349 plt.close()
350 #shutil.copyfile(filename, '/home/adc/articles/own/3/graphics/' + filename)
351
352 ## dilation vs. rate
353 '''
354 fig = plt.figure(figsize=(3.5,2.5))
355 ax1 = plt.subplot(111)
356 CS = ax1.scatter(friction, dilation[idx],
357 #tau_nonzero[idxfit2]/N_nonzero[idxfit2],
358 shearstrainrate_nonzero[idxfit2],
359 c=shearstrain_nonzero[idxfit2], linewidth=0.05,
360 #c=shearstrain_nonzero, linewidth=0.05,
361 cmap=matplotlib.cm.get_cmap('afmhot'))
362
363 ## plastic limit
364 x = [0.3, 0.3]
365 y = ax1.get_ylim()
366 limitcolor = '#333333'
367 ax1.plot(x, y, '--', linewidth=2, color=limitcolor)
368 ax1.text(x[0]+0.03, 2.0e-4,
369 'Yield strength', rotation=90, color=limitcolor,
370 bbox={'fc':'#ffffff', 'pad':3, 'alpha':0.7})
371 '''
372
373 ## Fit
374 ax1.plot(friction_fit, strainrate_fit)
375 #ax1.plot(friction_fit2, strainrate_fit2)
376 ax1.annotate('$\\dot{\\gamma} = (\\tau/N)^{6.4}$',
377 xy = (friction_fit[40], strainrate_fit[40]),
378 xytext = (0.32+0.05, 2.0e-9),
379 arrowprops=dict(facecolor='blue', edgecolor='blue', shrink=0.1,
380 width=1, headwidth=4, frac=0.2),)
381 #xytext = (friction_fit[50]+0.15, strainrate_fit[50]-1.0e-5))#,
382 #arrowprops=dict(facecolor='black', shrink=0.05),)
383 '''
384 '''
385
386 ax1.set_yscale('log')
387 ax1.set_xlim([x_min, x_max])
388 y_min = numpy.min(shearstrainrate_nonzero)*0.5
389 y_max = numpy.max(shearstrainrate_nonzero)*2.0
390 ax1.set_ylim([y_min, y_max])
391
392 cb = plt.colorbar(CS)
393 cb.set_label('Shear strain $\\gamma$ [-]')
394
395 #ax1.set_xlabel('Effective normal stress [kPa]')
396 ax1.set_xlabel('Friction $\\tau/N$ [-]')
397 #ax1.set_ylabel('Shear velocity [m/s]')
398 ax1.set_ylabel('Dilation $\\Delta h$ [m]')
399
400 plt.tight_layout()
401 filename = sid + '-rate-dependence-dilation.' + outformat
402 plt.savefig(filename)
403 plt.close()
404 shutil.copyfile(filename, '/home/adc/articles/own/3/graphics/' + filename)
405 print(filename)
406
407
408
409 fig = plt.figure(figsize=(3.5,2.5))
410 ax1 = plt.subplot(111)
411 #ax1.semilogy(N/1000., v)
412 #ax1.semilogy(tau_nonzero/N_nonzero, v_nonzero, '+k')
413 #ax1.plot(tau/N, v, '.')
414 #CS = ax1.scatter(friction, v_nonzero, c=shearstrain_nonzero,
415 #linewidth=0)
416 CS = ax1.plot(cycle, A_list, '-k')
417
418 ax1.set_xlabel('Forcing cycle [-]')
419 #ax1.set_ylabel('Shear velocity [m/s]')
420 ax1.set_ylabel('$A$ [s$^{-1}$]')
421
422 ax2 = ax1.twinx()
423 ax2color = 'blue'
424 lns = ax2.plot(cycle, n_list, ':', color=ax2color)
425 ax2.set_ylabel('$n$ [-]')
426 for tl in ax2.get_yticklabels():
427 tl.set_color(ax2color)
428
429 plt.tight_layout()
430 filename = sid + '-rate-dependence-A-n.' + outformat
431 plt.savefig(filename)
432 print(filename)
433 plt.close()
434 #shutil.copyfile(filename, '/home/adc/articles/own/3/graphics/' + filename)