tshear-results-velfac.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-velfac.py (9548B)
---
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 sys
11 import numpy
12 import sphere
13 from permeabilitycalculator import *
14 import matplotlib.pyplot as plt
15
16 smoothed_results = False
17 contact_forces = False
18 pressures = False
19 zflow = True
20
21 #sigma0_list = numpy.array([1.0e3, 2.0e3, 4.0e3, 10.0e3, 20.0e3, 40.0e3])
22 #sigma0 = 10.0e3
23 sigma0 = float(sys.argv[1])
24 #cvals = [1.0, 0.1]
25 #cvals = [1.0, 0.1, 0.01]
26 #cvals = [1.0]
27 fluid = False
28 if int(sys.argv[2]) == 1:
29 print('fluid = True')
30 fluid = True
31 c_v = float(sys.argv[3])
32 c_a = float(sys.argv[4])
33 velfacvals = [0.5, 1.0, 2.0]
34
35 # return a smoothed version of in. The returned array is smaller than the
36 # original input array
37 '''
38 def smooth(in_arr, plus_minus_steps):
39 out_arr = numpy.zeros(in_arr.size - 2*plus_minus_steps + 1)
40 s = 0
41 for i in numpy.arange(in_arr.size):
42 if i >= plus_minus_steps and i < plus_minus_steps:
43 for i in numpy.arange(-plus_minus_steps, plus_minus_steps+1):
44 out_arr[s] += in_arr[s+i]/(2.0*plus_minus_steps)
45 s += 1
46 '''
47
48 def smooth(x, window_len=10, window='hanning'):
49 """smooth the data using a window with requested size.
50
51 This method is based on the convolution of a scaled window with the signal.
52 The signal is prepared by introducing reflected copies of the signal
53 (with the window size) in both ends so that transient parts are minimized
54 in the begining and end part of the output signal.
55
56 input:
57 x: the input signal
58 window_len: the dimension of the smoothing window
59 window: the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'
60 flat window will produce a moving average smoothing.
61
62 output:
63 the smoothed signal
64
65 example:
66
67 import numpy as np
68 t = np.linspace(-2,2,0.1)
69 x = np.sin(t)+np.random.randn(len(t))*0.1
70 y = smooth(x)
71
72 see also:
73
74 numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve
75 scipy.signal.lfilter
76
77 TODO: the window parameter could be the window itself if an array instead of a string
78 """
79
80 if x.ndim != 1:
81 raise ValueError, "smooth only accepts 1 dimension arrays."
82
83 if x.size < window_len:
84 raise ValueError, "Input vector needs to be bigger than window size."
85
86 if window_len < 3:
87 return x
88
89 if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
90 raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"
91
92 s=numpy.r_[2*x[0]-x[window_len:1:-1], x, 2*x[-1]-x[-1:-window_len:-1]]
93 #print(len(s))
94
95 if window == 'flat': #moving average
96 w = numpy.ones(window_len,'d')
97 else:
98 w = getattr(numpy, window)(window_len)
99 y = numpy.convolve(w/w.sum(), s, mode='same')
100 return y[window_len-1:-window_len+1]
101
102
103 smooth_window = 10
104
105 shear_strain = [[], [], [], []]
106 friction = [[], [], [], []]
107 if smoothed_results:
108 friction_smooth = [[], [], [], []]
109 dilation = [[], [], [], []]
110 p_min = [[], [], [], []]
111 p_mean = [[], [], [], []]
112 p_max = [[], [], [], []]
113 f_n_mean = [[], [], [], []]
114 f_n_max = [[], [], [], []]
115 v_f_z_mean = [[], [], [], []]
116
117 for c in numpy.arange(0,len(velfacvals)):
118
119 velfac = velfacvals[c]
120
121 #sid = 'shear-sigma0=' + str(sigma0) + '-c_phi=' + \
122 #str(c_phi) + '-c_grad_p=' + str(c_grad_p) + \
123 #'-hi_mu-lo_visc-hw'
124 if fluid:
125 sid = 'halfshear-sigma0=' + str(sigma0) + '-c_v=' + str(c_v) + \
126 '-c_a=' + str(c_a) + '-velfac=' + str(velfac) + '-shear'
127 else:
128 sid = 'halfshear-sigma0=' + str(sigma0) + '-velfac=' + str(velfac) + \
129 '-shear'
130 if os.path.isfile('../output/' + sid + '.status.dat'):
131
132 sim = sphere.sim(sid, fluid=fluid)
133 shear_strain[c] = numpy.zeros(sim.status())
134 friction[c] = numpy.zeros_like(shear_strain[c])
135 dilation[c] = numpy.zeros_like(shear_strain[c])
136 if smoothed_results:
137 friction_smooth[c] = numpy.zeros_like(shear_strain[c])
138
139 sim.readlast(verbose=False)
140 sim.visualize('shear')
141 shear_strain[c] = sim.shear_strain
142 #shear_strain[c] = numpy.arange(sim.status()+1)
143 friction[c] = sim.tau/sim.sigma_eff
144 dilation[c] = sim.dilation
145 if smoothed_results:
146 friction_smooth[c] = smooth(friction[c], smooth_window)
147
148 # fluid pressures and particle forces
149 if pressures or contact_forces:
150 p_mean[c] = numpy.zeros_like(shear_strain[c])
151 p_min[c] = numpy.zeros_like(shear_strain[c])
152 p_max[c] = numpy.zeros_like(shear_strain[c])
153 f_n_mean[c] = numpy.zeros_like(shear_strain[c])
154 f_n_max[c] = numpy.zeros_like(shear_strain[c])
155 for i in numpy.arange(sim.status()):
156 if pressures and fluid:
157 sim.readstep(i, verbose=False)
158 iz_top = int(sim.w_x[0]/(sim.L[2]/sim.num[2]))-1
159 p_mean[c][i] = numpy.mean(sim.p_f[:,:,0:iz_top])/1000
160 p_min[c][i] = numpy.min(sim.p_f[:,:,0:iz_top])/1000
161 p_max[c][i] = numpy.max(sim.p_f[:,:,0:iz_top])/1000
162
163 if contact_forces:
164 sim.findNormalForces()
165 f_n_mean[c][i] = numpy.mean(sim.f_n_magn)
166 f_n_max[c][i] = numpy.max(sim.f_n_magn)
167
168 if zflow and fluid:
169 v_f_z_mean[c] = numpy.zeros_like(shear_strain[c])
170 for i in numpy.arange(sim.status()):
171 v_f_z_mean[c][i] = numpy.mean(sim.v_f[:,:,:,2])
172
173 else:
174 print(sid + ' not found')
175
176 # produce VTK files
177 #for sid in sids:
178 #sim = sphere.sim(sid, fluid=True)
179 #sim.writeVTKall()
180 c += 1
181
182
183 if zflow and fluid:
184 fig = plt.figure(figsize=(8,10))
185 else:
186 fig = plt.figure(figsize=(8,8)) # (w,h)
187 #fig = plt.figure(figsize=(8,12))
188 #fig = plt.figure(figsize=(8,16))
189 fig.subplots_adjust(hspace=0.0)
190
191 #plt.subplot(3,1,1)
192 #plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
193
194 if zflow and fluid:
195 ax1 = plt.subplot(311)
196 ax2 = plt.subplot(312, sharex=ax1)
197 ax3 = plt.subplot(313, sharex=ax1)
198 else:
199 ax1 = plt.subplot(211)
200 ax2 = plt.subplot(212, sharex=ax1)
201 #ax3 = plt.subplot(413, sharex=ax1)
202 #ax4 = plt.subplot(414, sharex=ax1)
203
204 color = ['b','g','r','c']
205 for c in numpy.arange(len(velfacvals)):
206
207 if smoothed_results:
208 ax1.plot(shear_strain[c][1:], friction_smooth[c][1:], \
209 label='velfac = %.1f' % (velfacvals[c]), linewidth=1, alpha=0.5)
210 else:
211 ax1.plot(shear_strain[c][1:], friction[c][1:], \
212 label='velfac = %.1f' % (velfacvals[c]), linewidth=1, alpha=0.5)
213
214 ax2.plot(shear_strain[c][1:], dilation[c][1:], \
215 label='velfac = %.1f' % (velfacvals[c]), linewidth=1, alpha=0.5)
216
217 if zflow and fluid:
218 ax3.plot(shear_strain[c][1:], v_f_z_mean[c][1:],
219 label='velfac = %.1f' % (velfacvals[c]), linewidth=1, alpha=0.5)
220
221
222 '''
223 alpha = 0.5
224 ax3.plot(shear_strain[c][1:], p_max[c][1:], '-' + color[c], alpha=alpha)
225 ax3.plot(shear_strain[c][1:], p_mean[c][1:], '-' + color[c], \
226 label='$c$ = %.2f' % (cvals[c-1]), linewidth=2)
227 ax3.plot(shear_strain[c][1:], p_min[c][1:], '-' + color[c], alpha=alpha)
228
229 ax3.fill_between(shear_strain[c][1:], p_min[c][1:], p_max[c][1:],
230 where=p_min[c][1:]<=p_max[c][1:], facecolor=color[c],
231 interpolate=True, alpha=alpha)
232
233 ax4.plot(shear_strain[c][1:], f_n_mean[c][1:], '-' + color[c],
234 label='$c$ = %.2f' % (cvals[c-1]), linewidth=2)
235 ax4.plot(shear_strain[c][1:], f_n_max[c][1:], '--' + color[c])
236 #label='$c$ = %.2f' % (cvals[c-1]), linewidth=2)
237 '''
238
239 #ax4.set_xlabel('Shear strain $\\gamma$ [-]')
240
241 ax1.set_ylabel('Shear friction $\\tau/\\sigma\'$ [-]')
242 ax2.set_ylabel('Dilation $\\Delta h/(2r)$ [-]')
243 if zflow and fluid:
244 ax3.set_ylabel('$\\boldsymbol{v}_\\text{f}^z h$ [ms$^{-1}$]')
245 #ax3.set_ylabel('Fluid pressure $p_\\text{f}$ [kPa]')
246 #ax4.set_ylabel('Particle contact force $||\\boldsymbol{f}_\\text{p}||$ [N]')
247
248 #ax1.set_xlim([200,300])
249 #ax3.set_ylim([595,608])
250
251 plt.setp(ax1.get_xticklabels(), visible=False)
252 #plt.setp(ax2.get_xticklabels(), visible=False)
253 #plt.setp(ax3.get_xticklabels(), visible=False)
254
255 ax1.grid()
256 ax2.grid()
257 if zflow and fluid:
258 ax3.grid()
259 #ax4.grid()
260
261 legend_alpha=0.5
262 ax1.legend(loc='lower right', prop={'size':18}, fancybox=True,
263 framealpha=legend_alpha)
264 ax2.legend(loc='lower right', prop={'size':18}, fancybox=True,
265 framealpha=legend_alpha)
266 if zflow and fluid:
267 ax3.legend(loc='lower right', prop={'size':18}, fancybox=True,
268 framealpha=legend_alpha)
269 #ax4.legend(loc='best', prop={'size':18}, fancybox=True,
270 #framealpha=legend_alpha)
271
272 plt.tight_layout()
273 plt.subplots_adjust(hspace=0.05)
274 filename = 'shear-' + str(int(sigma0/1000.0)) + 'kPa-stress-dilation-velfac-dry.pdf'
275 if fluid:
276 filename = 'shear-' + str(int(sigma0/1000.0)) + \
277 'kPa-stress-dilation-velfac-c_v=' + str(c_v) + \
278 '-c_a=' + str(c_a) + '.pdf'
279 #print(os.getcwd() + '/' + filename)
280 plt.savefig(filename)
281 shutil.copyfile(filename, '/home/adc/articles/own/2-org/' + filename)
282 print(filename)