thalfshear-darcy-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
---
thalfshear-darcy-internals.py (9695B)
---
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 import subprocess
16
17 import seaborn as sns
18 #sns.set(style='ticks', palette='Set2')
19 #sns.set(style='ticks', palette='colorblind')
20 #sns.set(style='ticks', palette='Set2')
21 sns.set(style='ticks', palette='Blues')
22
23 #sigma0 = float(sys.argv[1])
24 sigma0_list = [20000.0, 80000.0]
25 #k_c = 3.5e-13
26 #k_c = float(sys.argv[1])
27 k_c_list = [3.5e-13, 3.5e-14, 3.5e-15]
28
29 #if k_c == 3.5e-15:
30 # steps = [1232, 1332, 1433, 1534, 1635]
31 #elif k_c == 3.5e-13:
32 # steps = [100, 200, 300, 410, 515]
33 #else:
34 # steps = [10, 50, 100, 1000, 1999]
35 nsteps_avg = 1 # no. of steps to average over
36 #nsteps_avg = 100 # no. of steps to average over
37
38 steps = [1, 51, 101, 152, 204]
39
40 for sigma0 in sigma0_list:
41 for k_c in k_c_list:
42
43 sid = 'halfshear-darcy-sigma0=' + str(sigma0) + '-k_c=' + str(k_c) + \
44 '-mu=1.797e-06-velfac=1.0-shear'
45 sim = sphere.sim(sid, fluid=True)
46 sim.readfirst(verbose=False)
47
48 # particle z positions
49 zpos_p = numpy.zeros((len(steps), sim.np))
50
51 # cell midpoint cell positions
52 zpos_c = numpy.zeros((len(steps), sim.num[2]))
53 dz = sim.L[2]/sim.num[2]
54 for i in numpy.arange(sim.num[2]):
55 zpos_c[:,i] = i*dz + 0.5*dz
56
57 # particle x displacements
58 xdisp = numpy.zeros((len(steps), sim.np))
59
60 # particle z velocity
61 v_z_p = numpy.zeros((len(steps), sim.np))
62
63 # fluid permeability
64 k = numpy.zeros((len(steps), sim.num[0], sim.num[1], sim.num[2]))
65 k_bar = numpy.zeros((len(steps), sim.num[2]))
66
67 # pressure
68 p = numpy.zeros((len(steps), sim.num[2]))
69
70 # mean per-particle values
71 v_z_p_bar = numpy.zeros((len(steps), sim.num[2]))
72 v_z_f_bar = numpy.zeros((len(steps), sim.num[2]))
73
74 # particle-fluid force per particle
75 f_pf = numpy.zeros_like(xdisp)
76
77 # pressure - hydrostatic pressure
78 #dev_p = numpy.zeros((len(steps), sim.num[2]))
79
80 # mean porosity
81 phi_bar = numpy.zeros((len(steps), sim.num[2]))
82
83 # mean porosity change
84 dphi_bar = numpy.zeros((len(steps), sim.num[2]))
85
86 # mean per-particle values
87 xdisp_mean = numpy.zeros((len(steps), sim.num[2]))
88 f_pf_mean = numpy.zeros((len(steps), sim.num[2]))
89
90 shear_strain_start = numpy.zeros(len(steps))
91 shear_strain_end = numpy.zeros(len(steps))
92
93 #fig = plt.figure(figsize=(8,4*(len(steps))+1))
94 #fig = plt.figure(figsize=(8,4.5))
95 fig = plt.figure(figsize=(3.74*2,3.00))
96 ax = []
97 n = 4
98 ax.append(plt.subplot(1, n, 1)) # 0: xdisp
99 ax.append(plt.subplot(1, n, 2, sharey=ax[0])) # 3: k
100 ax.append(plt.subplot(1, n, 3, sharey=ax[0])) # 5: p_f
101 ax.append(plt.subplot(1, n, 4, sharey=ax[0])) # 6: f_pf_z
102
103 s = 0
104 for step_str in steps:
105
106 step = int(step_str)
107
108 if os.path.isfile('../output/' + sid + '.status.dat'):
109
110 for substep in numpy.arange(nsteps_avg):
111
112 if step + substep > sim.status():
113 raise Exception(
114 'Simulation step %d not available (sim.status = %d).'
115 % (step + substep, sim.status()))
116
117 sim.readstep(step + substep, verbose=False)
118
119 zpos_p[s,:] += sim.x[:,2]/nsteps_avg
120
121 xdisp[s,:] += sim.xyzsum[:,0]/nsteps_avg
122
123 '''
124 for i in numpy.arange(sim.np):
125 f_pf[s,i] += \
126 sim.f_sum[i].dot(sim.f_sum[i])/nsteps_avg
127 '''
128 f_pf[s,:] += sim.f_p[:,2]
129
130 dz = sim.L[2]/sim.num[2]
131 wall0_iz = int(sim.w_x[0]/dz)
132
133 p[s,:] += numpy.average(numpy.average(sim.p_f[:,:,:], axis=0),\
134 axis=0)/nsteps_avg
135
136 sim.findPermeabilities()
137 k[s,:] += sim.k[:,:,:]/nsteps_avg
138
139 k_bar[s,:] += \
140 numpy.average(numpy.average(sim.k[:,:,:], axis=0), axis=0)\
141 /nsteps_avg
142
143 if substep == 0:
144 shear_strain_start[s] = sim.shearStrain()
145 else:
146 shear_strain_end[s] = sim.shearStrain()
147
148 # calculate mean values of xdisp and f_pf
149 for iz in numpy.arange(sim.num[2]):
150 z_bot = iz*dz
151 z_top = (iz+1)*dz
152 I = numpy.nonzero((zpos_p[s,:] >= z_bot) & (zpos_p[s,:] < z_top))
153 if len(I) > 0:
154 xdisp_mean[s,iz] = numpy.mean(xdisp[s,I])
155 f_pf_mean[s,iz] = numpy.mean(f_pf[s,I])
156
157 k_bar[s][0] = k_bar[s][1]
158
159
160
161 #ax[0].plot(xdisp[s], zpos_p[s], ',', color = '#888888')
162 ax[0].plot(xdisp_mean[s], zpos_c[s], label='$\gamma$ = %.3f' %
163 (shear_strain_start[s]))
164
165 ax[1].semilogx(k_bar[s], zpos_c[s], label='$\gamma$ = %.3f' %
166 (shear_strain_start[s]))
167
168 ax[2].plot(p[s]/1000.0, zpos_c[s], label='$\gamma$ = %.3f' %
169 (shear_strain_start[s]))
170
171 # remove particles with 0.0 pressure force
172 I = numpy.nonzero(numpy.abs(f_pf[s]) > .01)
173 f_pf_nonzero = f_pf[s][I]
174 zpos_p_nonzero = zpos_p[s][I]
175 I = numpy.nonzero(numpy.abs(f_pf_mean[s]) > .01)
176 f_pf_mean_nonzero = f_pf_mean[s][I]
177 zpos_c_nonzero = zpos_c[s][I]
178 #f_pf_mean_nonzero = numpy.append(f_pf_mean_nonzero, 0.0)
179 #zpos_c_nonzero = numpy.append(zpos_c_nonzero, zpos_c[s][zpos_c_nonzero.size])
180
181 #ax[3].plot(f_pf_nonzero, zpos_p_nonzero, ',', alpha=0.5,
182 #color='#888888')
183 #ax[3].plot(f_pf_mean_nonzero, zpos_c_nonzero, label='$\gamma$ = %.3f' %
184 #(shear_strain_start[s]))
185 ax[3].plot(f_pf_mean_nonzero/(4.0*numpy.pi*sim.radius[0]**2)/1e3,
186 zpos_c_nonzero, label='$\gamma$ = %.3f' %
187 (shear_strain_start[s]))
188
189 else:
190 print(sid + ' not found')
191 s += 1
192
193
194 #max_z = numpy.max(zpos_p)
195 max_z = 0.5
196 #max_z = numpy.max(zpos_c[-2])
197 ax[0].set_ylim([0, max_z])
198 #ax[0].set_xlim([0, 0.5])
199 ax[0].set_xlim([0, 0.05])
200
201 if k_c == 3.5e-15:
202 #ax[1].set_xlim([1e-14, 1e-12])
203 ax[1].set_xlim([1e-16, 1e-14])
204 elif k_c == 3.5e-14:
205 ax[1].set_xlim([1e-15, 1e-13])
206 elif k_c == 3.5e-13:
207 #ax[1].set_xlim([1e-12, 1e-10])
208 ax[1].set_xlim([1e-14, 1e-12])
209
210 ax[0].set_ylabel('Vertical position $z$ [m]')
211 ax[0].set_xlabel('$\\bar{\\boldsymbol{x}}^x_\\text{p}$ [m]')
212 ax[1].set_xlabel('$\\bar{k}$ [m$^{2}$]')
213 ax[2].set_xlabel('$\\bar{p_\\text{f}}$ [kPa]')
214 #ax[3].set_xlabel('$\\boldsymbol{f}^z_\\text{i}$ [N]')
215 ax[3].set_xlabel('$\\bar{\boldsymbol{\sigma}}^z_\\text{i}$ [kPa]')
216
217 # align x labels
218 #labely = -0.3
219 #ax[0].xaxis.set_label_coords(0.5, labely)
220 #ax[1].xaxis.set_label_coords(0.5, labely)
221 #ax[2].xaxis.set_label_coords(0.5, labely)
222 #ax[3].xaxis.set_label_coords(0.5, labely)
223
224 plt.setp(ax[1].get_yticklabels(), visible=False)
225 plt.setp(ax[2].get_yticklabels(), visible=False)
226 plt.setp(ax[3].get_yticklabels(), visible=False)
227
228 plt.setp(ax[0].xaxis.get_majorticklabels(), rotation=90)
229 plt.setp(ax[1].xaxis.get_majorticklabels(), rotation=90)
230 plt.setp(ax[2].xaxis.get_majorticklabels(), rotation=90)
231 plt.setp(ax[3].xaxis.get_majorticklabels(), rotation=90)
232
233 '''
234 ax[0].grid()
235 ax[1].grid()
236 ax[2].grid()
237 ax[3].grid()
238 '''
239
240 sns.despine() # remove chartjunk
241
242 for i in range(4):
243 # vertical grid lines
244 ax[i].get_xaxis().grid(True, linestyle=':', linewidth=0.5)
245 # horizontal grid lines
246 ax[i].get_yaxis().grid(True, linestyle=':', linewidth=0.5)
247
248 if i>0:
249 ax[i].spines['left'].set_visible(False)
250 ax[i].get_yaxis().set_ticks_position('none')
251
252 legend_alpha=0.5
253 #ax[0].legend(loc='lower center', prop={'size':12}, fancybox=True,
254 #framealpha=legend_alpha)
255 ax[0].legend(loc='lower right')
256
257 #plt.subplots_adjust(wspace = .05) # doesn't work with tight_layout()
258 #plt.MaxNLocator(nbins=1) # doesn't work?
259 ax[0].locator_params(nbins=4)
260 ax[2].locator_params(nbins=5)
261 ax[3].locator_params(nbins=5)
262
263 plt.tight_layout()
264
265 filename = 'halfshear-darcy-internals-k_c=%.0e.pdf' % (k_c)
266 if sigma0 == 80000.0:
267 filename = 'halfshear-darcy-internals-k_c=%.0e-N80.pdf' % (k_c)
268
269 plt.savefig(filename)
270 #subprocess.call('pdfcrop ' + filename, shell=True)
271 shutil.copyfile(filename, '/home/adc/articles/own/2/graphics/' + filename)
272 print(filename)