tshear-results-pressures.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-pressures.py (5006B)
---
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 matplotlib.rcParams['image.cmap'] = 'bwr'
17
18 sigma0 = float(sys.argv[1])
19 #c_grad_p = 1.0
20 c_grad_p = [1.0, 0.1]
21 #c_grad_p = [1.0, 0.1, 0.01, 1e-07]
22 c_phi = 1.0
23
24
25 #sid = 'shear-sigma0=' + str(sigma0) + '-c_phi=' + \
26 # str(c_phi) + '-c_grad_p=' + str(c_grad_p) + '-hi_mu-lo_visc'
27 sid = 'halfshear-sigma0=' + str(sigma0) + '-c=' + str(c_grad_p[0]) + '-shear'
28 sim = sphere.sim(sid, fluid=True)
29 sim.readfirst(verbose=False)
30
31 # cell midpoint cell positions
32 zpos_c = numpy.zeros((len(c_grad_p), sim.num[2]))
33 dz = sim.L[2]/sim.num[2]
34 for c in numpy.arange(len(c_grad_p)):
35 for i in numpy.arange(sim.num[2]):
36 zpos_c[c,i] = i*dz + 0.5*dz
37
38
39 shear_strain = [[], [], [], []]
40 dev_pres = [[], [], [], []]
41 pres_static = [[], [], [], []]
42 pres = [[], [], [], []]
43
44 for c in numpy.arange(len(c_grad_p)):
45 sim.sid = 'halfshear-sigma0=' + str(sigma0) + '-c=' + str(c_grad_p[c]) \
46 + '-shear'
47
48 shear_strain[c] = numpy.zeros(sim.status())
49 dev_pres[c] = numpy.zeros((sim.num[2], sim.status()))
50 pres_static[c] = numpy.ones_like(dev_pres[c])*sim.p_f[0,0,-1]
51 pres[c] = numpy.zeros_like(dev_pres[c])
52
53 for i in numpy.arange(sim.status()):
54
55 sim.readstep(i, verbose=False)
56
57 pres[c][:,i] = numpy.average(numpy.average(sim.p_f, axis=0), axis=0)
58
59 dz = sim.L[2]/sim.num[2]
60 wall0_iz = int(sim.w_x[0]/dz)
61 for z in numpy.arange(0, wall0_iz+1):
62 pres_static[c][z,i] = \
63 (wall0_iz*dz - zpos_c[c,z] + 0.5*dz)\
64 *sim.rho_f*numpy.abs(sim.g[2])\
65 + sim.p_f[0,0,-1]
66
67 shear_strain[c][i] = sim.shearStrain()
68
69 dev_pres[c] = pres[c] - pres_static[c]
70
71 #fig = plt.figure(figsize=(8,6))
72 #fig = plt.figure(figsize=(8,12))
73 fig = plt.figure(figsize=(8,5*len(c_grad_p)+2))
74
75
76 #cmap = matplotlib.colors.ListedColormap(['b', 'w', 'r'])
77 #bounds = [min_p, 0, max_p]
78 #norm = matplotlib.colors.BoundaryNorm(bounds, cmap.N)
79
80 ax = []
81 for c in numpy.arange(len(c_grad_p)):
82
83 ax.append(plt.subplot(len(c_grad_p), 1, c+1))
84
85 max_p_dev = numpy.max((numpy.abs(numpy.min(dev_pres[c])),
86 numpy.max(dev_pres[c])))
87 #max_p = numpy.min(dev_pres)
88 min_p = -max_p_dev/1000.0
89 max_p = max_p_dev/1000.0
90 #min_p = -5.0
91 #max_p = 5.0
92
93 im1 = ax[c].pcolormesh(shear_strain[c], zpos_c[c], dev_pres[c]/1000.0,
94 vmin=min_p, vmax=max_p, rasterized=True)
95 #im1 = ax[c].pcolormesh(shear_strain[c], zpos_c[c], dev_pres[c]/1000.0,
96 #rasterized=True)
97 #im1 = ax[c].pcolormesh(shear_strain[c], zpos_c[c], pres[c]/1000.0,
98 #rasterized=True)
99 #if c == 0:
100 ax[c].set_xlim([0, numpy.max(shear_strain[c])])
101 ax[c].set_ylim([zpos_c[0,0], sim.w_x[0]])
102 ax[c].set_xlabel('Shear strain $\\gamma$ [-]')
103 ax[c].set_ylabel('Vertical position $z$ [m]')
104
105 #plt.text(0.0, 0.15, '$c = %.2f$' % (c_grad_p[c]),\
106 # horizontalalignment='left', fontsize=22,
107 # transform=ax[c].transAxes)
108 ax[c].set_title('$c = %.2f$' % (c_grad_p[c]))
109
110 #cb = plt.colorbar(im1, orientation='horizontal')
111 cb = plt.colorbar(im1)
112 cb.set_label('$p_\\text{f} - p^\\text{hyd}_\\text{f}$ [kPa]')
113 cb.solids.set_rasterized(True)
114
115 # annotate plot
116 #ax1.text(0.02, 0.15, 'compressive',
117 #bbox={'facecolor':'white', 'alpha':0.5, 'pad':10})
118
119 #ax1.text(0.12, 0.25, 'dilative',
120 #bbox={'facecolor':'white', 'alpha':0.5, 'pad':10})
121
122 #cb = plt.colorbar(im1, orientation='horizontal')
123 #cb.set_label('$p_\\text{f} - p^\\text{hyd}_\\text{f}$ [kPa]')
124 #cb.solids.set_rasterized(True)
125 '''
126 ax2 = plt.subplot(312)
127 im2 = ax2.pcolormesh(shear_strain, zpos_c, pres/1000.0, rasterized=True)
128 #ax2.set_xlim([0, shear_strain[-1]])
129 #ax2.set_ylim([zpos_c[0], sim.w_x[0]])
130 ax2.set_xlabel('Shear strain $\\gamma$ [-]')
131 ax2.set_ylabel('Vertical position $z$ [m]')
132 cb2 = plt.colorbar(im2)
133 cb2.set_label('Pressure $p_\\text{f}$ [kPa]')
134 cb2.solids.set_rasterized(True)
135
136 ax3 = plt.subplot(313)
137 im3 = ax3.pcolormesh(shear_strain, zpos_c, pres_static/1000.0, rasterized=True)
138 #ax3.set_xlim([0, shear_strain[-1]])
139 #ax3.set_ylim([zpos_c[0], sim.w_x[0]])
140 ax3.set_xlabel('Shear strain $\\gamma$ [-]')
141 ax3.set_ylabel('Vertical position $z$ [m]')
142 cb3 = plt.colorbar(im3)
143 cb3.set_label('Static Pressure $p_\\text{f}$ [kPa]')
144 cb3.solids.set_rasterized(True)
145 '''
146
147
148 #plt.MaxNLocator(nbins=4)
149 plt.tight_layout()
150 plt.subplots_adjust(wspace = .05)
151 #plt.MaxNLocator(nbins=4)
152
153 filename = 'shear-' + str(int(sigma0/1000.0)) + 'kPa-pressures.pdf'
154 plt.savefig(filename)
155 shutil.copyfile(filename, '/home/adc/articles/own/2-org/' + filename)
156 print(filename)