tshear-results-forces.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-forces.py (7512B)
---
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
21 sigma0 = float(sys.argv[1])
22 #c_grad_p = 1.0
23 c_grad_p = float(sys.argv[2])
24 c_phi = 1.0
25
26 #sid = 'shear-sigma0=' + str(sigma0) + '-c_phi=' + \
27 # str(c_phi) + '-c_grad_p=' + str(c_grad_p) + '-hi_mu-lo_visc'
28 sid = 'halfshear-sigma0=' + str(sigma0) + '-c=' + str(c_grad_p) + '-shear'
29 sim = sphere.sim(sid, fluid=True)
30 sim.readfirst(verbose=False)
31
32 # particle z positions
33 zpos_p = numpy.zeros((len(steps), sim.np))
34
35 # cell midpoint cell positions
36 zpos_c = numpy.zeros((len(steps), sim.num[2]))
37 dz = sim.L[2]/sim.num[2]
38 for i in numpy.arange(sim.num[2]):
39 zpos_c[:,i] = i*dz + 0.5*dz
40
41 # particle x displacements
42 xdisp = numpy.zeros((len(steps), sim.np))
43
44 # particle-fluid force per particle
45 f_pf = numpy.zeros_like(xdisp)
46
47 # pressure - hydrostatic pressure
48 #dev_p = numpy.zeros((len(steps), sim.num[2]))
49
50 # mean porosity
51 phi_bar = numpy.zeros((len(steps), sim.num[2]))
52
53 # mean porosity change
54 dphi_bar = numpy.zeros((len(steps), sim.num[2]))
55
56 # mean per-particle values
57 xdisp_mean = numpy.zeros((len(steps), sim.num[2]))
58 f_pf_mean = numpy.zeros((len(steps), sim.num[2]))
59
60 shear_strain = numpy.zeros(len(steps))
61
62 s = 0
63 for step_str in steps:
64
65 step = int(step_str)
66
67 if os.path.isfile('../output/' + sid + '.status.dat'):
68
69 for substep in numpy.arange(nsteps_avg):
70
71 if step + substep > sim.status():
72 raise Exception(
73 'Simulation step %d not available (sim.status = %d).'
74 % (step, sim.status()))
75
76 sim.readstep(step + substep, verbose=False)
77
78 zpos_p[s,:] += sim.x[:,2]/nsteps_avg
79
80 xdisp[s,:] += sim.xyzsum[:,0]/nsteps_avg
81
82 '''
83 for i in numpy.arange(sim.np):
84 f_pf[s,i] += \
85 sim.f_sum[i].dot(sim.f_sum[i])/nsteps_avg
86 '''
87 f_pf[s,:] += sim.f_sum[:,2]
88
89 #dev_p[s,:] += \
90 #numpy.average(numpy.average(sim.p_f, axis=0), axis=0)\
91 #/nsteps_avg
92
93 phi_bar[s,:] += \
94 numpy.average(numpy.average(sim.phi, axis=0), axis=0)\
95 /nsteps_avg
96
97 dphi_bar[s,:] += \
98 numpy.average(numpy.average(sim.dphi, axis=0), axis=0)\
99 /nsteps_avg/sim.time_dt
100
101 shear_strain[s] += sim.shearStrain()/nsteps_avg
102
103 # calculate mean values of xdisp and f_pf
104 for iz in numpy.arange(sim.num[2]):
105 z_bot = iz*dz
106 z_top = (iz+1)*dz
107 I = numpy.nonzero((zpos_p[s,:] >= z_bot) & (zpos_p[s,:] < z_top))
108 if len(I) > 0:
109 xdisp_mean[s,iz] = numpy.mean(xdisp[s,I])
110 f_pf_mean[s,iz] = numpy.mean(f_pf[s,I])
111
112 else:
113 print(sid + ' not found')
114 s += 1
115
116
117 #fig = plt.figure(figsize=(8,4*(len(steps))+1))
118 fig = plt.figure(figsize=(8,5*(len(steps))+1))
119
120 ax = []
121 for s in numpy.arange(len(steps)):
122
123
124 ax.append(plt.subplot(len(steps)*100 + 31 + s*3))
125 ax.append(plt.subplot(len(steps)*100 + 32 + s*3, sharey=ax[s*4+0]))
126 ax.append(plt.subplot(len(steps)*100 + 33 + s*3, sharey=ax[s*4+0]))
127 ax.append(ax[s*4+2].twiny())
128
129 ax[s*4+0].plot(xdisp[s], zpos_p[s], ',', color = '#888888')
130 ax[s*4+0].plot(xdisp_mean[s], zpos_c[s], color = 'k')
131
132 # remove particles with 0.0 pressure force
133 I = numpy.nonzero(numpy.abs(f_pf[s]) > .01)
134 f_pf_nonzero = f_pf[s][I]
135 zpos_p_nonzero = zpos_p[s][I]
136 I = numpy.nonzero(numpy.abs(f_pf_mean[s]) > .01)
137 f_pf_mean_nonzero = f_pf_mean[s][I]
138 zpos_c_nonzero = zpos_c[s][I]
139
140 #ax[s*4+1].plot(f_pf[s], zpos_p[s], ',', color = '#888888')
141 ax[s*4+1].plot(f_pf_nonzero, zpos_p_nonzero, ',', color = '#888888')
142 #ax[s*4+1].plot(f_pf_mean[s][1:-2], zpos_c[s][1:-2], color = 'k')
143 ax[s*4+1].plot(f_pf_mean_nonzero, zpos_c_nonzero, color = 'k')
144 #ax[s*4+1].plot([0.0, 0.0], [0.0, sim.L[2]], '--', color='k')
145
146 #ax[s*4+2].plot(dev_p[s]/1000.0, zpos_c[s], 'k')
147 #ax[s*4+2].plot(phi_bar[s,1:], zpos_c[s,1:], '-k', linewidth=3)
148 ax[s*4+2].plot(phi_bar[s,1:], zpos_c[s,1:], '-k')
149
150 #phicolor = '#888888'
151 #ax[s*4+3].plot(phi_bar[s], zpos_c[s], '-', color = phicolor)
152 #for tl in ax[s*4+3].get_xticklabels():
153 #tl.set_color(phicolor)
154 ax[s*4+3].plot(dphi_bar[s,1:], zpos_c[s,1:], '--k')
155 #ax[s*4+3].plot(dphi_bar[s,1:], zpos_c[s,1:], '-k', linewidth=3)
156 #ax[s*4+3].plot(dphi_bar[s,1:], zpos_c[s,1:], '-w', linewidth=2)
157
158 max_z = numpy.max(zpos_p)
159 ax[s*4+0].set_ylim([0, max_z])
160
161 #ax[s*4+1].set_xlim([0.15, 0.46]) # f_pf
162 ax[s*4+1].set_xlim([0.235, 0.409]) # f_pf
163 ax[s*4+1].set_ylim([0, max_z])
164
165 ax[s*4+2].set_ylim([0, max_z])
166 ax[s*4+2].set_xlim([0.33, 0.6]) # phi
167 ax[s*4+3].set_xlim([-0.09, 0.024]) # dphi/dt
168
169 #plt.plot(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
170 #plt.semilogx(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
171 #plt.semilogy(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
172 #plt.loglog(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
173
174 ax[s*4+0].set_ylabel('Vertical position $z$ [m]')
175 ax[s*4+0].set_xlabel('$\\boldsymbol{x}^x_\\text{p}$ [m]')
176 ax[s*4+1].set_xlabel('$\\boldsymbol{f}^z_\\text{pf}$ [N]')
177 #ax[s*4+2].set_xlabel('$\\bar{p_\\text{f}}$ [kPa]')
178 #ax[s*4+3].set_xlabel('$\\bar{\\phi}$ [-]', color=phicolor)
179 ax[s*4+2].set_xlabel('$\\bar{\\phi}$ [-] (solid)')
180 ax[s*4+3].set_xlabel('$\\delta \\bar{\\phi}/\\delta t$ [-] (dashed)')
181 plt.setp(ax[s*4+1].get_yticklabels(), visible=False)
182 plt.setp(ax[s*4+2].get_yticklabels(), visible=False)
183
184 ax[s*4+0].get_xaxis().set_major_locator(MaxNLocator(nbins=5))
185 ax[s*4+1].get_xaxis().set_major_locator(MaxNLocator(nbins=5))
186 ax[s*4+2].get_xaxis().set_major_locator(MaxNLocator(nbins=5))
187
188 plt.setp(ax[s*4+0].xaxis.get_majorticklabels(), rotation=90)
189 plt.setp(ax[s*4+1].xaxis.get_majorticklabels(), rotation=90)
190 plt.setp(ax[s*4+2].xaxis.get_majorticklabels(), rotation=90)
191 plt.setp(ax[s*4+3].xaxis.get_majorticklabels(), rotation=90)
192
193 #if s == 0:
194 #y = 0.95
195 #if s == 1:
196 #y = 0.55
197
198 strain_str = 'Shear strain $\\gamma = %.3f$' % (shear_strain[s])
199 #fig.text(0.1, y, strain_str, horizontalalignment='left', fontsize=22)
200 #ax[s*4+0].annotate(strain_str, xytext=(0,1.1), textcoords='figure fraction',
201 #horizontalalignment='left', fontsize=22)
202 plt.text(0.05, 1.06, strain_str, horizontalalignment='left', fontsize=22,
203 transform=ax[s*4+0].transAxes)
204 #ax[s*4+0].set_title(strain_str)
205
206 ax[s*4+0].grid()
207 ax[s*4+1].grid()
208 ax[s*4+2].grid()
209 #ax1.legend(loc='lower right', prop={'size':18})
210 #ax2.legend(loc='lower right', prop={'size':18})
211
212 plt.tight_layout()
213 plt.subplots_adjust(wspace = .05)
214 plt.MaxNLocator(nbins=4)
215
216 filename = 'shear-' + str(int(sigma0/1000.0)) + 'kPa-forces.pdf'
217 plt.savefig(filename)
218 shutil.copyfile(filename, '/home/adc/articles/own/2-org/' + filename)
219 print(filename)