plotting.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
---
plotting.py (43500B)
---
1 import math
2 import subprocess
3 import numpy
4 from .common import matplotlib, plt, py_mpl
5
6
7 class SimPlotting:
8 'matplotlib and gnuplot figures of simulation output.'
9
10 def forcechains(self, lc=200.0, uc=650.0, outformat='png', disp='2d'):
11 '''
12 Visualizes the force chains in the system from the magnitude of the
13 normal contact forces, and produces an image of them. Warning: Will
14 segfault if no contacts are found.
15
16 :param lc: Lower cutoff of contact forces. Contacts below are not
17 visualized
18 :type lc: float
19 :param uc: Upper cutoff of contact forces. Contacts above are
20 visualized with this value
21 :type uc: float
22 :param outformat: Format of output image. Possible values are
23 'interactive', 'png', 'epslatex', 'epslatex-color'
24 :type outformat: str
25 :param disp: Display forcechains in '2d' or '3d'
26 :type disp: str
27 '''
28
29 self.writebin(verbose=False)
30
31 nd = ''
32 if disp == '2d':
33 nd = '-2d '
34
35 subprocess.call("cd .. && ./forcechains " + nd + "-f " + outformat
36 + " -lc " + str(lc) + " -uc " + str(uc)
37 + " input/" + self.sid + ".bin > python/tmp.gp",
38 shell=True)
39 subprocess.call("gnuplot tmp.gp && rm tmp.gp", shell=True)
40
41 def forcechainsRose(self, lower_limit=0.25, graphics_format='pdf'):
42 '''
43 Visualize trend and plunge angles of the strongest force chains in a
44 rose plot. The plots are saved in the current folder with the name
45 'fc-<simulation id>-rose.pdf'.
46
47 :param lower_limit: Do not visualize force chains below this relative
48 contact force magnitude, in ]0;1[
49 :type lower_limit: float
50 :param graphics_format: Save the plot in this format
51 :type graphics_format: str
52 '''
53 self.writebin(verbose=False)
54
55 subprocess.call("cd .. && ./forcechains -f txt input/" + self.sid \
56 + ".bin > python/fc-tmp.txt", shell=True)
57
58 # data will have the shape (numcontacts, 7)
59 data = numpy.loadtxt("fc-tmp.txt", skiprows=1)
60
61 # find the max. value of the normal force
62 f_n_max = numpy.amax(data[:, 6])
63
64 # specify the lower limit of force chains to do statistics on
65 f_n_lim = lower_limit * f_n_max * 0.6
66
67 # find the indexes of these contacts
68 I = numpy.nonzero(data[:, 6] > f_n_lim)
69
70 # loop through these contacts and find the strike and dip of the
71 # contacts
72 strikelist = [] # strike direction of the normal vector, [0:360[
73 diplist = [] # dip of the normal vector, [0:90]
74 for i in I[0]:
75
76 x1 = data[i, 0]
77 y1 = data[i, 1]
78 z1 = data[i, 2]
79 x2 = data[i, 3]
80 y2 = data[i, 4]
81 z2 = data[i, 5]
82
83 if z1 < z2:
84 xlower = x1; ylower = y1; zlower = z1
85 xupper = x2; yupper = y2; zupper = z2
86 else:
87 xlower = x2; ylower = y2; zlower = z2
88 xupper = x1; yupper = y1; zupper = z1
89
90 # Vector pointing downwards
91 dx = xlower - xupper
92 dy = ylower - yupper
93 dhoriz = numpy.sqrt(dx**2 + dy**2)
94
95 # Find dip angle
96 diplist.append(math.degrees(math.atan((zupper - zlower)/dhoriz)))
97
98 # Find strike angle
99 if ylower >= yupper: # in first two quadrants
100 strikelist.append(math.acos(dx/dhoriz))
101 else:
102 strikelist.append(2.0*numpy.pi - math.acos(dx/dhoriz))
103
104
105 plt.figure(figsize=[4, 4])
106 ax = plt.subplot(111, polar=True)
107 ax.scatter(strikelist, diplist, c='k', marker='+')
108 ax.set_rmax(90)
109 ax.set_rticks([])
110 plt.savefig('fc-' + self.sid + '-rose.' + graphics_format,\
111 transparent=True)
112
113 subprocess.call('rm fc-tmp.txt', shell=True)
114
115 def bondsRose(self, graphics_format='pdf'):
116 '''
117 Visualize the trend and plunge angles of the bond pairs in a rose plot.
118 The plot is saved in the current folder as
119 'bonds-<simulation id>-rose.<graphics_format>'.
120
121 :param graphics_format: Save the plot in this format
122 :type graphics_format: str
123 '''
124 if not py_mpl:
125 print('Error: matplotlib module not found, cannot bondsRose.')
126 return
127 # loop through these contacts and find the strike and dip of the
128 # contacts
129 strikelist = [] # strike direction of the normal vector, [0:360[
130 diplist = [] # dip of the normal vector, [0:90]
131 for n in numpy.arange(self.nb0):
132
133 i = self.bonds[n, 0]
134 j = self.bonds[n, 1]
135
136 x1 = self.x[i, 0]
137 y1 = self.x[i, 1]
138 z1 = self.x[i, 2]
139 x2 = self.x[j, 0]
140 y2 = self.x[j, 1]
141 z2 = self.x[j, 2]
142
143 if z1 < z2:
144 xlower = x1; ylower = y1; zlower = z1
145 xupper = x2; yupper = y2; zupper = z2
146 else:
147 xlower = x2; ylower = y2; zlower = z2
148 xupper = x1; yupper = y1; zupper = z1
149
150 # Vector pointing downwards
151 dx = xlower - xupper
152 dy = ylower - yupper
153 dhoriz = numpy.sqrt(dx**2 + dy**2)
154
155 # Find dip angle
156 diplist.append(math.degrees(math.atan((zupper - zlower)/dhoriz)))
157
158 # Find strike angle
159 if ylower >= yupper: # in first two quadrants
160 strikelist.append(math.acos(dx/dhoriz))
161 else:
162 strikelist.append(2.0*numpy.pi - math.acos(dx/dhoriz))
163
164 plt.figure(figsize=[4, 4])
165 ax = plt.subplot(111, polar=True)
166 ax.scatter(strikelist, diplist, c='k', marker='+')
167 ax.set_rmax(90)
168 ax.set_rticks([])
169 plt.savefig('bonds-' + self.sid + '-rose.' + graphics_format,\
170 transparent=True)
171
172 def sheardisp(self, graphics_format='pdf', zslices=32):
173 '''
174 Plot the particle x-axis displacement against the original vertical
175 particle position. The plot is saved in the current directory with the
176 file name '<simulation id>-sheardisp.<graphics_format>'.
177
178 :param graphics_format: Save the plot in this format
179 :type graphics_format: str
180 '''
181 if not py_mpl:
182 print('Error: matplotlib module not found, cannot sheardisp.')
183 return
184
185 # Bin data and error bars for alternative visualization
186 h_total = numpy.max(self.x[:, 2]) - numpy.min(self.x[:, 2])
187 h_slice = h_total / zslices
188
189 zpos = numpy.zeros(zslices)
190 xdisp = numpy.zeros(zslices)
191 err = numpy.zeros(zslices)
192
193 for iz in range(zslices):
194
195 # Find upper and lower boundaries of bin
196 zlower = iz * h_slice
197 zupper = zlower + h_slice
198
199 # Save depth
200 zpos[iz] = zlower + 0.5*h_slice
201
202 # Find particle indexes within that slice
203 I = numpy.nonzero((self.x[:, 2] > zlower) & (self.x[:, 2] < zupper))
204
205 # Save mean x displacement
206 xdisp[iz] = numpy.mean(self.xyzsum[I, 0])
207
208 # Save x displacement standard deviation
209 err[iz] = numpy.std(self.xyzsum[I, 0])
210
211 plt.figure(figsize=[4, 4])
212 ax = plt.subplot(111)
213 ax.scatter(self.xyzsum[:, 0], self.x[:, 2], c='gray', marker='+')
214 ax.errorbar(xdisp, zpos, xerr=err,
215 c='black', linestyle='-', linewidth=1.4)
216 ax.set_xlabel("Horizontal particle displacement, [m]")
217 ax.set_ylabel("Vertical position, [m]")
218 plt.savefig(self.sid + '-sheardisp.' + graphics_format,
219 transparent=True)
220
221 def porosities(self, graphics_format='pdf', zslices=16):
222 '''
223 Plot the averaged porosities with depth. The plot is saved in the format
224 '<simulation id>-porosity.<graphics_format>'.
225
226 :param graphics_format: Save the plot in this format
227 :type graphics_format: str
228 :param zslices: The number of points along the vertical axis to sample
229 the porosity in
230 :type zslices: int
231 '''
232 if not py_mpl:
233 print('Error: matplotlib module not found, cannot sheardisp.')
234 return
235
236 porosity, depth = self.porosity(zslices)
237
238 plt.figure(figsize=[4, 4])
239 ax = plt.subplot(111)
240 ax.plot(porosity, depth, c='black', linestyle='-', linewidth=1.4)
241 ax.set_xlabel('Horizontally averaged porosity, [-]')
242 ax.set_ylabel('Vertical position, [m]')
243 plt.savefig(self.sid + '-porositiy.' + graphics_format,
244 transparent=True)
245
246 def thinsection_x1x3(self, x2='center', graphics_format='png', cbmax=None,
247 arrowscale=0.01, velarrowscale=1.0, slipscale=1.0,
248 verbose=False):
249 '''
250 Produce a 2D image of particles on a x1,x3 plane, intersecting the
251 second axis at x2. Output is saved as '<sid>-ts-x1x3.txt' in the
252 current folder.
253
254 An upper limit to the pressure color bar range can be set by the
255 cbmax parameter.
256
257 The data can be plotted in gnuplot with:
258 gnuplot> set size ratio -1
259 gnuplot> set palette defined (0 "blue", 0.5 "gray", 1 "red")
260 gnuplot> plot '<sid>-ts-x1x3.txt' with circles palette fs \
261 transparent solid 0.4 noborder
262
263 This function also saves a plot of the inter-particle slip angles.
264
265 :param x2: The position along the second axis of the intersecting plane
266 :type x2: foat
267 :param graphics_format: Save the slip angle plot in this format
268 :type graphics_format: str
269 :param cbmax: The maximal value of the pressure color bar range
270 :type cbmax: float
271 :param arrowscale: Scale the rotational arrows by this value
272 :type arrowscale: float
273 :param velarrowscale: Scale the translational arrows by this value
274 :type velarrowscale: float
275 :param slipscale: Scale the slip arrows by this value
276 :type slipscale: float
277 :param verbose: Show function output during calculations
278 :type verbose: bool
279 '''
280
281 if not py_mpl:
282 print('Error: matplotlib module not found (thinsection_x1x3).')
283 return
284
285 if x2 == 'center':
286 x2 = (self.L[1] - self.origo[1]) / 2.0
287
288 # Initialize plot circle positionsr, radii and pressures
289 ilist = []
290 xlist = []
291 ylist = []
292 rlist = []
293 plist = []
294 pmax = 0.0
295 rmax = 0.0
296 axlist = []
297 aylist = []
298 daxlist = []
299 daylist = []
300 dvxlist = []
301 dvylist = []
302 # Black circle at periphery of particles with angvel[:, 1] > 0.0
303 cxlist = []
304 cylist = []
305 crlist = []
306
307 # Loop over all particles, find intersections
308 for i in range(self.np):
309
310 delta = abs(self.x[i, 1] - x2) # distance between centre and plane
311
312 if delta < self.radius[i]: # if the sphere intersects the plane
313
314 # Store particle index
315 ilist.append(i)
316
317 # Store position on plane
318 xlist.append(self.x[i, 0])
319 ylist.append(self.x[i, 2])
320
321 # Store radius of intersection
322 r_circ = math.sqrt(self.radius[i]**2 - delta**2)
323 if r_circ > rmax:
324 rmax = r_circ
325 rlist.append(r_circ)
326
327 # Store pos. and radius if it is spinning around pos. y
328 if self.angvel[i, 1] > 0.0:
329 cxlist.append(self.x[i, 0])
330 cylist.append(self.x[i, 2])
331 crlist.append(r_circ)
332
333 # Store pressure
334 pval = self.p[i]
335 if cbmax != None:
336 if pval > cbmax:
337 pval = cbmax
338 plist.append(pval)
339
340 # Store rotational velocity data for arrows
341 # Save two arrows per particle
342 axlist.append(self.x[i, 0]) # x starting point of arrow
343 axlist.append(self.x[i, 0]) # x starting point of arrow
344
345 # y starting point of arrow
346 aylist.append(self.x[i, 2] + r_circ*0.5)
347
348 # y starting point of arrow
349 aylist.append(self.x[i, 2] - r_circ*0.5)
350
351 # delta x for arrow end point
352 daxlist.append(self.angvel[i, 1]*arrowscale)
353
354 # delta x for arrow end point
355 daxlist.append(-self.angvel[i, 1]*arrowscale)
356 daylist.append(0.0) # delta y for arrow end point
357 daylist.append(0.0) # delta y for arrow end point
358
359 # Store linear velocity data
360
361 # delta x for arrow end point
362 dvxlist.append(self.vel[i, 0]*velarrowscale)
363
364 # delta y for arrow end point
365 dvylist.append(self.vel[i, 2]*velarrowscale)
366
367 if r_circ > self.radius[i]:
368 raise Exception("Error, circle radius is larger than the "
369 "particle radius")
370 if self.p[i] > pmax:
371 pmax = self.p[i]
372
373 if verbose:
374 print("Max. pressure of intersecting spheres: " + str(pmax) + " Pa")
375 if cbmax != None:
376 print("Value limited to: " + str(cbmax) + " Pa")
377
378 # Save circle data
379 filename = '../gnuplot/data/' + self.sid + '-ts-x1x3.txt'
380 fh = None
381 try:
382 fh = open(filename, 'w')
383
384 for (x, y, r, p) in zip(xlist, ylist, rlist, plist):
385 fh.write("{}\t{}\t{}\t{}\n".format(x, y, r, p))
386
387 finally:
388 if fh is not None:
389 fh.close()
390
391 # Save circle data for articles spinning with pos. y
392 filename = '../gnuplot/data/' + self.sid + '-ts-x1x3-circ.txt'
393 fh = None
394 try:
395 fh = open(filename, 'w')
396
397 for (x, y, r) in zip(cxlist, cylist, crlist):
398 fh.write("{}\t{}\t{}\n".format(x, y, r))
399
400 finally:
401 if fh is not None:
402 fh.close()
403
404 # Save angular velocity data. The arrow lengths are normalized to max.
405 # radius
406 # Output format: x, y, deltax, deltay
407 # gnuplot> plot '-' using 1:2:3:4 with vectors head filled lt 2
408 filename = '../gnuplot/data/' + self.sid + '-ts-x1x3-arrows.txt'
409 fh = None
410 try:
411 fh = open(filename, 'w')
412
413 for (ax, ay, dax, day) in zip(axlist, aylist, daxlist, daylist):
414 fh.write("{}\t{}\t{}\t{}\n".format(ax, ay, dax, day))
415
416 finally:
417 if fh is not None:
418 fh.close()
419
420 # Save linear velocity data
421 # Output format: x, y, deltax, deltay
422 # gnuplot> plot '-' using 1:2:3:4 with vectors head filled lt 2
423 filename = '../gnuplot/data/' + self.sid + '-ts-x1x3-velarrows.txt'
424 fh = None
425 try:
426 fh = open(filename, 'w')
427
428 for (x, y, dvx, dvy) in zip(xlist, ylist, dvxlist, dvylist):
429 fh.write("{}\t{}\t{}\t{}\n".format(x, y, dvx, dvy))
430
431 finally:
432 if fh is not None:
433 fh.close()
434
435 # Check whether there are slips between the particles intersecting the
436 # plane
437 sxlist = []
438 sylist = []
439 dsxlist = []
440 dsylist = []
441 anglelist = [] # angle of the slip vector
442 slipvellist = [] # velocity of the slip
443 for i in ilist:
444
445 # Loop through other particles, and check whether they are in
446 # contact
447 for j in ilist:
448 #if i < j:
449 if i != j:
450
451 # positions
452 x_i = self.x[i, :]
453 x_j = self.x[j, :]
454
455 # radii
456 r_i = self.radius[i]
457 r_j = self.radius[j]
458
459 # Inter-particle vector
460 x_ij = x_i - x_j
461 x_ij_length = numpy.sqrt(x_ij.dot(x_ij))
462
463 # Check for overlap
464 if x_ij_length - (r_i + r_j) < 0.0:
465
466 # contact plane normal vector
467 n_ij = x_ij / x_ij_length
468
469 vel_i = self.vel[i, :]
470 vel_j = self.vel[j, :]
471 angvel_i = self.angvel[i, :]
472 angvel_j = self.angvel[j, :]
473
474 # Determine the tangential contact surface velocity in
475 # the x,z plane
476 dot_delta = (vel_i - vel_j) \
477 + r_i * numpy.cross(n_ij, angvel_i) \
478 + r_j * numpy.cross(n_ij, angvel_j)
479
480 # Subtract normal component to get tangential velocity
481 dot_delta_n = n_ij * numpy.dot(dot_delta, n_ij)
482 dot_delta_t = dot_delta - dot_delta_n
483
484 # Save slip velocity data for gnuplot
485 if dot_delta_t[0] != 0.0 or dot_delta_t[2] != 0.0:
486
487 # Center position of the contact
488 cpos = x_i - x_ij * 0.5
489
490 sxlist.append(cpos[0])
491 sylist.append(cpos[2])
492 dsxlist.append(dot_delta_t[0] * slipscale)
493 dsylist.append(dot_delta_t[2] * slipscale)
494 #anglelist.append(math.degrees(\
495 #math.atan(dot_delta_t[2]/dot_delta_t[0])))
496 anglelist.append(\
497 math.atan(dot_delta_t[2]/dot_delta_t[0]))
498 slipvellist.append(\
499 numpy.sqrt(dot_delta_t.dot(dot_delta_t)))
500
501
502 # Write slip lines to text file
503 filename = '../gnuplot/data/' + self.sid + '-ts-x1x3-slips.txt'
504 fh = None
505 try:
506 fh = open(filename, 'w')
507
508 for (sx, sy, dsx, dsy) in zip(sxlist, sylist, dsxlist, dsylist):
509 fh.write("{}\t{}\t{}\t{}\n".format(sx, sy, dsx, dsy))
510
511 finally:
512 if fh is not None:
513 fh.close()
514
515 # Plot thinsection with gnuplot script
516 gamma = self.shearStrain()
517 subprocess.call('''cd ../gnuplot/scripts && gnuplot -e "sid='{}'; ''' \
518 + '''gamma='{:.4}'; xmin='{}'; xmax='{}'; ymin='{}'; ''' \
519 + '''ymax='{}'" plotts.gp'''.format(\
520 self.sid, self.shearStrain(), self.origo[0], self.L[0], \
521 self.origo[2], self.L[2]), shell=True)
522
523 # Find all particles who have a slip velocity higher than slipvel
524 slipvellimit = 0.01
525 slipvels = numpy.nonzero(numpy.array(slipvellist) > slipvellimit)
526
527 # Bin slip angle data for histogram
528 binno = 36/2
529 hist_ang, bins_ang = numpy.histogram(numpy.array(anglelist)[slipvels],\
530 bins=binno, density=False)
531 center_ang = (bins_ang[:-1] + bins_ang[1:]) / 2.0
532
533 center_ang_mirr = numpy.concatenate((center_ang, center_ang + math.pi))
534 hist_ang_mirr = numpy.tile(hist_ang, 2)
535
536 # Write slip angles to text file
537 #numpy.savetxt(self.sid + '-ts-x1x3-slipangles.txt', zip(center_ang,\
538 #hist_ang), fmt="%f\t%f")
539
540 fig = plt.figure()
541 ax = fig.add_subplot(111, polar=True)
542 ax.bar(center_ang_mirr, hist_ang_mirr, width=30.0/180.0)
543 fig.savefig('../img_out/' + self.sid + '-ts-x1x3-slipangles.' +
544 graphics_format)
545 fig.clf()
546
547 def plotContacts(self, graphics_format='png', figsize=[4, 4], title=None,
548 lower_limit=0.0, upper_limit=1.0, alpha=1.0,
549 return_data=False, outfolder='.',
550 f_min=None, f_max=None, histogram=True,
551 forcechains=True):
552 '''
553 Plot current contact orientations on polar plot
554
555 :param lower_limit: Do not visualize force chains below this relative
556 contact force magnitude, in ]0;1[
557 :type lower_limit: float
558 :param upper_limit: Visualize force chains above this relative
559 contact force magnitude but cap color bar range, in ]0;1[
560 :type upper_limit: float
561 :param graphics_format: Save the plot in this format
562 :type graphics_format: str
563 '''
564
565 if not py_mpl:
566 print('Error: matplotlib module not found (plotContacts).')
567 return
568
569 self.writebin(verbose=False)
570
571 subprocess.call("cd .. && ./forcechains -f txt input/" + self.sid \
572 + ".bin > python/contacts-tmp.txt", shell=True)
573
574 # data will have the shape (numcontacts, 7)
575 data = numpy.loadtxt('contacts-tmp.txt', skiprows=1)
576
577 # find the max. value of the normal force
578 f_n_max = numpy.amax(data[:, 6])
579
580 # specify the lower limit of force chains to do statistics on
581 f_n_lim = lower_limit * f_n_max
582
583 if f_min:
584 f_n_lim = f_min
585 if f_max:
586 f_n_max = f_max
587
588 # find the indexes of these contacts
589 I = numpy.nonzero(data[:, 6] >= f_n_lim)
590
591 # loop through these contacts and find the strike and dip of the
592 # contacts
593
594 # strike direction of the normal vector, [0:360[
595 strikelist = numpy.empty(len(I[0]))
596 diplist = numpy.empty(len(I[0])) # dip of the normal vector, [0:90]
597 forcemagnitude = data[I, 6]
598 j = 0
599 for i in I[0]:
600
601 x1 = data[i, 0]
602 y1 = data[i, 1]
603 z1 = data[i, 2]
604 x2 = data[i, 3]
605 y2 = data[i, 4]
606 z2 = data[i, 5]
607
608 if z1 < z2:
609 xlower = x1; ylower = y1; zlower = z1
610 xupper = x2; yupper = y2; zupper = z2
611 else:
612 xlower = x2; ylower = y2; zlower = z2
613 xupper = x1; yupper = y1; zupper = z1
614
615 # Vector pointing downwards
616 dx = xlower - xupper
617 dy = ylower - yupper
618 dhoriz = numpy.sqrt(dx**2 + dy**2)
619
620 # Find dip angle
621 diplist[j] = numpy.degrees(numpy.arctan((zupper - zlower)/dhoriz))
622
623 # Find strike angle
624 if ylower >= yupper: # in first two quadrants
625 strikelist[j] = numpy.arccos(dx/dhoriz)
626 else:
627 strikelist[j] = 2.0*numpy.pi - numpy.arccos(dx/dhoriz)
628
629 j += 1
630
631 fig = plt.figure(figsize=figsize)
632 ax = plt.subplot(111, polar=True)
633 cs = ax.scatter(strikelist, 90. - diplist, marker='o',
634 c=forcemagnitude,
635 s=forcemagnitude/f_n_max*40.,
636 alpha=alpha,
637 edgecolors='none',
638 vmin=f_n_max*lower_limit,
639 vmax=f_n_max*upper_limit,
640 cmap=matplotlib.cm.get_cmap('afmhot_r'))
641 plt.colorbar(cs, extend='max')
642
643 # plot defined max compressive stress from tau/N ratio
644 ax.scatter(0., # prescribed stress
645 numpy.degrees(numpy.arctan(self.shearStress('defined')/
646 self.currentNormalStress('defined'))),
647 marker='o', c='none', edgecolor='blue', s=300)
648 ax.scatter(0., # actual stress
649 numpy.degrees(numpy.arctan(self.shearStress('effective')/
650 self.currentNormalStress('effective'))),
651 marker='+', color='blue', s=300)
652
653 ax.set_rmax(90)
654 ax.set_rticks([])
655
656 if title:
657 plt.title(title)
658 else:
659 plt.title('t={:.2f} s'.format(self.currentTime()))
660
661 #plt.tight_layout()
662 plt.savefig(outfolder + '/contacts-' + self.sid + '-' + \
663 str(self.time_step_count[0]) + '.' + \
664 graphics_format,\
665 transparent=False)
666
667 subprocess.call('rm contacts-tmp.txt', shell=True)
668
669 fig.clf()
670 if histogram:
671 #hist, bins = numpy.histogram(datadata[:, 6], bins=10)
672 _, _, _ = plt.hist(data[:, 6], alpha=0.75, facecolor='gray')
673 #plt.xlabel('$\\boldsymbol{f}_\text{n}$ [N]')
674 plt.yscale('log', nonposy='clip')
675 plt.xlabel('Contact load [N]')
676 plt.ylabel('Count $N$')
677 plt.grid(True)
678 plt.savefig(outfolder + '/contacts-hist-' + self.sid + '-' + \
679 str(self.time_step_count[0]) + '.' + \
680 graphics_format,\
681 transparent=False)
682 plt.clf()
683
684 # angle: 0 when vertical, 90 when horizontal
685 #hist, bins = numpy.histogram(datadata[:, 6], bins=10)
686 _, _, _ = plt.hist(90. - diplist, bins=range(0, 100, 10),
687 alpha=0.75, facecolor='gray')
688 theta_sigma1 = numpy.degrees(numpy.arctan(
689 self.currentNormalStress('defined')/\
690 self.shearStress('defined')))
691 plt.axvline(90. - theta_sigma1, color='k', linestyle='dashed',
692 linewidth=1)
693 plt.xlim([0, 90.])
694 plt.ylim([0, self.np/10])
695 #plt.xlabel('$\\boldsymbol{f}_\text{n}$ [N]')
696 plt.xlabel('Contact angle [deg]')
697 plt.ylabel('Count $N$')
698 plt.grid(True)
699 plt.savefig(outfolder + '/dip-' + self.sid + '-' + \
700 str(self.time_step_count[0]) + '.' + \
701 graphics_format,\
702 transparent=False)
703 plt.clf()
704
705 if forcechains:
706
707 #color = matplotlib.cm.spectral(data[:, 6]/f_n_max)
708 for i in I[0]:
709
710 x1 = data[i, 0]
711 #y1 = data[i, 1]
712 z1 = data[i, 2]
713 x2 = data[i, 3]
714 #y2 = data[i, 4]
715 z2 = data[i, 5]
716 f_n = data[i, 6]
717
718 lw_max = 1.0
719 if f_n >= f_n_max:
720 lw = lw_max
721 else:
722 lw = (f_n - f_n_lim)/(f_n_max - f_n_lim)*lw_max
723
724 plt.plot([x1, x2], [z1, z2], '-k', linewidth=lw)
725
726 axfc1 = plt.gca()
727 axfc1.spines['right'].set_visible(False)
728 axfc1.spines['left'].set_visible(False)
729 # Only show ticks on the left and bottom spines
730 axfc1.xaxis.set_ticks_position('none')
731 axfc1.yaxis.set_ticks_position('none')
732 #axfc1.set_xticklabels([])
733 #axfc1.set_yticklabels([])
734 axfc1.set_xlim([self.origo[0], self.L[0]])
735 axfc1.set_ylim([self.origo[2], self.L[2]])
736 axfc1.set_aspect('equal')
737
738 plt.xlabel('$x$ [m]')
739 plt.ylabel('$z$ [m]')
740 plt.grid(False)
741 plt.savefig(outfolder + '/fc-' + self.sid + '-' + \
742 str(self.time_step_count[0]) + '.' + \
743 graphics_format,\
744 transparent=False)
745
746 plt.close()
747
748 if return_data:
749 return data, strikelist, diplist, forcemagnitude, alpha, f_n_max
750
751 def plotFluidPressuresY(self, y=-1, graphics_format='png', verbose=True):
752 '''
753 Plot fluid pressures in a plane normal to the second axis.
754 The plot is saved in the current folder with the format
755 'p_f-<simulation id>-y<y value>.<graphics_format>'.
756
757 :param y: Plot pressures in fluid cells with these y axis values. If
758 this value is -1, the center y position is used.
759 :type y: int
760 :param graphics_format: Save the plot in this format
761 :type graphics_format: str
762 :param verbose: Print output filename after saving
763 :type verbose: bool
764
765 See also: :func:`writeFluidVTK()` and :func:`plotFluidPressuresZ()`
766 '''
767
768 if not py_mpl:
769 print('Error: matplotlib module not found (plotFluidPressuresY).')
770 return
771
772 if y == -1:
773 y = self.num[1]/2
774
775 plt.figure(figsize=[8, 8])
776 plt.title('Fluid pressures')
777 imgplt = plt.imshow(self.p_f[:, y, :].T, origin='lower')
778 imgplt.set_interpolation('nearest')
779 #imgplt.set_interpolation('bicubic')
780 #imgplt.set_cmap('hot')
781 plt.xlabel('$x_1$')
782 plt.ylabel('$x_3$')
783 plt.colorbar()
784 filename = 'p_f-' + self.sid + '-y' + str(y) + '.' + graphics_format
785 plt.savefig(filename, transparent=False)
786 if verbose:
787 print('saved to ' + filename)
788 plt.clf()
789 plt.close()
790
791 def plotFluidPressuresZ(self, z=-1, graphics_format='png', verbose=True):
792 '''
793 Plot fluid pressures in a plane normal to the third axis.
794 The plot is saved in the current folder with the format
795 'p_f-<simulation id>-z<z value>.<graphics_format>'.
796
797 :param z: Plot pressures in fluid cells with these z axis values. If
798 this value is -1, the center z position is used.
799 :type z: int
800 :param graphics_format: Save the plot in this format
801 :type graphics_format: str
802 :param verbose: Print output filename after saving
803 :type verbose: bool
804
805 See also: :func:`writeFluidVTK()` and :func:`plotFluidPressuresY()`
806 '''
807
808 if not py_mpl:
809 print('Error: matplotlib module not found (plotFluidPressuresZ).')
810 return
811
812 if z == -1:
813 z = self.num[2]/2
814
815 plt.figure(figsize=[8, 8])
816 plt.title('Fluid pressures')
817 imgplt = plt.imshow(self.p_f[:, :, z].T, origin='lower')
818 imgplt.set_interpolation('nearest')
819 #imgplt.set_interpolation('bicubic')
820 #imgplt.set_cmap('hot')
821 plt.xlabel('$x_1$')
822 plt.ylabel('$x_2$')
823 plt.colorbar()
824 filename = 'p_f-' + self.sid + '-z' + str(z) + '.' + graphics_format
825 plt.savefig(filename, transparent=False)
826 if verbose:
827 print('saved to ' + filename)
828 plt.clf()
829 plt.close()
830
831 def plotFluidVelocitiesY(self, y=-1, graphics_format='png', verbose=True):
832 '''
833 Plot fluid velocities in a plane normal to the second axis.
834 The plot is saved in the current folder with the format
835 'v_f-<simulation id>-z<z value>.<graphics_format>'.
836
837 :param y: Plot velocities in fluid cells with these y axis values. If
838 this value is -1, the center y position is used.
839 :type y: int
840 :param graphics_format: Save the plot in this format
841 :type graphics_format: str
842 :param verbose: Print output filename after saving
843 :type verbose: bool
844
845 See also: :func:`writeFluidVTK()` and :func:`plotFluidVelocitiesZ()`
846 '''
847
848 if not py_mpl:
849 print('Error: matplotlib module not found (plotFluidVelocitiesY).')
850 return
851
852 if y == -1:
853 y = self.num[1]/2
854
855 plt.title('Fluid velocities')
856 plt.figure(figsize=[8, 8])
857
858 plt.subplot(131)
859 imgplt = plt.imshow(self.v_f[:, y, :, 0].T, origin='lower')
860 imgplt.set_interpolation('nearest')
861 #imgplt.set_interpolation('bicubic')
862 #imgplt.set_cmap('hot')
863 plt.title("$v_1$")
864 plt.xlabel('$x_1$')
865 plt.ylabel('$x_3$')
866 plt.colorbar(orientation='horizontal')
867
868 plt.subplot(132)
869 imgplt = plt.imshow(self.v_f[:, y, :, 1].T, origin='lower')
870 imgplt.set_interpolation('nearest')
871 #imgplt.set_interpolation('bicubic')
872 #imgplt.set_cmap('hot')
873 plt.title("$v_2$")
874 plt.xlabel('$x_1$')
875 plt.ylabel('$x_3$')
876 plt.colorbar(orientation='horizontal')
877
878 plt.subplot(133)
879 imgplt = plt.imshow(self.v_f[:, y, :, 2].T, origin='lower')
880 imgplt.set_interpolation('nearest')
881 #imgplt.set_interpolation('bicubic')
882 #imgplt.set_cmap('hot')
883 plt.title("$v_3$")
884 plt.xlabel('$x_1$')
885 plt.ylabel('$x_3$')
886 plt.colorbar(orientation='horizontal')
887
888 filename = 'v_f-' + self.sid + '-y' + str(y) + '.' + graphics_format
889 plt.savefig(filename, transparent=False)
890 if verbose:
891 print('saved to ' + filename)
892 plt.clf()
893 plt.close()
894
895 def plotFluidVelocitiesZ(self, z=-1, graphics_format='png', verbose=True):
896 '''
897 Plot fluid velocities in a plane normal to the third axis.
898 The plot is saved in the current folder with the format
899 'v_f-<simulation id>-z<z value>.<graphics_format>'.
900
901 :param z: Plot velocities in fluid cells with these z axis values. If
902 this value is -1, the center z position is used.
903 :type z: int
904 :param graphics_format: Save the plot in this format
905 :type graphics_format: str
906 :param verbose: Print output filename after saving
907 :type verbose: bool
908
909 See also: :func:`writeFluidVTK()` and :func:`plotFluidVelocitiesY()`
910 '''
911 if not py_mpl:
912 print('Error: matplotlib module not found (plotFluidVelocitiesZ).')
913 return
914
915 if z == -1:
916 z = self.num[2]/2
917
918 plt.title("Fluid velocities")
919 plt.figure(figsize=[8, 8])
920
921 plt.subplot(131)
922 imgplt = plt.imshow(self.v_f[:, :, z, 0].T, origin='lower')
923 imgplt.set_interpolation('nearest')
924 #imgplt.set_interpolation('bicubic')
925 #imgplt.set_cmap('hot')
926 plt.title("$v_1$")
927 plt.xlabel('$x_1$')
928 plt.ylabel('$x_2$')
929 plt.colorbar(orientation='horizontal')
930
931 plt.subplot(132)
932 imgplt = plt.imshow(self.v_f[:, :, z, 1].T, origin='lower')
933 imgplt.set_interpolation('nearest')
934 #imgplt.set_interpolation('bicubic')
935 #imgplt.set_cmap('hot')
936 plt.title("$v_2$")
937 plt.xlabel('$x_1$')
938 plt.ylabel('$x_2$')
939 plt.colorbar(orientation='horizontal')
940
941 plt.subplot(133)
942 imgplt = plt.imshow(self.v_f[:, :, z, 2].T, origin='lower')
943 imgplt.set_interpolation('nearest')
944 #imgplt.set_interpolation('bicubic')
945 #imgplt.set_cmap('hot')
946 plt.title("$v_3$")
947 plt.xlabel('$x_1$')
948 plt.ylabel('$x_2$')
949 plt.colorbar(orientation='horizontal')
950
951 filename = 'v_f-' + self.sid + '-z' + str(z) + '.' + graphics_format
952 plt.savefig(filename, transparent=False)
953 if verbose:
954 print('saved to ' + filename)
955 plt.clf()
956 plt.close()
957
958 def plotFluidDiffAdvPresZ(self, graphics_format='png', verbose=True):
959 '''
960 Compare contributions to the velocity from diffusion and advection,
961 assuming the flow is 1D along the z-axis, phi=1, and dphi=0. This
962 solution is analog to the predicted velocity and not constrained by the
963 conservation of mass. The plot is saved in the output folder with the
964 name format '<simulation id>-diff_adv-t=<current time>s-mu=<dynamic
965 viscosity>Pa-s.<graphics_format>'.
966
967 :param graphics_format: Save the plot in this format
968 :type graphics_format: str
969 :param verbose: Print output filename after saving
970 :type verbose: bool
971 '''
972 if not py_mpl:
973 print('Error: matplotlib module not found (plotFluidDiffAdvPresZ).')
974 return
975
976 # The v_z values are read from self.v_f[0, 0, :, 2]
977 dz = self.L[2]/self.num[2]
978 rho = self.rho_f
979
980 # Central difference gradients
981 dvz_dz = (self.v_f[0, 0, 1:, 2] - self.v_f[0, 0, :-1, 2])/(2.0*dz)
982 dvzvz_dz = (self.v_f[0, 0, 1:, 2]**2 - self.v_f[0, 0, :-1, 2]**2)\
983 /(2.0*dz)
984
985 # Diffusive contribution to velocity change
986 dvz_diff = 2.0*self.mu/rho*dvz_dz*self.time_dt
987
988 # Advective contribution to velocity change
989 dvz_adv = dvzvz_dz*self.time_dt
990
991 # Pressure gradient
992 dp_dz = (self.p_f[0, 0, 1:] - self.p_f[0, 0, :-1])/(2.0*dz)
993
994 cellno = numpy.arange(1, self.num[2])
995
996 fig = plt.figure()
997 titlesize = 12
998
999 plt.subplot(1, 3, 1)
1000 plt.title('Pressure', fontsize=titlesize)
1001 plt.ylabel('$i_z$')
1002 plt.xlabel('$p_z$')
1003 plt.plot(self.p_f[0, 0, :], numpy.arange(self.num[2]))
1004 plt.grid()
1005
1006 plt.subplot(1, 3, 2)
1007 plt.title('Pressure gradient', fontsize=titlesize)
1008 plt.ylabel('$i_z$')
1009 plt.xlabel('$\Delta p_z$')
1010 plt.plot(dp_dz, cellno)
1011 plt.grid()
1012
1013 plt.subplot(1, 3, 3)
1014 plt.title('Velocity prediction terms', fontsize=titlesize)
1015 plt.ylabel('$i_z$')
1016 plt.xlabel('$\Delta v_z$')
1017 plt.plot(dvz_diff, cellno, label='Diffusion')
1018 plt.plot(dvz_adv, cellno, label='Advection')
1019 plt.plot(dvz_diff+dvz_adv, cellno, '--', label='Sum')
1020 leg = plt.legend(loc='best', prop={'size':8})
1021 leg.get_frame().set_alpha(0.5)
1022 plt.grid()
1023
1024 plt.tight_layout()
1025 filename = '../output/{}-diff_adv-t={:.2e}s-mu={:.2e}Pa-s.{}'\
1026 .format(self.sid, self.time_current[0], self.mu[0],
1027 graphics_format)
1028 plt.savefig(filename)
1029 if verbose:
1030 print('saved to ' + filename)
1031 plt.clf()
1032 plt.close(fig)
1033
1034 def plotLoadCurve(self, graphics_format='png', verbose=True):
1035 '''
1036 Plot the load curve (log time vs. upper wall movement). The plot is
1037 saved in the current folder with the file name
1038 '<simulation id>-loadcurve.<graphics_format>'.
1039 The consolidation coefficient calculations are done on the base of
1040 Bowles 1992, p. 129--139, using the "Casagrande" method.
1041 It is assumed that the consolidation has stopped at the end of the
1042 simulation (i.e. flat curve).
1043
1044 :param graphics_format: Save the plot in this format
1045 :type graphics_format: str
1046 :param verbose: Print output filename after saving
1047 :type verbose: bool
1048 '''
1049 if not py_mpl:
1050 print('Error: matplotlib module not found (plotLoadCurve).')
1051 return
1052
1053 t = numpy.empty(self.status())
1054 H = numpy.empty_like(t)
1055 from .core import sim
1056 sb = sim(self.sid, fluid=self.fluid)
1057 sb.readfirst(verbose=False)
1058 for i in numpy.arange(1, self.status()+1):
1059 sb.readstep(i, verbose=False)
1060 if i == 0:
1061 load = sb.w_sigma0[0]
1062 t[i-1] = sb.time_current[0]
1063 H[i-1] = sb.w_x[0]
1064
1065 # find consolidation parameters
1066 H0 = H[0]
1067 H100 = H[-1]
1068 H50 = (H0 + H100)/2.0
1069 T50 = 0.197 # case I
1070
1071 # find the time where 50% of the consolidation (H50) has happened by
1072 # linear interpolation. The values in H are expected to be
1073 # monotonically decreasing. See Numerical Recipies p. 115
1074 i_lower = 0
1075 i_upper = self.status()-1
1076 while i_upper - i_lower > 1:
1077 i_midpoint = int((i_upper + i_lower)/2)
1078 if H50 < H[i_midpoint]:
1079 i_lower = i_midpoint
1080 else:
1081 i_upper = i_midpoint
1082 t50 = t[i_lower] + (t[i_upper] - t[i_lower]) * \
1083 (H50 - H[i_lower])/(H[i_upper] - H[i_lower])
1084
1085 c_coeff = T50*H50**2.0/(t50)
1086 if self.fluid:
1087 e = numpy.mean(sb.phi[:, :, 3:-8]) # ignore boundaries
1088 else:
1089 e = sb.voidRatio()
1090
1091 phi_bar = e
1092 fig = plt.figure()
1093 plt.xlabel('Time [s]')
1094 plt.ylabel('Height [m]')
1095 plt.title('$c_v$=%.2e m$^2$ s$^{-1}$ at %.1f kPa and $e$=%.2f' \
1096 % (c_coeff, sb.w_sigma0[0]/1000.0, e))
1097 plt.semilogx(t, H, '+-')
1098 plt.axhline(y=H0, color='gray')
1099 plt.axhline(y=H50, color='gray')
1100 plt.axhline(y=H100, color='gray')
1101 plt.axvline(x=t50, color='red')
1102 plt.grid()
1103 filename = self.sid + '-loadcurve.' + graphics_format
1104 plt.savefig(filename)
1105 if verbose:
1106 print('saved to ' + filename)
1107 plt.clf()
1108 plt.close(fig)
1109
1110 def plotConvergence(self, graphics_format='png', verbose=True):
1111 '''
1112 Plot the convergence evolution in the CFD solver. The plot is saved
1113 in the output folder with the file name
1114 '<simulation id>-conv.<graphics_format>'.
1115
1116 :param graphics_format: Save the plot in this format
1117 :type graphics_format: str
1118 :param verbose: Print output filename after saving
1119 :type verbose: bool
1120
1121 See also: :func:`convergence()`
1122 '''
1123 if not py_mpl:
1124 print('Error: matplotlib module not found (plotConvergence).')
1125 return
1126
1127 fig = plt.figure()
1128 conv = self.convergence()
1129
1130 plt.title('Convergence evolution in CFD solver in "' + self.sid + '"')
1131 plt.xlabel('Time step')
1132 plt.ylabel('Jacobi iterations')
1133 plt.plot(conv[:, 0], conv[:, 1])
1134 plt.grid()
1135 filename = self.sid + '-conv.' + graphics_format
1136 plt.savefig(filename)
1137 if verbose:
1138 print('saved to ' + filename)
1139 plt.clf()
1140 plt.close(fig)
1141
1142 def plotSinFunction(self, baseval, A, f, phi=0.0, xlabel='$t$ [s]',
1143 ylabel='$y$', plotstyle='.', outformat='png',
1144 verbose=True):
1145 '''
1146 Plot the values of a sinusoidal modulated base value. Saves the output
1147 as a plot in the current folder.
1148 The time values will range from `self.time_current` to
1149 `self.time_total`.
1150
1151 :param baseval: The center value which the sinusoidal fluctuations are
1152 modulating
1153 :type baseval: float
1154 :param A: The fluctuation amplitude
1155 :type A: float
1156 :param phi: The phase shift [s]
1157 :type phi: float
1158 :param xlabel: The label for the x axis
1159 :type xlabel: str
1160 :param ylabel: The label for the y axis
1161 :type ylabel: str
1162 :param plotstyle: Matplotlib-string for specifying plotting style
1163 :type plotstyle: str
1164 :param outformat: File format of the output plot
1165 :type outformat: str
1166 :param verbose: Print output filename after saving
1167 :type verbose: bool
1168 '''
1169 if not py_mpl:
1170 print('Error: matplotlib module not found (plotSinFunction).')
1171 return
1172
1173 fig = plt.figure(figsize=[8, 6])
1174 steps_left = (self.time_total[0] - self.time_current[0]) \
1175 /self.time_file_dt[0]
1176 t = numpy.linspace(self.time_current[0], self.time_total[0], steps_left)
1177 f = baseval + A*numpy.sin(2.0*numpy.pi*f*t + phi)
1178 plt.plot(t, f, plotstyle)
1179 plt.grid()
1180 plt.xlabel(xlabel)
1181 plt.ylabel(ylabel)
1182 plt.tight_layout()
1183 filename = self.sid + '-sin.' + outformat
1184 plt.savefig(filename)
1185 if verbose:
1186 print(filename)
1187 plt.clf()
1188 plt.close(fig)