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