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