tshear-results-internals.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
---
tshear-results-internals.py (14651B)
---
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 numpy
11 import sphere
12 from permeabilitycalculator import *
13 import matplotlib.pyplot as plt
14 from matplotlib.ticker import MaxNLocator
15
16 #steps = [5, 10, 100]
17 #steps = [5, 10]
18 steps = sys.argv[3:]
19 #nsteps_avg = 5 # no. of steps to average over
20 nsteps_avg = 100 # no. of steps to average over
21
22 sigma0 = float(sys.argv[1])
23 #c_grad_p = 1.0
24 c_grad_p = float(sys.argv[2])
25 c_phi = 1.0
26
27 #sid = 'shear-sigma0=' + str(sigma0) + '-c_phi=' + \
28 # str(c_phi) + '-c_grad_p=' + str(c_grad_p) + '-hi_mu-lo_visc'
29 sid = 'halfshear-sigma0=' + str(sigma0) + '-c=' + str(c_grad_p) + '-shear'
30 sim = sphere.sim(sid, fluid=True)
31 sim.readfirst(verbose=False)
32
33 # particle z positions
34 zpos_p = numpy.zeros((len(steps), sim.np))
35
36 # cell midpoint cell positions
37 zpos_c = numpy.zeros((len(steps), sim.num[2]))
38 dz = sim.L[2]/sim.num[2]
39 for i in numpy.arange(sim.num[2]):
40 zpos_c[:,i] = i*dz + 0.5*dz
41
42 # particle x displacements
43 xdisp = numpy.zeros((len(steps), sim.np))
44
45 # particle z velocity
46 v_z_p = numpy.zeros((len(steps), sim.np))
47
48 # fluid z velocity
49 v_z_f = numpy.zeros((len(steps), sim.num[0], sim.num[1], sim.num[2]))
50
51 # pressure - hydrostatic pressure
52 dev_p = numpy.zeros((len(steps), sim.num[2]))
53 p = numpy.zeros((len(steps), sim.num[2]))
54
55 # mean per-particle values
56 v_z_p_bar = numpy.zeros((len(steps), sim.num[2]))
57 v_z_f_bar = numpy.zeros((len(steps), sim.num[2]))
58
59 # particle-fluid force per particle
60 f_pf = numpy.zeros_like(xdisp)
61
62 # pressure - hydrostatic pressure
63 #dev_p = numpy.zeros((len(steps), sim.num[2]))
64
65 # mean porosity
66 phi_bar = numpy.zeros((len(steps), sim.num[2]))
67
68 # mean porosity change
69 dphi_bar = numpy.zeros((len(steps), sim.num[2]))
70
71 # mean per-particle values
72 xdisp_mean = numpy.zeros((len(steps), sim.num[2]))
73 f_pf_mean = numpy.zeros((len(steps), sim.num[2]))
74
75 shear_strain_start = numpy.zeros(len(steps))
76 shear_strain_end = numpy.zeros(len(steps))
77
78 s = 0
79 for step_str in steps:
80
81 step = int(step_str)
82
83 if os.path.isfile('../output/' + sid + '.status.dat'):
84
85 for substep in numpy.arange(nsteps_avg):
86
87 if step + substep > sim.status():
88 raise Exception(
89 'Simulation step %d not available (sim.status = %d).'
90 % (step + substep, sim.status()))
91
92 sim.readstep(step + substep, verbose=False)
93
94 zpos_p[s,:] += sim.x[:,2]/nsteps_avg
95
96 xdisp[s,:] += sim.xyzsum[:,0]/nsteps_avg
97 v_z_p[s,:] += sim.vel[:,2]/nsteps_avg
98
99 '''
100 for i in numpy.arange(sim.np):
101 f_pf[s,i] += \
102 sim.f_sum[i].dot(sim.f_sum[i])/nsteps_avg
103 '''
104 f_pf[s,:] += sim.f_sum[:,2]
105
106 dz = sim.L[2]/sim.num[2]
107 wall0_iz = int(sim.w_x[0]/dz)
108 for z in numpy.arange(0, wall0_iz+1):
109 dev_p[s,z] += \
110 (numpy.average(sim.p_f[:,:,z]) -
111 ((wall0_iz*dz - zpos_c[s][z] + 0.5*dz)
112 *sim.rho_f*numpy.abs(sim.g[2])\
113 + sim.p_f[0,0,-1])) \
114 /nsteps_avg
115
116 p[s,:] += numpy.average(numpy.average(sim.p_f[:,:,:], axis=0),\
117 axis=0)/nsteps_avg
118
119 v_z_f[s,:] += sim.v_f[:,:,:,2]/nsteps_avg
120
121 v_z_f_bar[s,:] += \
122 numpy.average(numpy.average(sim.v_f[:,:,:,2], axis=0), axis=0)\
123 /nsteps_avg
124
125 phi_bar[s,:] += \
126 numpy.average(numpy.average(sim.phi, axis=0), axis=0)\
127 /nsteps_avg
128
129 dphi_bar[s,:] += \
130 numpy.average(numpy.average(sim.dphi, axis=0), axis=0)\
131 /nsteps_avg/sim.time_dt
132
133
134 if substep == 0:
135 shear_strain_start[s] = sim.shearStrain()
136 else:
137 shear_strain_end[s] = sim.shearStrain()
138
139 # calculate mean values of xdisp and f_pf
140 for iz in numpy.arange(sim.num[2]):
141 z_bot = iz*dz
142 z_top = (iz+1)*dz
143 I = numpy.nonzero((zpos_p[s,:] >= z_bot) & (zpos_p[s,:] < z_top))
144 if len(I) > 0:
145 xdisp_mean[s,iz] = numpy.mean(xdisp[s,I])
146 v_z_p_bar[s,iz] = numpy.mean(v_z_p[s,I])
147 f_pf_mean[s,iz] = numpy.mean(f_pf[s,I])
148
149 else:
150 print(sid + ' not found')
151 s += 1
152
153 #fig = plt.figure(figsize=(8,4*(len(steps))+1))
154 #fig = plt.figure(figsize=(8,5*(len(steps))+1))
155 fig = plt.figure(figsize=(16,5*(len(steps))+1))
156
157 def color(c):
158 return 'black'
159 if c == 1.0:
160 return 'green'
161 elif c == 0.1:
162 return 'red'
163 elif c == 0.01:
164 return 'cyan'
165 else:
166 return 'blue'
167
168 ax = []
169 for s in numpy.arange(len(steps)):
170
171 #strain_str = 'Shear strain\n $\\gamma = %.3f$' % (shear_strain[s])
172
173 if s == 0:
174 strain_str = 'Dilating state\n$\\gamma = %.2f$ to $%.2f$\n$c = %.2f$' % \
175 (shear_strain_start[s], shear_strain_end[s], c_grad_p)
176 else:
177 strain_str = 'Steady state\n$\\gamma = %.2f$ to $%.2f$\n$c = %.2f$' % \
178 (shear_strain_start[s], shear_strain_end[s], c_grad_p)
179
180 n = 7
181 if s == 0:
182 ax.append(plt.subplot(len(steps), n-1, s*(n-1)+1)) # 0: xdisp
183 ax.append(plt.subplot(len(steps), n-1, s*(n-1)+2, sharey=ax[s*n+0])) # 1: phi
184 ax.append(ax[s*n+1].twiny()) # 2: dphi/dt
185 ax.append(plt.subplot(len(steps), n-1, s*(n-1)+3, sharey=ax[s*n+0])) # 3: v_z^p
186 ax.append(plt.subplot(len(steps), n-1, s*(n-1)+4, sharey=ax[s*n+0])) # 4: p_f
187 ax.append(plt.subplot(len(steps), n-1, s*(n-1)+5, sharey=ax[s*n+0])) # 5: f_pf_z
188 ax.append(plt.subplot(len(steps), n-1, s*(n-1)+6, sharey=ax[s*n+0])) # 6: v_z^f
189 else:
190 ax.append(plt.subplot(len(steps), n-1, s*(n-1)+1, sharex=ax[0])) # 0: xdisp
191 ax.append(plt.subplot(len(steps), n-1, s*(n-1)+2, sharey=ax[s*n+0],
192 sharex=ax[1])) # 1: phi
193 ax.append(ax[s*n+1].twiny()) # 2: dphi/dt
194 ax.append(plt.subplot(len(steps), n-1, s*(n-1)+3, sharey=ax[s*n+0],
195 sharex=ax[3])) # 3: v_z^p
196 ax.append(plt.subplot(len(steps), n-1, s*(n-1)+4, sharey=ax[s*n+0],
197 sharex=ax[4])) # 4: p_f
198 ax.append(plt.subplot(len(steps), n-1, s*(n-1)+5, sharey=ax[s*n+0],
199 sharex=ax[5])) # 5: f_pf_z
200 ax.append(plt.subplot(len(steps), n-1, s*(n-1)+6, sharey=ax[s*n+0],
201 sharex=ax[6])) # 6: v_z^f
202
203 ax[s*n+0].plot(xdisp[s], zpos_p[s], ',', color = '#888888')
204 ax[s*n+0].plot(xdisp_mean[s], zpos_c[s], color=color(c_grad_p))
205
206 #ax[s*4+2].plot(dev_p[s]/1000.0, zpos_c[s], 'k')
207 #ax[s*4+2].plot(phi_bar[s,1:], zpos_c[s,1:], '-k', linewidth=3)
208 ax[s*n+1].plot(phi_bar[s,1:], zpos_c[s,1:], '-', color=color(c_grad_p))
209
210 #phicolor = '#888888'
211 #ax[s*4+3].plot(phi_bar[s], zpos_c[s], '-', color = phicolor)
212 #for tl in ax[s*4+3].get_xticklabels():
213 #tl.set_color(phicolor)
214 ax[s*n+2].plot(dphi_bar[s,1:], zpos_c[s,1:], '--', color=color(c_grad_p))
215 #ax[s*4+3].plot(dphi_bar[s,1:], zpos_c[s,1:], '-k', linewidth=3)
216 #ax[s*4+3].plot(dphi_bar[s,1:], zpos_c[s,1:], '-w', linewidth=2)
217
218 ax[s*n+3].plot(v_z_p[s]*100.0, zpos_p[s], ',', alpha=0.5,
219 color='#888888')
220 #color=color(c_grad_p))
221 ax[s*n+3].plot(v_z_p_bar[s]*100.0, zpos_c[s], color=color(c_grad_p))
222 #ax[s*n+0].plot([0.0,0.0], [0.0, sim.L[2]], '--', color='k')
223
224 # hydrostatic pressure distribution
225 #ax[s*n+4].plot(dev_p[s]/1000.0, zpos_c[s], color=color(c_grad_p))
226 ax[s*n+4].plot(p[s]/1000.0, zpos_c[s], color=color(c_grad_p))
227 dz = sim.L[2]/sim.num[2]
228 wall0_iz = int(sim.w_x[0]/dz)
229 y_top = wall0_iz*dz + 0.5*dz
230 x_top = sim.p_f[0,0,-1]
231 y_bot = 0.0
232 x_bot = x_top + (wall0_iz*dz - zpos_c[s][0] + 0.5*dz)*sim.rho_f*numpy.abs(sim.g[2])
233 ax[s*n+4].plot([x_top/1000.0, x_bot/1000.0], [y_top, y_bot], '--', color='k')
234 #ax[s*n+1].set_title(strain_str)
235 #ax[s*n+1].set_title(' ')
236
237 # remove particles with 0.0 pressure force
238 I = numpy.nonzero(numpy.abs(f_pf[s]) > .01)
239 f_pf_nonzero = f_pf[s][I]
240 zpos_p_nonzero = zpos_p[s][I]
241 I = numpy.nonzero(numpy.abs(f_pf_mean[s]) > .01)
242 f_pf_mean_nonzero = f_pf_mean[s][I]
243 zpos_c_nonzero = zpos_c[s][I]
244
245 ax[s*n+5].plot(f_pf_nonzero, zpos_p_nonzero, ',', alpha=0.5,
246 color='#888888')
247 #color=color(c_grad_p))
248 #ax[s*4+1].plot(f_pf_mean[s][1:-2], zpos_c[s][1:-2], color = 'k')
249 ax[s*n+5].plot(f_pf_mean_nonzero, zpos_c_nonzero, color=color(c_grad_p))
250 #ax[s*4+1].plot([0.0, 0.0], [0.0, sim.L[2]], '--', color='k')
251
252 ax[s*n+6].plot(v_z_f_bar[s]*100.0, zpos_c[s], color=color(c_grad_p))
253 #ax[s*n+2].plot([0.0,0.0], [0.0, sim.L[2]], '--', color='k')
254
255
256 #plt.plot(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
257 #plt.semilogx(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
258 #plt.semilogy(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
259 #plt.loglog(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
260
261 #for x in numpy.arange(sim.num[0]):
262 #for y in numpy.arange(sim.num[1]):
263 #ax[s*n+2].plot(v_z_f[s,x,y,:], zpos_c[s], ',', color = '#888888')
264
265
266 #phicolor = '#888888'
267 #ax[s*n+3].plot(phi_bar[s], zpos_c[s], '-', color = phicolor)
268 #for tl in ax[s*n+3].get_xticklabels():
269 #tl.set_color(phicolor)
270
271 max_z = numpy.max(zpos_p)
272 ax[s*n+0].set_ylim([0, max_z])
273 #ax[s*n+1].set_ylim([0, max_z])
274 #ax[s*n+2].set_ylim([0, max_z])
275
276 #ax[s*n+0].set_xlim([-0.01,0.01])
277 #ax[s*n+0].set_xlim([-0.005,0.005])
278 #ax[s*n+0].set_xlim([-0.25,0.75])
279 #ax[s*n+4].set_xlim([595,625]) # p_f
280 #ax[s*n+2].set_xlim([-0.0005,0.0005])
281 #ax[s*n+2].set_xlim([-0.08,0.08])
282
283 #ax[s*4+1].set_xlim([0.15, 0.46]) # f_pf
284 #ax[s*n+1].set_xlim([0.235, 0.409]) # f_pf
285
286 ax[s*n+1].set_xlim([0.33, 0.6]) # phi
287 ax[s*n+2].set_xlim([-0.09, 0.035]) # dphi/dt
288 ax[s*n+3].set_xlim([-1.50, 1.50]) # v_z_p
289 ax[s*n+5].set_xlim([5.0, 8.0]) # f_z_pf
290
291
292 #plt.plot(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
293 #plt.semilogx(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
294 #plt.semilogy(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
295 #plt.loglog(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
296
297 ax[s*n+0].set_ylabel('Vertical position $z$ [m]')
298 ax[s*n+0].set_xlabel('$\\boldsymbol{x}^x_\\text{p}$ [m]')
299 ax[s*n+1].set_xlabel('$\\bar{\\phi}$ [-] (solid)')
300 ax[s*n+2].set_xlabel('$\\delta \\bar{\\phi}/\\delta t$ [-] (dashed)')
301 ax[s*n+3].set_xlabel('$\\boldsymbol{v}^z_\\text{p}$ [cms$^{-1}$]')
302 ax[s*n+4].set_xlabel('$\\bar{p_\\text{f}}$ [kPa]')
303 #ax[s*n+4].set_xlabel('$\\bar{p_\\text{f}} - p_\\text{hyd}$ [kPa]')
304 ax[s*n+5].set_xlabel('$\\boldsymbol{f}^z_\\text{pf}$ [N]')
305 ax[s*n+6].set_xlabel('$\\bar{\\boldsymbol{v}}^z_\\text{f}$ [cms$^{-1}$]')
306
307 # align x labels
308 labely = -0.3
309 ax[s*n+0].xaxis.set_label_coords(0.5, labely)
310 ax[s*n+1].xaxis.set_label_coords(0.5, labely)
311 #ax[s*n+2].xaxis.set_label_coords(0.5, labely)
312 ax[s*n+3].xaxis.set_label_coords(0.5, labely)
313 ax[s*n+4].xaxis.set_label_coords(0.5, labely)
314 ax[s*n+5].xaxis.set_label_coords(0.5, labely)
315 ax[s*n+6].xaxis.set_label_coords(0.5, labely)
316
317 plt.setp(ax[s*n+1].get_yticklabels(), visible=False)
318 plt.setp(ax[s*n+2].get_yticklabels(), visible=False)
319 plt.setp(ax[s*n+3].get_yticklabels(), visible=False)
320 plt.setp(ax[s*n+4].get_yticklabels(), visible=False)
321 plt.setp(ax[s*n+5].get_yticklabels(), visible=False)
322 plt.setp(ax[s*n+6].get_yticklabels(), visible=False)
323
324 #nbins = 4
325 #ax[s*n+0].get_xaxis().set_major_locator(MaxNLocator(nbins=nbins))
326 #ax[s*n+1].get_xaxis().set_major_locator(MaxNLocator(nbins=nbins))
327 #ax[s*n+2].get_xaxis().set_major_locator(MaxNLocator(nbins=nbins))
328
329 plt.setp(ax[s*n+0].xaxis.get_majorticklabels(), rotation=90)
330 plt.setp(ax[s*n+1].xaxis.get_majorticklabels(), rotation=90)
331 plt.setp(ax[s*n+2].xaxis.get_majorticklabels(), rotation=90)
332 plt.setp(ax[s*n+3].xaxis.get_majorticklabels(), rotation=90)
333 plt.setp(ax[s*n+4].xaxis.get_majorticklabels(), rotation=90)
334 plt.setp(ax[s*n+5].xaxis.get_majorticklabels(), rotation=90)
335 plt.setp(ax[s*n+6].xaxis.get_majorticklabels(), rotation=90)
336
337 #if s == 0:
338 #y = 0.95
339 #if s == 1:
340 #y = 0.55
341
342 #plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
343 #ax[s*n+0].ticklabel_format(style='sci', axis='x', scilimits=(-3,3))
344 #ax[s*n+1].ticklabel_format(style='sci', axis='x', scilimits=(-3,3))
345 #ax[s*n+2].ticklabel_format(style='sci', axis='x', scilimits=(-3,3))
346
347 #fig.text(0.1, y, strain_str, horizontalalignment='left', fontsize=22)
348 #ax[s*4+0].annotate(strain_str, xytext=(0,1.1), textcoords='figure fraction',
349 #horizontalalignment='left', fontsize=22)
350 #plt.text(0.05, 1.06, strain_str, horizontalalignment='left', fontsize=22,
351 #transform=ax[s*n+0].transAxes)
352 #ax[s*4+0].set_title(strain_str)
353
354 #ax[s*n+0].set_title('a')
355 #ax[s*n+1].set_title('b')
356 #ax[s*n+3].set_title('c')
357 #ax[s*n+4].set_title('d')
358 #ax[s*n+5].set_title('e')
359 #ax[s*n+6].set_title('f')
360
361 ax[s*n+0].grid()
362 ax[s*n+1].grid()
363 #ax[s*n+2].grid()
364 ax[s*n+3].grid()
365 ax[s*n+4].grid()
366 ax[s*n+5].grid()
367 ax[s*n+6].grid()
368 #ax1.legend(loc='lower right', prop={'size':18})
369 #ax2.legend(loc='lower right', prop={'size':18})
370
371 #strain_str = 'Shear strain $\\gamma = %.3f$' % (shear_strain[s])
372 #fig.text(0.1, y, strain_str, horizontalalignment='left', fontsize=22)
373 #ax[s*4+0].annotate(strain_str, xytext=(0,1.1), textcoords='figure fraction',
374 #horizontalalignment='left', fontsize=22)
375 plt.text(-0.38, 1.10, strain_str, horizontalalignment='left', fontsize=22,
376 transform=ax[s*n+0].transAxes)
377
378 #plt.title(' ')
379 plt.MaxNLocator(nbins=4)
380 #ax1.legend(loc='lower right', prop={'size':18})
381 #ax2.legend(loc='lower right', prop={'size':18})
382
383 #plt.title(' ')
384 #plt.MaxNLocator(nbins=4)
385 #plt.subplots_adjust(wspace = .05)
386 #plt.subplots_adjust(hspace = 1.05)
387 plt.tight_layout()
388 #plt.MaxNLocator(nbins=4)
389
390 filename = 'shear-' + str(int(sigma0/1000.0)) + 'kPa-internals.pdf'
391 plt.savefig(filename)
392 shutil.copyfile(filename, '/home/adc/articles/own/2-org/' + filename)
393 print(filename)