thalfshear-darcy-strength-dilation-rate.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-strength-dilation-rate.py (17044B)
---
1 #!/usr/bin/env python
2 import matplotlib
3 matplotlib.use('Agg')
4 matplotlib.rcParams.update({'font.size': 18, 'font.family': 'serif'})
5 matplotlib.rc('text', usetex=True)
6 matplotlib.rcParams['text.latex.preamble']=[r"\usepackage{amsmath}"]
7 import shutil
8
9 import os
10 import sys
11 import numpy
12 import sphere
13 from permeabilitycalculator import *
14 import matplotlib.pyplot as plt
15
16 import seaborn as sns
17 #sns.set(style='ticks', palette='Set2')
18 sns.set(style='ticks', palette='colorblind')
19 #sns.set(style='ticks', palette='muted')
20 #sns.set(style='ticks', palette='pastel')
21 sns.despine() # remove right and top spines
22
23 pressures = True
24 zflow = False
25 contact_forces = False
26 smooth_friction = True
27 smooth_window = 100
28
29 #sigma0_list = numpy.array([1.0e3, 2.0e3, 4.0e3, 10.0e3, 20.0e3, 40.0e3])
30 sigma0_list = [20000.0, 80000.0]
31 #k_c_vals = [3.5e-13, 3.5e-15]
32 k_c = 3.5e-15
33 #k_c = 3.5e-13
34
35 # 5.0e-8 results present
36 #mu_f_vals = [1.797e-06, 1.204e-06, 5.0e-8, 1.797e-08]
37 #mu_f_vals = [1.797e-06, 1.204e-06, 3.594e-07, 1.797e-08]
38 mu_f_vals = [1.797e-06, 1.797e-07, 1.797e-08]
39 #mu_f_vals = [1.797e-06, 1.204e-06, 1.797e-08]
40 #mu_f_vals = [1.797e-06, 3.594e-07, 1.797e-08]
41 velfac = 1.0
42
43 # return a smoothed version of in. The returned array is smaller than the
44 # original input array
45 def smooth(x, window_len=10, window='hanning'):
46 """smooth the data using a window with requested size.
47
48 This method is based on the convolution of a scaled window with the signal.
49 The signal is prepared by introducing reflected copies of the signal
50 (with the window size) in both ends so that transient parts are minimized
51 in the begining and end part of the output signal.
52
53 input:
54 x: the input signal
55 window_len: the dimension of the smoothing window
56 window: the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'
57 flat window will produce a moving average smoothing.
58
59 output:
60 the smoothed signal
61
62 example:
63
64 import numpy as np
65 t = np.linspace(-2,2,0.1)
66 x = np.sin(t)+np.random.randn(len(t))*0.1
67 y = smooth(x)
68
69 see also:
70
71 numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve
72 scipy.signal.lfilter
73
74 TODO: the window parameter could be the window itself if an array instead of a string
75 """
76
77 if x.ndim != 1:
78 raise ValueError, "smooth only accepts 1 dimension arrays."
79
80 if x.size < window_len:
81 raise ValueError, "Input vector needs to be bigger than window size."
82
83 if window_len < 3:
84 return x
85
86 if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
87 raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"
88
89 s=numpy.r_[2*x[0]-x[window_len:1:-1], x, 2*x[-1]-x[-1:-window_len:-1]]
90 #print(len(s))
91
92 if window == 'flat': #moving average
93 w = numpy.ones(window_len,'d')
94 else:
95 w = getattr(numpy, window)(window_len)
96 y = numpy.convolve(w/w.sum(), s, mode='same')
97 return y[window_len-1:-window_len+1]
98
99 for sigma0 in sigma0_list:
100 shear_strain = [[], [], [], []]
101 friction = [[], [], [], []]
102 dilation = [[], [], [], []]
103 p_min = [[], [], [], []]
104 p_mean = [[], [], [], []]
105 p_max = [[], [], [], []]
106 f_n_mean = [[], [], [], []]
107 f_n_max = [[], [], [], []]
108 v_f_z_mean = [[], [], [], []]
109
110 fluid=True
111
112 # wet shear
113 for c in numpy.arange(0,len(mu_f_vals)):
114 #for c in numpy.arange(len(mu_f_vals)-1, -1, -1):
115 mu_f = mu_f_vals[c]
116
117 # halfshear-darcy-sigma0=20000.0-k_c=3.5e-13-mu=1.797e-06-velfac=1.0-shear
118 sid = 'halfshear-darcy-sigma0=' + str(sigma0) + '-k_c=' + str(k_c) + \
119 '-mu=' + str(mu_f) + '-velfac=' + str(velfac) + '-shear'
120 #sid = 'halfshear-sigma0=' + str(sigma0) + '-c_v=' + str(c_v) +\
121 #'-c_a=0.0-velfac=1.0-shear'
122 if os.path.isfile('../output/' + sid + '.status.dat'):
123
124 sim = sphere.sim(sid, fluid=fluid)
125 n = sim.status()
126 #n = 20
127 shear_strain[c] = numpy.zeros(n)
128 friction[c] = numpy.zeros_like(shear_strain[c])
129 dilation[c] = numpy.zeros_like(shear_strain[c])
130
131 '''
132 sim.readlast(verbose=False)
133 #sim.visualize('shear')
134 shear_strain[c] = sim.shear_strain
135 #shear_strain[c] = numpy.arange(sim.status()+1)
136 #friction[c] = sim.tau/1000.0#/sim.sigma_eff
137 friction[c] = sim.shearStress('effective')/sim.currentNormalStress('defined')
138 dilation[c] = sim.dilation
139 '''
140
141 # fluid pressures and particle forces
142 p_mean[c] = numpy.zeros_like(shear_strain[c])
143 p_min[c] = numpy.zeros_like(shear_strain[c])
144 p_max[c] = numpy.zeros_like(shear_strain[c])
145 f_n_mean[c] = numpy.zeros_like(shear_strain[c])
146 f_n_max[c] = numpy.zeros_like(shear_strain[c])
147
148 for i in numpy.arange(n):
149
150 sim.readstep(i, verbose=False)
151
152 shear_strain[c][i] = sim.shearStrain()
153 friction[c][i] = sim.shearStress('effective')/sim.currentNormalStress('defined')
154 dilation[c][i] = sim.w_x[0]
155
156 if pressures:
157 iz_top = int(sim.w_x[0]/(sim.L[2]/sim.num[2]))-1
158 p_mean[c][i] = numpy.mean(sim.p_f[:,:,0:iz_top])/1000
159 p_min[c][i] = numpy.min(sim.p_f[:,:,0:iz_top])/1000
160 p_max[c][i] = numpy.max(sim.p_f[:,:,0:iz_top])/1000
161
162 if contact_forces:
163 sim.findNormalForces()
164 f_n_mean[c][i] = numpy.mean(sim.f_n_magn)
165 f_n_max[c][i] = numpy.max(sim.f_n_magn)
166
167 if zflow:
168 v_f_z_mean[c] = numpy.zeros_like(shear_strain[c])
169 for i in numpy.arange(n):
170 v_f_z_mean[c][i] = numpy.mean(sim.v_f[:,:,:,2])
171
172 dilation[c] =\
173 (dilation[c] - dilation[c][0])/(numpy.mean(sim.radius)*2.0)
174
175 else:
176 print(sid + ' not found')
177
178 # produce VTK files
179 #for sid in sids:
180 #sim = sphere.sim(sid, fluid=True)
181 #sim.writeVTKall()
182
183
184 if zflow or pressures:
185 #fig = plt.figure(figsize=(8,10))
186 #fig = plt.figure(figsize=(3.74, 2*3.74))
187 fig = plt.figure(figsize=(2*3.74, 2*3.74))
188 else:
189 fig = plt.figure(figsize=(8,8)) # (w,h)
190 #fig = plt.figure(figsize=(8,12))
191 #fig = plt.figure(figsize=(8,16))
192
193 #plt.subplot(3,1,1)
194 #plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
195
196 for c in numpy.arange(0,len(mu_f_vals)):
197
198 if zflow or pressures:
199 ax1 = plt.subplot(3, len(mu_f_vals), 1+c)
200 ax2 = plt.subplot(3, len(mu_f_vals), 4+c, sharex=ax1)
201 ax3 = plt.subplot(3, len(mu_f_vals), 7+c, sharex=ax1)
202 else:
203 ax1 = plt.subplot(211)
204 ax2 = plt.subplot(212, sharex=ax1)
205 #ax3 = plt.subplot(413, sharex=ax1)
206 #ax4 = plt.subplot(414, sharex=ax1)
207 #alpha = 0.5
208 alpha = 1.0
209 #ax1.plot(shear_strain[0], friction[0], label='dry', linewidth=1, alpha=alpha)
210 #ax2.plot(shear_strain[0], dilation[0], label='dry', linewidth=1)
211 #ax4.plot(shear_strain[0], f_n_mean[0], '-', label='dry', color='blue')
212 #ax4.plot(shear_strain[0], f_n_max[0], '--', color='blue')
213
214 #color = ['b','g','r','c']
215 #color = ['g','r','c']
216 color = sns.color_palette()
217 #for c, mu_f in enumerate(mu_f_vals):
218 #for c in numpy.arange(len(mu_f_vals)-1, -1, -1):
219 mu_f = mu_f_vals[c]
220
221 if numpy.isclose(mu_f, 1.797e-6):
222 label = 'ref.\\ shear velocity'
223 #elif numpy.isclose(mu_f, 1.204e-6):
224 #label = 'ref.\\ shear velocity$\\times$0.67'
225 #elif numpy.isclose(mu_f, 1.797e-8):
226 #label = 'ref.\\ shear velocity$\\times$0.01'
227 else:
228 #label = '$\\mu_\\text{{f}}$ = {:.3e} Pa s'.format(mu_f)
229 label = 'shear velocity$\\times${:.2}'.format(mu_f/mu_f_vals[0])
230
231 # unsmoothed
232 ax1.plot(shear_strain[c][1:], friction[c][1:], \
233 label=label, linewidth=1,
234 alpha=0.3, color=color[c], clip_on=False)
235 #alpha=0.2, color='gray')
236 #alpha=alpha, color=color[c])
237
238 # smoothed
239 ax1.plot(shear_strain[c][1:-50], smooth(friction[c], smooth_window)[1:-50], \
240 label=label, linewidth=1,
241 alpha=alpha, color=color[c])
242
243
244 ax2.plot(shear_strain[c], dilation[c], \
245 label=label, linewidth=1,
246 color=color[c])
247
248 if zflow:
249 ax3.plot(shear_strain[c], v_f_z_mean[c],
250 label=label, linewidth=1)
251
252 if pressures:
253 ax3.plot(shear_strain[c], p_max[c], ':', color=color[c], alpha=0.5,
254 linewidth=0.5)
255
256 ax3.plot(shear_strain[c], p_mean[c], '-', color=color[c], \
257 label=label, linewidth=1)
258
259 ax3.plot(shear_strain[c], p_min[c], ':', color=color[c], alpha=0.5,
260 linewidth=0.5)
261
262
263 #ax3.fill_between(shear_strain[c], p_min[c], p_max[c],
264 # where=p_min[c]<=p_max[c], facecolor=color[c], edgecolor='None',
265 # interpolate=True, alpha=0.3)
266
267 #ax4.plot(shear_strain[c][1:], f_n_mean[c][1:], '-' + color[c],
268 #label='$c$ = %.2f' % (cvals[c-1]), linewidth=2)
269 #ax4.plot(shear_strain[c][1:], f_n_max[c][1:], '--' + color[c])
270 #label='$c$ = %.2f' % (cvals[c-1]), linewidth=2)
271
272 #ax4.set_xlabel('Shear strain $\\gamma$ [-]')
273 if zflow or pressures:
274 ax3.set_xlabel('Shear strain $\\gamma$ [-]')
275 else:
276 ax2.set_xlabel('Shear strain $\\gamma$ [-]')
277
278 if c == 0:
279 ax1.set_ylabel('Shear friction $\\tau/\\sigma_0$ [-]')
280 #ax1.set_ylabel('Shear stress $\\tau$ [kPa]')
281 ax2.set_ylabel('Dilation $\\Delta h/(2r)$ [-]')
282 if zflow:
283 ax3.set_ylabel('$\\boldsymbol{v}_\\text{f}^z h$ [ms$^{-1}$]')
284 if pressures:
285 ax3.set_ylabel('Fluid pressure $\\bar{p}_\\text{f}$ [kPa]')
286 #ax4.set_ylabel('Particle contact force $||\\boldsymbol{f}_\\text{p}||$ [N]')
287
288 #ax1.set_xlim([200,300])
289 #ax3.set_ylim([595,608])
290
291 plt.setp(ax1.get_xticklabels(), visible=False)
292 if zflow or pressures:
293 plt.setp(ax2.get_xticklabels(), visible=False)
294 #plt.setp(ax2.get_xticklabels(), visible=False)
295 #plt.setp(ax3.get_xticklabels(), visible=False)
296
297 '''
298 ax1.grid()
299 ax2.grid()
300 if zflow or pressures:
301 ax3.grid()
302 #ax4.grid()
303 '''
304
305 if c == 0: # left
306 # remove box at top and right
307 ax1.spines['top'].set_visible(False)
308 ax1.spines['bottom'].set_visible(False)
309 ax1.spines['right'].set_visible(False)
310 #ax1.spines['left'].set_visible(True)
311 # remove ticks at top and right
312 ax1.get_xaxis().set_ticks_position('none')
313 ax1.get_yaxis().set_ticks_position('none')
314 ax1.get_yaxis().tick_left()
315
316 # remove box at top and right
317 ax2.spines['top'].set_visible(False)
318 ax2.spines['right'].set_visible(False)
319 ax2.spines['bottom'].set_visible(False)
320 # remove ticks at top and right
321 ax2.get_xaxis().set_ticks_position('none')
322 ax2.get_yaxis().set_ticks_position('none')
323 ax2.get_yaxis().tick_left()
324
325 # remove box at top and right
326 ax3.spines['top'].set_visible(False)
327 ax3.spines['right'].set_visible(False)
328 # remove ticks at top and right
329 ax3.get_xaxis().set_ticks_position('none')
330 ax3.get_yaxis().set_ticks_position('none')
331 ax3.get_xaxis().tick_bottom()
332 ax3.get_yaxis().tick_left()
333
334 elif c == len(mu_f_vals)-1: # right
335 # remove box at top and right
336 ax1.spines['top'].set_visible(False)
337 ax1.spines['bottom'].set_visible(False)
338 ax1.spines['right'].set_visible(True)
339 ax1.spines['left'].set_visible(False)
340 # remove ticks at top and right
341 ax1.get_xaxis().set_ticks_position('none')
342 ax1.get_yaxis().set_ticks_position('none')
343 ax1.get_yaxis().tick_right()
344
345 # remove box at top and right
346 ax2.spines['top'].set_visible(False)
347 ax2.spines['right'].set_visible(True)
348 ax2.spines['bottom'].set_visible(False)
349 ax2.spines['left'].set_visible(False)
350 # remove ticks at top and right
351 ax2.get_xaxis().set_ticks_position('none')
352 ax2.get_yaxis().set_ticks_position('none')
353 #ax2.get_yaxis().tick_left()
354 ax2.get_yaxis().tick_right()
355
356 # remove box at top and right
357 ax3.spines['top'].set_visible(False)
358 ax3.spines['right'].set_visible(True)
359 ax3.spines['left'].set_visible(False)
360 # remove ticks at top and right
361 ax3.get_xaxis().set_ticks_position('none')
362 ax3.get_yaxis().set_ticks_position('none')
363 ax3.get_xaxis().tick_bottom()
364 #ax3.get_yaxis().tick_left()
365 ax3.get_yaxis().tick_right()
366
367 else: # middle
368 # remove box at top and right
369 ax1.spines['top'].set_visible(False)
370 ax1.spines['bottom'].set_visible(False)
371 ax1.spines['right'].set_visible(False)
372 ax1.spines['left'].set_visible(False)
373 # remove ticks at top and right
374 ax1.get_xaxis().set_ticks_position('none')
375 ax1.get_yaxis().set_ticks_position('none')
376 #ax1.get_yaxis().tick_left()
377 plt.setp(ax1.get_yticklabels(), visible=False)
378
379 # remove box at top and right
380 ax2.spines['top'].set_visible(False)
381 ax2.spines['right'].set_visible(False)
382 ax2.spines['bottom'].set_visible(False)
383 ax2.spines['left'].set_visible(False)
384 # remove ticks at top and right
385 ax2.get_xaxis().set_ticks_position('none')
386 ax2.get_yaxis().set_ticks_position('none')
387 #ax2.get_yaxis().tick_left()
388 plt.setp(ax2.get_yticklabels(), visible=False)
389
390 # remove box at top and right
391 ax3.spines['top'].set_visible(False)
392 ax3.spines['right'].set_visible(False)
393 ax3.spines['left'].set_visible(False)
394 # remove ticks at top and right
395 ax3.get_xaxis().set_ticks_position('none')
396 ax3.get_yaxis().set_ticks_position('none')
397 ax3.get_xaxis().tick_bottom()
398 #ax3.get_yaxis().tick_left()
399 plt.setp(ax3.get_yticklabels(), visible=False)
400
401
402 # vertical grid lines
403 ax1.get_xaxis().grid(True, linestyle=':', linewidth=0.5)
404 ax2.get_xaxis().grid(True, linestyle=':', linewidth=0.5)
405 ax3.get_xaxis().grid(True, linestyle=':', linewidth=0.5)
406
407
408 # horizontal grid lines
409 ax1.get_yaxis().grid(True, linestyle=':', linewidth=0.5)
410 ax2.get_yaxis().grid(True, linestyle=':', linewidth=0.5)
411 ax3.get_yaxis().grid(True, linestyle=':', linewidth=0.5)
412
413 ax1.set_title(label)
414 #ax1.legend(loc='best')
415 #legend_alpha=0.5
416 #ax1.legend(loc='upper right', prop={'size':18}, fancybox=True,
417 #framealpha=legend_alpha)
418 #ax2.legend(loc='lower right', prop={'size':18}, fancybox=True,
419 #framealpha=legend_alpha)
420 #if zflow or pressures:
421 #ax3.legend(loc='upper right', prop={'size':18}, fancybox=True,
422 #framealpha=legend_alpha)
423 #ax4.legend(loc='best', prop={'size':18}, fancybox=True,
424 #framealpha=legend_alpha)
425
426 #ax1.set_xlim([0.0, 0.09])
427 #ax2.set_xlim([0.0, 0.09])
428 #ax2.set_xlim([0.0, 0.2])
429
430 #ax1.set_ylim([-7, 45])
431 #ax1.set_xlim([0.0, 1.0])
432 ax1.set_xlim([0.0, 2.0])
433 ax1.set_ylim([0.0, 1.0])
434 ax2.set_ylim([0.0, 1.0])
435 ax3.set_ylim([-150., 150.])
436
437 #ax1.set_ylim([0.0, 1.0])
438 #if pressures:
439 #ax3.set_ylim([-1400, 900])
440 #ax3.set_ylim([-200, 200])
441 #ax3.set_xlim([0.0, 0.09])
442
443 #plt.tight_layout()
444 #plt.subplots_adjust(hspace=0.05)
445 plt.subplots_adjust(hspace=0.15)
446 #filename = 'shear-' + str(int(sigma0/1000.0)) + 'kPa-stress-dilation.pdf'
447 filename = 'halfshear-darcy-rate.pdf'
448 if sigma0 == 80000.0:
449 filename = 'halfshear-darcy-rate-N80.pdf'
450 #print(os.getcwd() + '/' + filename)
451 plt.savefig(filename)
452 shutil.copyfile(filename, '/home/adc/articles/own/2/graphics/' + filename)
453 plt.close()
454 print(filename)