thalfshear-darcy-combined.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-combined.py (12044B)
---
1 #!/usr/bin/env python
2 import matplotlib
3 matplotlib.use('Agg')
4 import shutil
5
6 import os
7 import sys
8 import numpy
9 import sphere
10 import matplotlib.pyplot as plt
11 import matplotlib.patches
12 import matplotlib.colors
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 #if len(sys.argv) > 1:
19 #sid = sys.argv[1]
20 #else:
21 #raise Exception('Missing sid')
22 sid = 'halfshear-darcy-sigma0=80000.0-k_c=3.5e-13-mu=1.04e-07-ss=10000.0-A=70000.0-f=0.2'
23 outformat = 'pdf'
24 fluid = True
25 #threshold = 100.0 # [N]
26 if len(sys.argv) > 3:
27 calculateforcechains = bool(int(sys.argv[3]))
28 else:
29 calculateforcechains = False
30 calculateforcechainhistory = False
31 legend_alpha=0.7
32 linewidth=0.5
33
34 rasterized=True # rasterize colored areas (pcolormesh and colorbar)
35 #izfactor = 4 # factor for vertical discretization in xvel
36 izfactor = 1 # factor for vertical discretization in poros
37
38
39 #if len(sys.argv) > 2:
40 # t_DEM_to_t_real = float(sys.argv[2])
41 #else:
42 # raise Exception('Missing t_DEM_to_t_real')
43
44
45 ###################
46 #### DATA READ ####
47 ###################
48
49 sim = sphere.sim(sid, fluid=fluid)
50 sim.readfirst(verbose=False)
51
52 #nsteps = 2
53 #nsteps = 10
54 #nsteps = 400
55 nsteps = sim.status()
56
57 t = numpy.empty(nsteps)
58
59 # Stress plot
60 sigma_def = numpy.empty_like(t)
61 sigma_eff = numpy.empty_like(t)
62 tau_def = numpy.empty_like(t)
63 tau_eff = numpy.empty_like(t)
64 p_f_bar = numpy.empty_like(t)
65 p_f_top = numpy.empty_like(t)
66
67 # shear velocity plot
68 v = numpy.empty_like(t)
69
70 # displacement and mean porosity plot
71 xdisp = numpy.empty_like(t)
72 phi_bar = numpy.empty_like(t)
73
74 # mean horizontal porosity plot
75 poros = numpy.empty((sim.num[2], nsteps))
76 xvel = numpy.zeros((sim.num[2]*izfactor, nsteps))
77 zpos_c = numpy.empty(sim.num[2]*izfactor)
78 dz = sim.L[2]/(sim.num[2]*izfactor)
79 for i in numpy.arange(sim.num[2]*izfactor):
80 zpos_c[i] = i*dz + 0.5*dz
81
82 # Contact statistics plot
83 n = numpy.empty(nsteps)
84 coordinationnumber = numpy.empty(nsteps)
85 nkept = numpy.empty(nsteps)
86
87
88 for i in numpy.arange(nsteps):
89
90 sim.readstep(i+1, verbose=False)
91
92 t[i] = sim.currentTime()
93
94 sigma_def[i] = sim.currentNormalStress('defined')
95 sigma_eff[i] = sim.currentNormalStress('effective')
96 tau_def[i] = sim.shearStress('defined')
97 tau_eff[i] = sim.shearStress('effective')
98
99 phi_bar[i] = numpy.mean(sim.phi[:,:,0:sim.wall0iz()])
100 p_f_bar[i] = numpy.mean(sim.p_f[:,:,0:sim.wall0iz()])
101 p_f_top[i] = sim.p_f[0,0,-1]
102
103 v[i] = sim.shearVelocity()
104 xdisp[i] = sim.shearDisplacement()
105
106 poros[:,i] = numpy.average(numpy.average(sim.phi, axis=0),axis=0)
107
108 # calculate mean values of xvel
109 #'''
110 dz = sim.L[2]/(sim.num[2]*izfactor)
111 for iz in numpy.arange(sim.num[2]*izfactor):
112 z_bot = iz*dz
113 z_top = (iz+1)*dz
114 I = numpy.nonzero((sim.x[:,2] >= z_bot) & (sim.x[:,2] < z_top))
115 if I[0].size > 0:
116 #xvel[iz,i] = numpy.mean(numpy.abs(sim.vel[I,0]))
117 xvel[iz,i] = numpy.mean(sim.vel[I,0])
118 #'''
119
120 if calculateforcechains:
121 if i > 0 and calculateforcechainhistory:
122 loaded_contacts_prev = numpy.copy(loaded_contacts)
123 pairs_prev = numpy.copy(sim.pairs)
124
125 loaded_contacts = sim.findLoadedContacts(
126 threshold = sim.currentNormalStress()*2.)
127 #sim.currentNormalStress()/1000.)
128 n[i] = numpy.size(loaded_contacts)
129
130 sim.findCoordinationNumber()
131 coordinationnumber[i] = sim.findMeanCoordinationNumber()
132
133 if calculateforcechainhistory:
134 nfound = 0
135 if i > 0:
136 for a in loaded_contacts[0]:
137 for b in loaded_contacts_prev[0]:
138 if (sim.pairs[:,a] == pairs_prev[:,b]).all():
139 nfound += 1;
140
141 nkept[i] = nfound
142 print nfound
143
144 print coordinationnumber[i]
145
146
147 if calculateforcechains:
148 numpy.savetxt(sid + '-fc.txt', (n, nkept, coordinationnumber))
149 else:
150 n, nkept, coordinationnumber = numpy.loadtxt(sid + '-fc.txt')
151
152 # Transform time from model time to real time [s]
153 t_DEM_to_t_real = sim.mu/1.787e-3
154 t = t/t_DEM_to_t_real
155
156 ## integrate velocities to displacement along x (xdispint)
157 # Taylor two term expansion
158 xdispint = numpy.zeros_like(t)
159 v_limit = 2.78e-3 # 1 m/hour (WIP)
160 dt = (t[1] - t[0])
161 dt2 = dt*2.
162 for i in numpy.arange(t.size):
163 if i > 0 and i < t.size-1:
164 acc = (numpy.min([v[i+1], v_limit]) - numpy.min([v[i-1], v_limit]))/dt2
165 xdispint[i] = xdispint[i-1] +\
166 numpy.min([v[i], v_limit])*dt + 0.5*acc*dt**2
167 elif i == t.size-1:
168 xdispint[i] = xdispint[i-1] + numpy.min([v[i], v_limit])*dt
169
170
171 ##################
172 #### PLOTTING ####
173 ##################
174 bbox_x = 0.03
175 bbox_y = 0.96
176 verticalalignment='top'
177 horizontalalignment='left'
178 fontweight='bold'
179 bbox={'facecolor':'white', 'alpha':1.0, 'pad':3}
180
181 # Time in days
182 t = t/(60.*60.*24.)
183
184 fig = plt.figure(figsize=[3.5,8])
185
186 ## ax1: N, tau, ax2: p_f
187 ax1 = plt.subplot(5, 1, 1)
188 lns0 = ax1.plot(t, sigma_def/1000., '-k', label="$N$",
189 linewidth=linewidth)
190 #lns1 = ax1.plot(t, sigma_eff/1000., '-k', label="$\\sigma'$")
191 #lns2 = ax1.plot(t, tau_def/1000., '-r', label="$\\tau$")
192 #ns2 = ax1.plot(t, tau_def/1000., '-r')
193 lns3 = ax1.plot(t, tau_eff/1000., '-r', label="$\\tau'$", linewidth=linewidth)
194
195 ax1.set_ylabel('Stress [kPa]')
196
197 ax2 = ax1.twinx()
198 ax2color = 'blue'
199 #lns4 = ax2.plot(t, p_f_top/1000.0 + 80.0, '-',
200 #color=ax2color,
201 #label='$p_\\text{f}^\\text{forcing}$')
202 lns5 = ax2.plot(t, p_f_bar/1000.0, ':',
203 color=ax2color,
204 label='$\\Delta\\bar{p}_\\text{f}$', linewidth=linewidth)
205 ax2.set_ylabel('Mean change in fluid pressure [kPa]')
206 ax2.yaxis.label.set_color(ax2color)
207 for tl in ax2.get_yticklabels():
208 tl.set_color(ax2color)
209 #ax2.legend(loc='upper right')
210 #lns = lns0+lns1+lns2+lns3+lns4+lns5
211 #lns = lns0+lns1+lns2+lns3+lns5
212 #lns = lns1+lns3+lns5
213 lns = lns0+lns3+lns5
214 labs = [l.get_label() for l in lns]
215 ax2.legend(lns, labs, loc='upper right', ncol=3,
216 fancybox=True, framealpha=legend_alpha)
217 ax1.set_ylim([-30, 200])
218 ax2.set_ylim([-115,125])
219
220 #ax1.text(bbox_x, bbox_y, 'A',
221 ax1.text(bbox_x, bbox_y, 'a',
222 horizontalalignment=horizontalalignment,
223 verticalalignment=verticalalignment,
224 fontweight=fontweight, bbox=bbox,
225 transform=ax1.transAxes)
226
227
228 ## ax3: v, ax4: unused
229 ax3 = plt.subplot(5, 1, 2, sharex=ax1)
230 #ax3.semilogy(t, v*t_DEM_to_t_real, 'k', linewidth=linewidth)
231 ax3.semilogy(t, v, 'k', linewidth=linewidth)
232 ax3.set_ylabel('Shear velocity [ms$^{-1}$]')
233 #ax3.set_ylabel('Shear velocity [ma$^{-1}$]')
234 #ax3.set_ylabel('Shear velocity [ma$^{-1}$]')
235 # shade stick periods
236 collection = matplotlib.collections.BrokenBarHCollection.span_where(
237 t, ymin=1.0e-11, ymax=1.0e4,
238 where=numpy.isclose(v, 0.0),
239 facecolor='black', alpha=0.2,
240 linewidth=0)
241 ax3.add_collection(collection)
242
243 #ax3.text(bbox_x, bbox_y, 'B',
244 ax3.text(bbox_x, bbox_y, 'b',
245 horizontalalignment=horizontalalignment,
246 verticalalignment=verticalalignment,
247 fontweight=fontweight, bbox=bbox,
248 transform=ax3.transAxes)
249
250
251 ## ax5: xdisp, ax6: mean(phi)
252 ax5 = plt.subplot(5, 1, 3, sharex=ax1)
253
254 #ax5.plot(t, xdisp, 'k', linewidth=linewidth)
255
256 # integrated displacement
257 #ax5.plot(t, xdispint, 'k', linewidth=linewidth)
258
259 # normalized displacement
260 #ax5.plot(t, xdisp/xdisp[-1], 'k', linewidth=linewidth)
261 ax5.plot(t, xdisp/xdisp[4500], 'k', linewidth=linewidth)
262
263 #ax5.set_ylabel('Shear displacement [m]')
264 ax5.set_ylabel('Normalized displacement [-]')
265 ax5.set_ylim([-0.02, 1.02])
266
267 ax6color='blue'
268 ax6 = ax5.twinx()
269 ax6.plot(t, phi_bar, color=ax6color, linewidth=linewidth)
270 ax6.set_ylabel('Mean porosity [-]')
271 ax6.yaxis.label.set_color(ax6color)
272 for tl in ax6.get_yticklabels():
273 tl.set_color(ax6color)
274
275 #ax6.text(bbox_x, bbox_y, 'C',
276 ax6.text(bbox_x, bbox_y, 'c',
277 horizontalalignment=horizontalalignment,
278 verticalalignment=verticalalignment,
279 fontweight=fontweight, bbox=bbox,
280 transform=ax6.transAxes)
281 #ax6.set_ylim([0.36, 0.39])
282 ax6.set_ylim([0.36, 0.40])
283
284
285 ## ax7: n_heavy, dn_heavy, ax8: z
286 ax7 = plt.subplot(5, 1, 4, sharex=ax1)
287 ax7.semilogy(t[:n.size], n, 'k', label='$n_\\text{heavy}$', linewidth=linewidth)
288 ax7.set_ylabel('Number of heavily loaded contacts [-]')
289 #ax7.semilogy(t, n - nkept, 'b', label='$\Delta n_\\text{heavy}$')
290 ax7.set_ylim([1.0e1, 2.0e4])
291
292 ax8 = ax7.twinx()
293 ax8color='green'
294 ax8.plot(t[:n.size], coordinationnumber, color=ax8color, linewidth=linewidth)
295 ax8.set_ylabel('Contacts per grain [-]')
296 ax8.yaxis.label.set_color(ax8color)
297 for tl in ax8.get_yticklabels():
298 tl.set_color(ax8color)
299 ax8.set_ylim([-0.2,9.8])
300
301 #ax7.text(bbox_x, bbox_y, 'D',
302 ax7.text(bbox_x, bbox_y, 'd',
303 horizontalalignment=horizontalalignment,
304 verticalalignment=verticalalignment,
305 fontweight=fontweight, bbox=bbox,
306 transform=ax7.transAxes)
307
308
309 ## ax9: porosity, ax10: unused
310 ax9 = plt.subplot(5, 1, 5, sharex=ax1)
311 poros_min = 0.375
312 poros_max = 0.45
313 poros[:,0] = poros[:,2] # remove erroneous porosity increase
314 cmap = matplotlib.cm.get_cmap('Blues_r')
315 #cmap = matplotlib.cm.get_cmap('afmhot')
316 #im9 = ax9.pcolormesh(t, zpos_c, poros,
317 #zpos_c = zpos_c[:-1]
318
319 xvel = xvel[:-1]
320 xvel[xvel < 0.0] = 0.0
321 im9 = ax9.pcolormesh(t, zpos_c, poros,
322 #im9 = ax9.pcolormesh(t, zpos_c, xvel,
323 cmap=cmap,
324 vmin=poros_min, vmax=poros_max,
325 #norm=matplotlib.colors.LogNorm(vmin=1.0e-8, vmax=xvel.max()),
326 shading='goraud',
327 rasterized=rasterized)
328 ax9.set_ylim([zpos_c[0], sim.w_x[0]])
329 ax9.set_ylabel('Vertical position [m]')
330
331 cbaxes = fig.add_axes([0.32, 0.1, 0.4, 0.01]) # x,y,w,h
332
333 #ax9.add_patch(matplotlib.patches.Rectangle(
334 #(3.0, 0.04), # x,y
335 #15., # dx
336 #.15, # dy
337 #fill=True,
338 #linewidth=1,
339 #facecolor='white'))
340 ax9.add_patch(matplotlib.patches.Rectangle(
341 (1.5, 0.04), # x,y
342 7., # dx
343 .15, # dy
344 fill=True,
345 linewidth=1,
346 facecolor='white',
347 alpha=legend_alpha))
348
349 cb9 = plt.colorbar(im9, cax=cbaxes,
350 #ticks=[poros_min, poros_min + 0.5*(poros_max-poros_min), poros_max],
351 #ticks=[xvel.min(), xvel.min() + 0.5*(xvel.max()-xvel.min()), xvel.max()],
352 orientation='horizontal',
353 extend='min',
354 cmap=cmap)
355 #cmap.set_under([8./255., 48./255., 107./255.]) # for poros
356 #cmap.set_under([1.0e-3, 1.0e-3, 1.0e-3]) # for xvel
357 from matplotlib import ticker
358 tick_locator = ticker.MaxNLocator(nbins=4)
359 cb9.locator = tick_locator
360 cb9.update_ticks()
361
362 cb9.set_label('Mean horizontal porosity [-]')
363 '''
364 ax9.text(0.5, 0.4, 'Mean horizontal porosity [-]\\\\',
365 horizontalalignment='center',
366 verticalalignment='center',
367 bbox={'facecolor':'white', 'alpha':1.0, 'pad':3})
368 '''
369 cb9.solids.set_rasterized(rasterized)
370
371 #ax9.text(bbox_x, bbox_y, 'E',
372 ax9.text(bbox_x, bbox_y, 'e',
373 horizontalalignment=horizontalalignment,
374 verticalalignment=verticalalignment,
375 fontweight=fontweight, bbox=bbox,
376 transform=ax9.transAxes)
377
378
379
380 plt.setp(ax1.get_xticklabels(), visible=False)
381 #plt.setp(ax2.get_xticklabels(), visible=False)
382 plt.setp(ax3.get_xticklabels(), visible=False)
383 #plt.setp(ax4.get_xticklabels(), visible=False)
384 plt.setp(ax5.get_xticklabels(), visible=False)
385 #plt.setp(ax6.get_xticklabels(), visible=False)
386 plt.setp(ax7.get_xticklabels(), visible=False)
387 #plt.setp(ax8.get_xticklabels(), visible=False)
388
389 ax1.set_xlim([0,9])
390 ax2.set_xlim([0,9])
391 ax3.set_xlim([0,9])
392 #ax4.set_xlim([0,9])
393 ax5.set_xlim([0,9])
394 ax6.set_xlim([0,9])
395 ax7.set_xlim([0,9])
396 ax8.set_xlim([0,9])
397 ax9.set_xlim([0,9])
398
399
400 ax9.set_xlabel('Time [d]')
401 fig.tight_layout()
402 plt.subplots_adjust(hspace=0.05)
403
404 filename = sid + '-combined.' + outformat
405 plt.savefig(filename)
406 plt.close()
407 shutil.copyfile(filename, '/home/adc/articles/own/3/graphics/' + filename)
408 print(filename)