visualize.py - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
HTML git clone git://src.adamsgaard.dk/sphere
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
visualize.py (38553B)
---
1 import subprocess
2 import pickle as pl
3 import numpy
4 from .common import FontProperties, legend_alpha, matplotlib, plt, py_mpl
5
6
7 class SimVisualize:
8 'Time-series visualization of simulation output.'
9
10 def visualize(self, method='energy', savefig=True, outformat='png',
11 figsize=False, pickle=False, xlim=False, firststep=0,
12 f_min=None, f_max=None, cmap=None, smoothing=0,
13 smoothing_window='hanning'):
14 '''
15 Visualize output from the simulation, where the temporal progress is
16 of interest. The output will be saved in the current folder with a name
17 combining the simulation id of the simulation, and the visualization
18 method.
19
20 :param method: The type of plot to render. Possible values are 'energy',
21 'walls', 'triaxial', 'inertia', 'mean-fluid-pressure',
22 'fluid-pressure', 'shear', 'shear-displacement', 'porosity',
23 'rate-dependence', 'contacts'
24 :type method: str
25 :param savefig: Save the image instead of showing it on screen
26 :type savefig: bool
27 :param outformat: The output format of the plot data. This can be an
28 image format, or in text ('txt').
29 :param figsize: Specify output figure size in inches
30 :type figsize: array
31 :param pickle: Save all figure content as a Python pickle file. It can
32 be opened later using `fig=pickle.load(open('file.pickle','rb'))`.
33 :type pickle: bool
34 :param xlim: Set custom limits to the x axis. If not specified, the x
35 range will correspond to the entire data interval.
36 :type xlim: array
37 :param firststep: The first output file step to read (default: 0)
38 :type firststep: int
39 :param cmap: Choose custom color map, e.g.
40 `cmap=matplotlib.cm.get_cmap('afmhot')`
41 :type cmap: matplotlib.colors.LinearSegmentedColormap
42 :param smoothing: Apply smoothing across a number of output files to the
43 `method='shear'` plot. A value of less than 3 means that no
44 smoothing occurs.
45 :type smoothing: int
46 :param smoothing_window: Type of smoothing to use when `smoothing >= 3`.
47 Valid values are 'flat', 'hanning' (default), 'hamming', 'bartlett',
48 and 'blackman'.
49 :type smoothing_window: str
50 '''
51
52 lastfile = self.status()
53 from .core import sim
54 sb = sim(sid=self.sid, np=self.np, nw=self.nw, fluid=self.fluid)
55
56 if not py_mpl:
57 print('Error: matplotlib module not found (visualize).')
58 return
59
60 ### Plotting
61 if outformat != 'txt':
62 if figsize:
63 fig = plt.figure(figsize=figsize)
64 else:
65 fig = plt.figure(figsize=(8, 8))
66
67 if method == 'energy':
68 if figsize:
69 fig = plt.figure(figsize=figsize)
70 else:
71 fig = plt.figure(figsize=(20, 8))
72
73 # Allocate arrays
74 t = numpy.zeros(lastfile-firststep + 1)
75 Epot = numpy.zeros_like(t)
76 Ekin = numpy.zeros_like(t)
77 Erot = numpy.zeros_like(t)
78 Es = numpy.zeros_like(t)
79 Ev = numpy.zeros_like(t)
80 Es_dot = numpy.zeros_like(t)
81 Ev_dot = numpy.zeros_like(t)
82 Ebondpot = numpy.zeros_like(t)
83 Esum = numpy.zeros_like(t)
84
85 # Read energy values from simulation binaries
86 for i in numpy.arange(firststep, lastfile+1):
87 sb.readstep(i, verbose=False)
88
89 Epot[i] = sb.energy("pot")
90 Ekin[i] = sb.energy("kin")
91 Erot[i] = sb.energy("rot")
92 Es[i] = sb.energy("shear")
93 Ev[i] = sb.energy("visc_n")
94 Es_dot[i] = sb.energy("shearrate")
95 Ev_dot[i] = sb.energy("visc_n_rate")
96 Ebondpot[i] = sb.energy("bondpot")
97 Esum[i] = Epot[i] + Ekin[i] + Erot[i] + Es[i] + Ev[i] +\
98 Ebondpot[i]
99 t[i] = sb.currentTime()
100
101
102 if outformat != 'txt':
103 # Potential energy
104 ax1 = plt.subplot2grid((2, 5), (0, 0))
105 ax1.set_xlabel('Time [s]')
106 ax1.set_ylabel('Total potential energy [J]')
107 ax1.plot(t, Epot, '+-')
108 ax1.grid()
109
110 # Kinetic energy
111 ax2 = plt.subplot2grid((2, 5), (0, 1))
112 ax2.set_xlabel('Time [s]')
113 ax2.set_ylabel('Total kinetic energy [J]')
114 ax2.plot(t, Ekin, '+-')
115 ax2.grid()
116
117 # Rotational energy
118 ax3 = plt.subplot2grid((2, 5), (0, 2))
119 ax3.set_xlabel('Time [s]')
120 ax3.set_ylabel('Total rotational energy [J]')
121 ax3.plot(t, Erot, '+-')
122 ax3.grid()
123
124 # Bond energy
125 ax4 = plt.subplot2grid((2, 5), (0, 3))
126 ax4.set_xlabel('Time [s]')
127 ax4.set_ylabel('Bond energy [J]')
128 ax4.plot(t, Ebondpot, '+-')
129 ax4.grid()
130
131 # Total energy
132 ax5 = plt.subplot2grid((2, 5), (0, 4))
133 ax5.set_xlabel('Time [s]')
134 ax5.set_ylabel('Total energy [J]')
135 ax5.plot(t, Esum, '+-')
136 ax5.grid()
137
138 # Shear energy rate
139 ax6 = plt.subplot2grid((2, 5), (1, 0))
140 ax6.set_xlabel('Time [s]')
141 ax6.set_ylabel('Frictional dissipation rate [W]')
142 ax6.plot(t, Es_dot, '+-')
143 ax6.grid()
144
145 # Shear energy
146 ax7 = plt.subplot2grid((2, 5), (1, 1))
147 ax7.set_xlabel('Time [s]')
148 ax7.set_ylabel('Total frictional dissipation [J]')
149 ax7.plot(t, Es, '+-')
150 ax7.grid()
151
152 # Visc_n energy rate
153 ax8 = plt.subplot2grid((2, 5), (1, 2))
154 ax8.set_xlabel('Time [s]')
155 ax8.set_ylabel('Viscous dissipation rate [W]')
156 ax8.plot(t, Ev_dot, '+-')
157 ax8.grid()
158
159 # Visc_n energy
160 ax9 = plt.subplot2grid((2, 5), (1, 3))
161 ax9.set_xlabel('Time [s]')
162 ax9.set_ylabel('Total viscous dissipation [J]')
163 ax9.plot(t, Ev, '+-')
164 ax9.grid()
165
166 # Combined view
167 ax10 = plt.subplot2grid((2, 5), (1, 4))
168 ax10.set_xlabel('Time [s]')
169 ax10.set_ylabel('Energy [J]')
170 ax10.plot(t, Epot, '+-g')
171 ax10.plot(t, Ekin, '+-b')
172 ax10.plot(t, Erot, '+-r')
173 ax10.legend(('$\sum E_{pot}$', '$\sum E_{kin}$',
174 '$\sum E_{rot}$'), 'upper right', shadow=True)
175 ax10.grid()
176
177 if xlim:
178 ax1.set_xlim(xlim)
179 ax2.set_xlim(xlim)
180 ax3.set_xlim(xlim)
181 ax4.set_xlim(xlim)
182 ax5.set_xlim(xlim)
183 ax6.set_xlim(xlim)
184 ax7.set_xlim(xlim)
185 ax8.set_xlim(xlim)
186 ax9.set_xlim(xlim)
187 ax10.set_xlim(xlim)
188
189 fig.tight_layout()
190
191 elif method == 'walls':
192
193 # Read energy values from simulation binaries
194 for i in numpy.arange(firststep, lastfile+1):
195 sb.readstep(i, verbose=False)
196
197 # Allocate arrays on first run
198 if i == firststep:
199 wforce = numpy.zeros((lastfile+1)*sb.nw,\
200 dtype=numpy.float64).reshape((lastfile+1), sb.nw)
201 wvel = numpy.zeros((lastfile+1)*sb.nw,\
202 dtype=numpy.float64).reshape((lastfile+1), sb.nw)
203 wpos = numpy.zeros((lastfile+1)*sb.nw,\
204 dtype=numpy.float64).reshape((lastfile+1), sb.nw)
205 wsigma0 = numpy.zeros((lastfile+1)*sb.nw,\
206 dtype=numpy.float64).reshape((lastfile+1), sb.nw)
207 maxpos = numpy.zeros((lastfile+1), dtype=numpy.float64)
208 logstress = numpy.zeros((lastfile+1), dtype=numpy.float64)
209 voidratio = numpy.zeros((lastfile+1), dtype=numpy.float64)
210
211 wforce[i] = sb.w_force[0]
212 wvel[i] = sb.w_vel[0]
213 wpos[i] = sb.w_x[0]
214 wsigma0[i] = sb.w_sigma0[0]
215 maxpos[i] = numpy.max(sb.x[:, 2]+sb.radius)
216 logstress[i] = numpy.log((sb.w_force[0]/(sb.L[0]*sb.L[1]))/1000.0)
217 voidratio[i] = sb.voidRatio()
218
219 t = numpy.linspace(0.0, sb.time_current, lastfile+1)
220
221 # Plotting
222 if outformat != 'txt':
223 # linear plot of time vs. wall position
224 ax1 = plt.subplot2grid((2, 2), (0, 0))
225 ax1.set_xlabel('Time [s]')
226 ax1.set_ylabel('Position [m]')
227 ax1.plot(t, wpos, '+-', label="upper wall")
228 ax1.plot(t, maxpos, '+-', label="heighest particle")
229 ax1.legend()
230 ax1.grid()
231
232 #ax2 = plt.subplot2grid((2, 2), (1, 0))
233 #ax2.set_xlabel('Time [s]')
234 #ax2.set_ylabel('Force [N]')
235 #ax2.plot(t, wforce, '+-')
236
237 # semilog plot of log stress vs. void ratio
238 ax2 = plt.subplot2grid((2, 2), (1, 0))
239 ax2.set_xlabel('log deviatoric stress [kPa]')
240 ax2.set_ylabel('Void ratio [-]')
241 ax2.plot(logstress, voidratio, '+-')
242 ax2.grid()
243
244 # linear plot of time vs. wall velocity
245 ax3 = plt.subplot2grid((2, 2), (0, 1))
246 ax3.set_xlabel('Time [s]')
247 ax3.set_ylabel('Velocity [m/s]')
248 ax3.plot(t, wvel, '+-')
249 ax3.grid()
250
251 # linear plot of time vs. deviatoric stress
252 ax4 = plt.subplot2grid((2, 2), (1, 1))
253 ax4.set_xlabel('Time [s]')
254 ax4.set_ylabel('Deviatoric stress [Pa]')
255 ax4.plot(t, wsigma0, '+-', label="$\sigma_0$")
256 ax4.plot(t, wforce/(sb.L[0]*sb.L[1]), '+-', label="$\sigma'$")
257 ax4.legend(loc=4)
258 ax4.grid()
259
260 if xlim:
261 ax1.set_xlim(xlim)
262 ax2.set_xlim(xlim)
263 ax3.set_xlim(xlim)
264 ax4.set_xlim(xlim)
265
266 elif method == 'triaxial':
267
268 # Read energy values from simulation binaries
269 for i in numpy.arange(firststep, lastfile+1):
270 sb.readstep(i, verbose=False)
271
272 vol = (sb.w_x[0]-sb.origo[2]) * (sb.w_x[1]-sb.w_x[2]) \
273 * (sb.w_x[3] - sb.w_x[4])
274
275 # Allocate arrays on first run
276 if i == firststep:
277 axial_strain = numpy.zeros(lastfile+1, dtype=numpy.float64)
278 deviatoric_stress =\
279 numpy.zeros(lastfile+1, dtype=numpy.float64)
280 volumetric_strain =\
281 numpy.zeros(lastfile+1, dtype=numpy.float64)
282
283 w0pos0 = sb.w_x[0]
284 vol0 = vol
285
286 sigma1 = sb.w_force[0]/\
287 ((sb.w_x[1]-sb.w_x[2])*(sb.w_x[3]-sb.w_x[4]))
288
289 axial_strain[i] = (w0pos0 - sb.w_x[0])/w0pos0
290 volumetric_strain[i] = (vol0-vol)/vol0
291 deviatoric_stress[i] = sigma1 / sb.w_sigma0[1]
292
293
294 # Plotting
295 if outformat != 'txt':
296
297 # linear plot of deviatoric stress
298 ax1 = plt.subplot2grid((2, 1), (0, 0))
299 ax1.set_xlabel('Axial strain, $\gamma_1$, [-]')
300 ax1.set_ylabel('Deviatoric stress, $\sigma_1 - \sigma_3$, [Pa]')
301 ax1.plot(axial_strain, deviatoric_stress, '+-')
302 #ax1.legend()
303 ax1.grid()
304
305 #ax2 = plt.subplot2grid((2, 2), (1, 0))
306 #ax2.set_xlabel('Time [s]')
307 #ax2.set_ylabel('Force [N]')
308 #ax2.plot(t, wforce, '+-')
309
310 # semilog plot of log stress vs. void ratio
311 ax2 = plt.subplot2grid((2, 1), (1, 0))
312 ax2.set_xlabel('Axial strain, $\gamma_1$ [-]')
313 ax2.set_ylabel('Volumetric strain, $\gamma_v$, [-]')
314 ax2.plot(axial_strain, volumetric_strain, '+-')
315 ax2.grid()
316
317 if xlim:
318 ax1.set_xlim(xlim)
319 ax2.set_xlim(xlim)
320
321 elif method == 'shear':
322
323 # Read stress values from simulation binaries
324 for i in numpy.arange(firststep, lastfile+1):
325 sb.readstep(i, verbose=False)
326
327 # First iteration: Allocate arrays and find constant values
328 if i == firststep:
329 # Shear displacement
330 xdisp = numpy.zeros(lastfile+1, dtype=numpy.float64)
331
332 # Normal stress
333 sigma_eff = numpy.zeros(lastfile+1, dtype=numpy.float64)
334
335 # Normal stress
336 sigma_def = numpy.zeros(lastfile+1, dtype=numpy.float64)
337
338 # Shear stress
339 tau = numpy.zeros(lastfile+1, dtype=numpy.float64)
340
341 # Upper wall position
342 dilation = numpy.zeros(lastfile+1, dtype=numpy.float64)
343
344 # Peak shear stress
345 tau_p = 0.0
346
347 # Shear strain value of peak sh. stress
348 tau_p_shearstrain = 0.0
349
350 fixvel = numpy.nonzero(sb.fixvel > 0.0)
351 #fixvel_upper = numpy.nonzero(sb.vel[fixvel, 0] > 0.0)
352 shearvel = sb.vel[fixvel, 0].max()
353 w_x0 = sb.w_x[0] # Original height
354 A = sb.L[0] * sb.L[1] # Upper surface area
355
356 if i == firststep+1:
357 w_x0 = sb.w_x[0] # Original height
358
359 # Summation of shear stress contributions
360 for j in fixvel[0]:
361 if sb.vel[j, 0] > 0.0:
362 tau[i] += -sb.force[j, 0]/A
363
364 if i > 0:
365 xdisp[i] = xdisp[i-1] + sb.time_file_dt[0]*shearvel
366 sigma_eff[i] = sb.w_force[0]/A
367 sigma_def[i] = sb.w_sigma0[0]
368
369 # dilation in meters
370 #dilation[i] = sb.w_x[0] - w_x0
371
372 # dilation in percent
373 #dilation[i] = (sb.w_x[0] - w_x0)/w_x0 * 100.0 # dilation in percent
374
375 # dilation in number of mean particle diameters
376 d_bar = numpy.mean(self.radius)*2.0
377 if numpy.isnan(d_bar):
378 print('No radii in self.radius, attempting to read first '
379 + 'file')
380 self.readfirst()
381 d_bar = numpy.mean(self.radius)*2.0
382 dilation[i] = (sb.w_x[0] - w_x0)/d_bar
383
384 # Test if this was the max. shear stress
385 if tau[i] > tau_p:
386 tau_p = tau[i]
387 tau_p_shearstrain = xdisp[i]/w_x0
388
389 shear_strain = xdisp/w_x0
390
391 # Copy values so they can be modified during smoothing
392 shear_strain_smooth = shear_strain
393 tau_smooth = tau
394 sigma_def_smooth = sigma_def
395
396 # Optionally smooth the shear stress
397 if smoothing > 2:
398
399 if smoothing_window not in ['flat', 'hanning', 'hamming',
400 'bartlett', 'blackman']:
401 raise ValueError
402
403 s = numpy.r_[2*tau[0]-tau[smoothing:1:-1], tau,
404 2*tau[-1]-tau[-1:-smoothing:-1]]
405
406 if smoothing_window == 'flat': # moving average
407 w = numpy.ones(smoothing, 'd')
408 else:
409 w = getattr(self.np, smoothing_window)(smoothing)
410 y = numpy.convolve(w/w.sum(), s, mode='same')
411 tau_smooth = y[smoothing-1:-smoothing+1]
412
413 # Plot stresses
414 if outformat != 'txt':
415 shearinfo = "$\\tau_p$={:.3} Pa at $\gamma$={:.3}".format(\
416 tau_p, tau_p_shearstrain)
417 fig.text(0.01, 0.01, shearinfo, horizontalalignment='left',
418 fontproperties=FontProperties(size=14))
419 ax1 = plt.subplot2grid((2, 1), (0, 0))
420 ax1.set_xlabel('Shear strain [-]')
421 ax1.set_ylabel('Shear friction $\\tau/\\sigma_0$ [-]')
422 if smoothing > 2:
423 ax1.plot(shear_strain_smooth[1:-(smoothing+1)/2],
424 tau_smooth[1:-(smoothing+1)/2] /
425 sigma_def_smooth[1:-(smoothing+1)/2],
426 '-', label="$\\tau/\\sigma_0$")
427 else:
428 ax1.plot(shear_strain[1:],\
429 tau[1:]/sigma_def[1:],\
430 '-', label="$\\tau/\\sigma_0$")
431 ax1.grid()
432
433 # Plot dilation
434 ax2 = plt.subplot2grid((2, 1), (1, 0))
435 ax2.set_xlabel('Shear strain [-]')
436 ax2.set_ylabel('Dilation, $\Delta h/(2\\bar{r})$ [m]')
437 if smoothing > 2:
438 ax2.plot(shear_strain_smooth[1:-(smoothing+1)/2],
439 dilation[1:-(smoothing+1)/2], '-')
440 else:
441 ax2.plot(shear_strain, dilation, '-')
442 ax2.grid()
443
444 if xlim:
445 ax1.set_xlim(xlim)
446 ax2.set_xlim(xlim)
447
448 fig.tight_layout()
449
450 else:
451 # Write values to textfile
452 filename = "shear-stresses-{0}.txt".format(self.sid)
453 fh = None
454 try:
455 fh = open(filename, "w")
456 for i in numpy.arange(firststep, lastfile+1):
457 # format: shear distance [mm], sigma [kPa], tau [kPa],
458 # Dilation [%]
459 fh.write("{0}\t{1}\t{2}\t{3}\n"
460 .format(xdisp[i], sigma_eff[i]/1000.0,
461 tau[i]/1000.0, dilation[i]))
462 finally:
463 if fh is not None:
464 fh.close()
465
466 elif method == 'shear-displacement':
467
468 time = numpy.zeros(lastfile+1, dtype=numpy.float64)
469 # Read stress values from simulation binaries
470 for i in numpy.arange(firststep, lastfile+1):
471 sb.readstep(i, verbose=False)
472
473 # First iteration: Allocate arrays and find constant values
474 if i == firststep:
475
476 # Shear displacement
477 xdisp = numpy.zeros(lastfile+1, dtype=numpy.float64)
478
479 # Normal stress
480 sigma_eff = numpy.zeros(lastfile+1, dtype=numpy.float64)
481
482 # Normal stress
483 sigma_def = numpy.zeros(lastfile+1, dtype=numpy.float64)
484
485 # Shear stress
486 tau_eff = numpy.zeros(lastfile+1, dtype=numpy.float64)
487
488 # Upper wall position
489 dilation = numpy.zeros(lastfile+1, dtype=numpy.float64)
490
491 # Mean porosity
492 phi_bar = numpy.zeros(lastfile+1, dtype=numpy.float64)
493
494 # Mean fluid pressure
495 p_f_bar = numpy.zeros(lastfile+1, dtype=numpy.float64)
496 p_f_top = numpy.zeros(lastfile+1, dtype=numpy.float64)
497
498 # Upper wall position
499 tau_p = 0.0 # Peak shear stress
500 # Shear strain value of peak sh. stress
501 tau_p_shearstrain = 0.0
502
503 fixvel = numpy.nonzero(sb.fixvel > 0.0)
504 #fixvel_upper=numpy.nonzero(sb.vel[fixvel, 0] > 0.0)
505 w_x0 = sb.w_x[0] # Original height
506 A = sb.L[0]*sb.L[1] # Upper surface area
507
508 d_bar = numpy.mean(sb.radius)*2.0
509
510 # Shear velocity
511 v = numpy.zeros(lastfile+1, dtype=numpy.float64)
512
513 time[i] = sb.time_current[0]
514
515 if i == firststep+1:
516 w_x0 = sb.w_x[0] # Original height
517
518 # Summation of shear stress contributions
519 for j in fixvel[0]:
520 if sb.vel[j, 0] > 0.0:
521 tau_eff[i] += -sb.force[j, 0]/A
522
523 if i > 0:
524 xdisp[i] = sb.xyzsum[fixvel, 0].max()
525 v[i] = sb.vel[fixvel, 0].max()
526
527 sigma_eff[i] = sb.w_force[0]/A
528 sigma_def[i] = sb.currentNormalStress()
529
530 # dilation in number of mean particle diameters
531 dilation[i] = (sb.w_x[0] - w_x0)/d_bar
532
533 wall0_iz = int(sb.w_x[0]/(sb.L[2]/sb.num[2]))
534
535 if self.fluid:
536 if i > 0:
537 phi_bar[i] = numpy.mean(sb.phi[:, :, 0:wall0_iz])
538 if i == firststep+1:
539 phi_bar[0] = phi_bar[1]
540 p_f_bar[i] = numpy.mean(sb.p_f[:, :, 0:wall0_iz])
541 p_f_top[i] = sb.p_f[0, 0, -1]
542
543 # Test if this was the max. shear stress
544 if tau_eff[i] > tau_p:
545 tau_p = tau_eff[i]
546 tau_p_shearstrain = xdisp[i]/w_x0
547
548 shear_strain = xdisp/w_x0
549
550 # Plot stresses
551 if outformat != 'txt':
552 if figsize:
553 fig = plt.figure(figsize=figsize)
554 else:
555 fig = plt.figure(figsize=(8, 12))
556
557 # Upper plot
558 ax1 = plt.subplot(3, 1, 1)
559 ax1.plot(time, xdisp, 'k', label='Displacement')
560 ax1.set_ylabel('Horizontal displacement [m]')
561
562 ax2 = ax1.twinx()
563
564 #ax2color = '#666666'
565 ax2color = 'blue'
566 if self.fluid:
567 ax2.plot(time, phi_bar, color=ax2color, label='Porosity')
568 ax2.set_ylabel('Mean porosity $\\bar{\\phi}$ [-]')
569 else:
570 ax2.plot(time, dilation, color=ax2color, label='Dilation')
571 ax2.set_ylabel('Dilation, $\Delta h/(2\\bar{r})$ [-]')
572 for tl in ax2.get_yticklabels():
573 tl.set_color(ax2color)
574
575 # Middle plot
576 ax5 = plt.subplot(3, 1, 2, sharex=ax1)
577 ax5.semilogy(time[1:], v[1:], label='Shear velocity')
578 ax5.set_ylabel('Shear velocity [ms$^{-1}$]')
579
580 # shade stick periods
581 collection = \
582 matplotlib.collections.BrokenBarHCollection.span_where(
583 time, ymin=1.0e-7, ymax=1.0,
584 where=numpy.isclose(v, 0.0),
585 facecolor='black', alpha=0.2,
586 linewidth=0)
587 ax5.add_collection(collection)
588
589 # Lower plot
590 ax3 = plt.subplot(3, 1, 3, sharex=ax1)
591 if sb.w_sigma0_A > 1.0e-3:
592 lns0 = ax3.plot(time, sigma_def/1000.0,
593 '-k', label="$\\sigma_0$")
594 lns1 = ax3.plot(time, sigma_eff/1000.0,
595 '--k', label="$\\sigma'$")
596 lns2 = ax3.plot(time, numpy.ones_like(time)*sb.w_tau_x/1000.0,
597 '-r', label="$\\tau$")
598 lns3 = ax3.plot(time, tau_eff/1000.0,
599 '--r', label="$\\tau'$")
600 ax3.set_ylabel('Stress [kPa]')
601 else:
602 ax3.plot(time, tau_eff/sb.w_sigma0[0],
603 '-k', label="$Shear friction$")
604 ax3.plot([0, time[-1]],
605 [sb.w_tau_x/sigma_def, sb.w_tau_x/sigma_def],
606 '--k', label="$Applied shear friction$")
607 ax3.set_ylabel('Shear friction $\\tau\'/\\sigma_0$ [-]')
608 # axis limits
609 ax3.set_ylim([sb.w_tau_x/sigma_def[0]*0.5,
610 sb.w_tau_x/sigma_def[0]*1.5])
611
612 if self.fluid:
613 ax4 = ax3.twinx()
614 #ax4color = '#666666'
615 ax4color = ax2color
616 lns4 = ax4.plot(time, p_f_top/1000.0, '-', color=ax4color,
617 label='$p_\\text{f}^\\text{forcing}$')
618 lns5 = ax4.plot(time, p_f_bar/1000.0, '--', color=ax4color,
619 label='$\\bar{p}_\\text{f}$')
620 ax4.set_ylabel('Mean fluid pressure '
621 + '$\\bar{p_\\text{f}}$ [kPa]')
622 for tl in ax4.get_yticklabels():
623 tl.set_color(ax4color)
624 if sb.w_sigma0_A > 1.0e-3:
625 #ax4.legend(loc='upper right')
626 lns = lns0+lns1+lns2+lns3+lns4+lns5
627 labs = [l.get_label() for l in lns]
628 ax4.legend(lns, labs, loc='upper right',
629 fancybox=True, framealpha=legend_alpha)
630 if xlim:
631 ax4.set_xlim(xlim)
632
633 # aesthetics
634 ax3.set_xlabel('Time [s]')
635
636 ax1.grid()
637 ax3.grid()
638 ax5.grid()
639
640 if xlim:
641 ax1.set_xlim(xlim)
642 ax2.set_xlim(xlim)
643 ax3.set_xlim(xlim)
644 ax5.set_xlim(xlim)
645
646 plt.setp(ax1.get_xticklabels(), visible=False)
647 plt.setp(ax5.get_xticklabels(), visible=False)
648 fig.tight_layout()
649 plt.subplots_adjust(hspace=0.05)
650
651 elif method == 'rate-dependence':
652
653 if figsize:
654 fig = plt.figure(figsize=figsize)
655 else:
656 fig = plt.figure(figsize=(8, 6))
657
658 tau = numpy.empty(sb.status())
659 N = numpy.empty(sb.status())
660 #v = numpy.empty(sb.status())
661 shearstrainrate = numpy.empty(sb.status())
662 shearstrain = numpy.empty(sb.status())
663 for i in numpy.arange(firststep, sb.status()):
664 sb.readstep(i+1, verbose=False)
665 #tau = sb.shearStress()
666 tau[i] = sb.w_tau_x # defined shear stress
667 N[i] = sb.currentNormalStress() # defined normal stress
668 shearstrainrate[i] = sb.shearStrainRate()
669 shearstrain[i] = sb.shearStrain()
670
671 # remove nonzero sliding velocities and their associated values
672 idx = numpy.nonzero(shearstrainrate)
673 shearstrainrate_nonzero = shearstrainrate[idx]
674 tau_nonzero = tau[idx]
675 N_nonzero = N[idx]
676 shearstrain_nonzero = shearstrain[idx]
677
678 ax1 = plt.subplot(111)
679 #ax1.semilogy(N/1000., v)
680 #ax1.semilogy(tau_nonzero/N_nonzero, v_nonzero, '+k')
681 #ax1.plot(tau/N, v, '.')
682 friction = tau_nonzero/N_nonzero
683 #CS = ax1.scatter(friction, v_nonzero, c=shearstrain_nonzero,
684 #linewidth=0)
685 if cmap:
686 CS = ax1.scatter(friction, shearstrainrate_nonzero,
687 c=shearstrain_nonzero, linewidth=0.1,
688 cmap=cmap)
689 else:
690 CS = ax1.scatter(friction, shearstrainrate_nonzero,
691 c=shearstrain_nonzero, linewidth=0.1,
692 cmap=matplotlib.cm.get_cmap('afmhot'))
693 ax1.set_yscale('log')
694 x_min = numpy.floor(numpy.min(friction))
695 x_max = numpy.ceil(numpy.max(friction))
696 ax1.set_xlim([x_min, x_max])
697 y_min = numpy.min(shearstrainrate_nonzero)*0.5
698 y_max = numpy.max(shearstrainrate_nonzero)*2.0
699 ax1.set_ylim([y_min, y_max])
700
701 cb = plt.colorbar(CS)
702 cb.set_label('Shear strain $\\gamma$ [-]')
703
704 ax1.set_xlabel('Friction $\\tau/N$ [-]')
705 ax1.set_ylabel('Shear strain rate $\\dot{\\gamma}$ [s$^{-1}$]')
706
707 elif method == 'inertia':
708
709 t = numpy.zeros(sb.status())
710 I = numpy.zeros(sb.status())
711
712 for i in numpy.arange(firststep, sb.status()):
713 sb.readstep(i, verbose=False)
714 t[i] = sb.currentTime()
715 I[i] = sb.inertiaParameterPlanarShear()
716
717 # Plotting
718 if outformat != 'txt':
719
720 if xlim:
721 ax1.set_xlim(xlim)
722
723 # linear plot of deviatoric stress
724 ax1 = plt.subplot2grid((1, 1), (0, 0))
725 ax1.set_xlabel('Time $t$ [s]')
726 ax1.set_ylabel('Inertia parameter $I$ [-]')
727 ax1.semilogy(t, I)
728 #ax1.legend()
729 ax1.grid()
730
731 elif method == 'mean-fluid-pressure':
732
733 # Read pressure values from simulation binaries
734 for i in numpy.arange(firststep, lastfile+1):
735 sb.readstep(i, verbose=False)
736
737 # Allocate arrays on first run
738 if i == firststep:
739 p_mean = numpy.zeros(lastfile+1, dtype=numpy.float64)
740
741 p_mean[i] = numpy.mean(sb.p_f)
742
743 t = numpy.linspace(0.0, sb.time_current, lastfile+1)
744
745 # Plotting
746 if outformat != 'txt':
747
748 if xlim:
749 ax1.set_xlim(xlim)
750
751 # linear plot of deviatoric stress
752 ax1 = plt.subplot2grid((1, 1), (0, 0))
753 ax1.set_xlabel('Time $t$, [s]')
754 ax1.set_ylabel('Mean fluid pressure, $\\bar{p}_f$, [kPa]')
755 ax1.plot(t, p_mean/1000.0, '+-')
756 #ax1.legend()
757 ax1.grid()
758
759 elif method == 'fluid-pressure':
760
761 if figsize:
762 fig = plt.figure(figsize=figsize)
763 else:
764 fig = plt.figure(figsize=(8, 6))
765
766 sb.readfirst(verbose=False)
767
768 # cell midpoint cell positions
769 zpos_c = numpy.zeros(sb.num[2])
770 dz = sb.L[2]/sb.num[2]
771 for i in numpy.arange(sb.num[2]):
772 zpos_c[i] = i*dz + 0.5*dz
773
774 shear_strain = numpy.zeros(sb.status())
775 pres = numpy.zeros((sb.num[2], sb.status()))
776
777 # Read pressure values from simulation binaries
778 for i in numpy.arange(firststep, sb.status()):
779 sb.readstep(i, verbose=False)
780 pres[:, i] = numpy.average(numpy.average(sb.p_f, axis=0), axis=0)
781 shear_strain[i] = sb.shearStrain()
782 t = numpy.linspace(0.0, sb.time_current, lastfile+1)
783
784 # Plotting
785 if outformat != 'txt':
786
787 ax = plt.subplot(1, 1, 1)
788
789 pres /= 1000.0 # Pa to kPa
790
791 if xlim:
792 sb.readstep(10, verbose=False)
793 gamma_per_i = sb.shearStrain()/10.0
794 i_min = int(xlim[0]/gamma_per_i)
795 i_max = int(xlim[1]/gamma_per_i)
796 pres = pres[:, i_min:i_max]
797 else:
798 i_min = 0
799 i_max = sb.status()
800 # use largest difference in p from 0 as +/- limit on colormap
801 p_ext = numpy.max(numpy.abs(pres))
802
803 if sb.wmode[0] == 3:
804 x = t
805 else:
806 x = shear_strain
807 if xlim:
808 x = x[i_min:i_max]
809 if cmap:
810 im1 = ax.pcolormesh(x, zpos_c, pres, cmap=cmap,
811 vmin=-p_ext, vmax=p_ext,
812 rasterized=True)
813 else:
814 im1 = ax.pcolormesh(x, zpos_c, pres,
815 cmap=matplotlib.cm.get_cmap('RdBu_r'),
816 vmin=-p_ext, vmax=p_ext,
817 rasterized=True)
818 ax.set_xlim([0, numpy.max(x)])
819 if sb.w_x[0] < sb.L[2]:
820 ax.set_ylim([zpos_c[0], sb.w_x[0]])
821 else:
822 ax.set_ylim([zpos_c[0], zpos_c[-1]])
823 if sb.wmode[0] == 3:
824 ax.set_xlabel('Time $t$ [s]')
825 else:
826 ax.set_xlabel('Shear strain $\\gamma$ [-]')
827 ax.set_ylabel('Vertical position $z$ [m]')
828
829 if xlim:
830 ax.set_xlim([x[0], x[-1]])
831
832 # for article2
833 ax.set_ylim([zpos_c[0], zpos_c[9]])
834
835 cb = plt.colorbar(im1)
836 cb.set_label('$p_\\text{f}$ [kPa]')
837 cb.solids.set_rasterized(True)
838 plt.tight_layout()
839
840 elif method == 'porosity':
841
842 sb.readfirst(verbose=False)
843 if not sb.fluid:
844 raise Exception('Porosities can only be visualized in wet ' +
845 'simulations')
846
847 wall0_iz = int(sb.w_x[0]/(sb.L[2]/sb.num[2]))
848
849 # cell midpoint cell positions
850 zpos_c = numpy.zeros(sb.num[2])
851 dz = sb.L[2]/sb.num[2]
852 for i in numpy.arange(firststep, sb.num[2]):
853 zpos_c[i] = i*dz + 0.5*dz
854
855 shear_strain = numpy.zeros(sb.status())
856 poros = numpy.zeros((sb.num[2], sb.status()))
857
858 # Read pressure values from simulation binaries
859 for i in numpy.arange(firststep, sb.status()):
860 sb.readstep(i, verbose=False)
861 poros[:, i] = numpy.average(numpy.average(sb.phi, axis=0), axis=0)
862 shear_strain[i] = sb.shearStrain()
863 t = numpy.linspace(0.0, sb.time_current, lastfile+1)
864
865 # Plotting
866 if outformat != 'txt':
867
868 ax = plt.subplot(1, 1, 1)
869
870 poros_max = numpy.max(poros[0:wall0_iz-1, 1:])
871 poros_min = numpy.min(poros)
872
873 if sb.wmode[0] == 3:
874 x = t
875 else:
876 x = shear_strain
877 if cmap:
878 im1 = ax.pcolormesh(x, zpos_c, poros,
879 cmap=cmap,
880 vmin=poros_min, vmax=poros_max,
881 rasterized=True)
882 else:
883 im1 = ax.pcolormesh(x, zpos_c, poros,
884 cmap=matplotlib.cm.get_cmap('Blues_r'),
885 vmin=poros_min, vmax=poros_max,
886 rasterized=True)
887 ax.set_xlim([0, numpy.max(x)])
888 if sb.w_x[0] < sb.L[2]:
889 ax.set_ylim([zpos_c[0], sb.w_x[0]])
890 else:
891 ax.set_ylim([zpos_c[0], zpos_c[-1]])
892 if sb.wmode[0] == 3:
893 ax.set_xlabel('Time $t$ [s]')
894 else:
895 ax.set_xlabel('Shear strain $\\gamma$ [-]')
896 ax.set_ylabel('Vertical position $z$ [m]')
897
898 if xlim:
899 ax.set_xlim(xlim)
900
901 cb = plt.colorbar(im1)
902 cb.set_label('Mean horizontal porosity $\\bar{\phi}$ [-]')
903 cb.solids.set_rasterized(True)
904 plt.tight_layout()
905 plt.subplots_adjust(wspace=.05)
906
907 elif method == 'contacts':
908
909 for i in numpy.arange(sb.status()+1):
910 fn = "../output/{0}.output{1:0=5}.bin".format(self.sid, i)
911 sb.sid = self.sid + ".{:0=5}".format(i)
912 sb.readbin(fn, verbose=True)
913 if f_min and f_max:
914 sb.plotContacts(lower_limit=0.25, upper_limit=0.75,
915 outfolder='../img_out/',
916 f_min=f_min, f_max=f_max,
917 title="t={:.2f} s, $N$={:.0f} kPa"
918 .format(sb.currentTime(),
919 sb.currentNormalStress('defined')
920 /1000.))
921 else:
922 sb.plotContacts(lower_limit=0.25, upper_limit=0.75,
923 title="t={:.2f} s, $N$={:.0f} kPa"
924 .format(sb.currentTime(),
925 sb.currentNormalStress('defined')
926 /1000.), outfolder='../img_out/')
927
928 # render images to movie
929 subprocess.call('cd ../img_out/ && ' +
930 'ffmpeg -sameq -i {}.%05d-contacts.png '
931 .format(self.sid) +
932 '{}-contacts.mp4'.format(self.sid),
933 shell=True)
934
935 else:
936 print("Visualization type '" + method + "' not understood")
937 return
938
939 # Optional save of figure content
940 filename = ''
941 if xlim:
942 filename = '{0}-{1}-{3}.{2}'.format(self.sid, method, outformat,
943 xlim[-1])
944 else:
945 filename = '{0}-{1}.{2}'.format(self.sid, method, outformat)
946 if pickle:
947 pl.dump(fig, open(filename + '.pickle', 'wb'))
948
949 # Optional save of figure
950 if outformat != 'txt':
951 if savefig:
952 fig.savefig(filename)
953 print(filename)
954 fig.clf()
955 plt.close()
956 else:
957 plt.show()