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