tvisualize-rs0.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
---
tvisualize-rs0.py (15103B)
---
1 #!/usr/bin/env python
2 import sphere
3 import numpy
4 import matplotlib
5 matplotlib.use('Agg')
6
7 import matplotlib.pyplot as plt
8 import matplotlib.patches
9 import matplotlib.colors
10 #import seaborn
11
12 jobname_prefix = 'rs0-'
13
14 verbose = True
15
16 # Visualization settings
17 outformat = 'pdf'
18 plotforcechains = False
19 calculateforcechains = False
20 # calculateforcechains = True
21 calculateforcechainhistory = False
22 legend_alpha = 0.7
23 linewidth = 0.5
24 smooth_friction = False
25 smooth_window = 10
26
27
28 # Simulation parameter values
29 effective_stresses = [10e3, 20e3, 100e3, 200e3, 1000e3, 2000e3]
30 velfacs = [0.1, 1.0, 10.0]
31 mu_s_vals = [0.5]
32 mu_d_vals = [0.5]
33
34 # return a smoothed version of in. The returned array is smaller than the
35 # original input array
36 def smooth(x, window_len=10, window='hanning'):
37 """smooth the data using a window with requested size.
38
39 This method is based on the convolution of a scaled window with the signal.
40 The signal is prepared by introducing reflected copies of the signal
41 (with the window size) in both ends so that transient parts are minimized
42 in the begining and end part of the output signal.
43
44 input:
45 x: the input signal
46 window_len: the dimension of the smoothing window
47 window: the type of window from 'flat', 'hanning', 'hamming',
48 'bartlett', 'blackman' flat window will produce a moving average
49 smoothing.
50
51 output:
52 the smoothed signal
53
54 example:
55
56 import numpy as np
57 t = np.linspace(-2,2,0.1)
58 x = np.sin(t)+np.random.randn(len(t))*0.1
59 y = smooth(x)
60
61 see also:
62
63 numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman,
64 numpy.convolve scipy.signal.lfilter
65
66 TODO: the window parameter could be the window itself if an array instead
67 of a string
68 """
69
70 if x.ndim != 1:
71 raise ValueError("smooth only accepts 1 dimension arrays.")
72
73 if x.size < window_len:
74 raise ValueError("Input vector needs to be bigger than window size.")
75
76 if window_len < 3:
77 return x
78
79 if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
80 raise ValueError(
81 "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'")
82
83 s = numpy.r_[2*x[0]-x[window_len:1:-1], x, 2*x[-1]-x[-1:-window_len:-1]]
84 # print(len(s))
85
86 if window == 'flat': # moving average
87 w = numpy.ones(window_len, 'd')
88 else:
89 w = getattr(numpy, window)(window_len)
90 y = numpy.convolve(w/w.sum(), s, mode='same')
91 return y[window_len-1:-window_len+1]
92
93
94 def rateStatePlot(sid):
95
96 matplotlib.rcParams.update({'font.size': 7, 'font.family': 'sans-serif'})
97 matplotlib.rc('text', usetex=True)
98 matplotlib.rcParams['text.latex.preamble'] = [r"\usepackage{amsmath}"]
99
100
101 rasterized = True # rasterize colored areas (pcolormesh and colorbar)
102 # izfactor = 4 # factor for vertical discretization in xvel
103 izfactor = 1 # factor for vertical discretization in poros
104
105
106 #############
107 # DATA READ #
108 #############
109
110 sim = sphere.sim(sid, fluid=False)
111 sim.readfirst()
112
113 # nsteps = 2
114 # nsteps = 10
115 # nsteps = 400
116 nsteps = sim.status()
117
118 t = numpy.empty(nsteps)
119
120 # Stress, pressure and friction
121 sigma_def = numpy.empty_like(t)
122 sigma_eff = numpy.empty_like(t)
123 tau_def = numpy.empty_like(t)
124 tau_eff = numpy.empty_like(t)
125 p_f_bar = numpy.empty_like(t)
126 p_f_top = numpy.empty_like(t)
127 mu = numpy.empty_like(t) # shear friction
128
129 # shear velocity plot
130 v = numpy.empty_like(t) # velocity
131 I = numpy.empty_like(t) # inertia number
132
133 # displacement and mean porosity plot
134 xdisp = numpy.empty_like(t)
135 phi_bar = numpy.empty_like(t)
136
137 # mean horizontal porosity plot
138 poros = numpy.empty((sim.num[2], nsteps))
139 xvel = numpy.zeros((sim.num[2]*izfactor, nsteps))
140 zpos_c = numpy.empty(sim.num[2]*izfactor)
141 dz = sim.L[2]/(sim.num[2]*izfactor)
142 for i in numpy.arange(sim.num[2]*izfactor):
143 zpos_c[i] = i*dz + 0.5*dz
144
145 # Contact statistics plot
146 n = numpy.empty(nsteps)
147 coordinationnumber = numpy.empty(nsteps)
148 nkept = numpy.empty(nsteps)
149
150
151 for i in numpy.arange(nsteps):
152
153 sim.readstep(i+1, verbose=verbose) # use step 1 to n
154
155 t[i] = sim.currentTime()
156
157 sigma_def[i] = sim.currentNormalStress('defined')
158 sigma_eff[i] = sim.currentNormalStress('effective')
159 tau_def[i] = sim.shearStress('defined')
160 tau_eff[i] = sim.shearStress('effective')
161 mu[i] = tau_eff[i]/sigma_eff[i]
162 #mu[i] = tau_eff[i]/sigma_def[i]
163
164 I[i] = sim.inertiaParameterPlanarShear()
165
166 v[i] = sim.shearVelocity()
167 xdisp[i] = sim.shearDisplacement()
168
169 #poros[:, i] = numpy.average(numpy.average(sim.phi, axis=0), axis=0)
170
171 # calculate mean values of xvel
172 dz = sim.L[2]/(sim.num[2]*izfactor)
173 for iz in numpy.arange(sim.num[2]*izfactor):
174 z_bot = iz*dz
175 z_top = (iz+1)*dz
176 idx = numpy.nonzero((sim.x[:, 2] >= z_bot) & (sim.x[:, 2] < z_top))
177 #import ipdb; ipdb.set_trace()
178 if idx[0].size > 0:
179 # xvel[iz,i] = numpy.mean(numpy.abs(sim.vel[I,0]))
180 # xvel[iz, i] = numpy.mean(sim.vel[idx, 0])
181 xvel[iz, i] = numpy.mean(sim.vel[idx, 0])/sim.shearVelocity()
182
183 loaded_contacts = 0
184 if calculateforcechains:
185 if i > 0 and calculateforcechainhistory:
186 loaded_contacts_prev = numpy.copy(loaded_contacts)
187 pairs_prev = numpy.copy(sim.pairs)
188
189 loaded_contacts = sim.findLoadedContacts(
190 threshold=sim.currentNormalStress()*2.)
191 # sim.currentNormalStress()/1000.)
192 n[i] = numpy.size(loaded_contacts)
193
194 sim.findCoordinationNumber()
195 coordinationnumber[i] = sim.findMeanCoordinationNumber()
196
197 if calculateforcechainhistory:
198 nfound = 0
199 if i > 0:
200 for a in loaded_contacts[0]:
201 for b in loaded_contacts_prev[0]:
202 if (sim.pairs[:, a] == pairs_prev[:, b]).all():
203 nfound += 1
204
205 nkept[i] = nfound
206 print nfound
207
208 print coordinationnumber[i]
209
210
211 if calculateforcechains:
212 numpy.savetxt(sid + '-fc.txt', (n, nkept, coordinationnumber))
213 else:
214 if plotforcechains:
215 n, nkept, coordinationnumber = numpy.loadtxt(sid + '-fc.txt')
216
217 # Transform time from model time to real time [s]
218 #t = t/t_DEM_to_t_real
219
220 # integrate velocities to displacement along x (xdispint)
221 # Taylor two term expansion
222 xdispint = numpy.zeros_like(t)
223 v_limit = 2.78e-3 # 1 m/hour (WIP)
224 dt = (t[1] - t[0])
225 dt2 = dt*2.
226 for i in numpy.arange(t.size):
227 if i > 0 and i < t.size-1:
228 acc = (numpy.min([v[i+1], v_limit]) - numpy.min([v[i-1], v_limit]))/dt2
229 xdispint[i] = xdispint[i-1] +\
230 numpy.min([v[i], v_limit])*dt + 0.5*acc*dt**2
231 elif i == t.size-1:
232 xdispint[i] = xdispint[i-1] + numpy.min([v[i], v_limit])*dt
233
234
235 ############
236 # PLOTTING #
237 ############
238 bbox_x = 0.03
239 bbox_y = 0.96
240 verticalalignment = 'top'
241 horizontalalignment = 'left'
242 fontweight = 'bold'
243 bbox = {'facecolor': 'white', 'alpha': 1.0, 'pad': 3}
244
245 # Time in days
246 #t = t/(60.*60.*24.)
247
248 nplots = 4
249 fig = plt.figure(figsize=[3.5, 8.0/4.0*nplots])
250
251 # ax1: v, ax2: I
252 ax1 = plt.subplot(nplots, 1, 1)
253 lns0 = ax1.plot(t, v, '-k', label="$v$",
254 linewidth=linewidth)
255 # lns1 = ax1.plot(t, sigma_eff/1000., '-k', label="$\\sigma'$")
256 # lns2 = ax1.plot(t, tau_def/1000., '-r', label="$\\tau$")
257 # ns2 = ax1.plot(t, tau_def/1000., '-r')
258 #lns3 = ax1.plot(t, tau_eff/1000., '-r', label="$\\tau'$", linewidth=linewidth)
259
260 ax1.set_ylabel('Shear velocity $v$ [ms$^{-1}$]')
261
262 ax2 = ax1.twinx()
263 ax2color = 'blue'
264 # lns4 = ax2.plot(t, p_f_top/1000.0 + 80.0, '-',
265 # color=ax2color,
266 # label='$p_\\text{f}^\\text{forcing}$')
267 # lns5 = ax2.semilogy(t, I, ':',
268 lns5 = ax2.semilogy(t, I, '--',
269 color=ax2color,
270 label='$I$', linewidth=linewidth)
271 ax2.set_ylabel('Inertia number $I$ [-]')
272 ax2.yaxis.label.set_color(ax2color)
273 for tl in ax2.get_yticklabels():
274 tl.set_color(ax2color)
275 #ax2.legend(loc='upper right')
276 #lns = lns0+lns1+lns2+lns3+lns4+lns5
277 #lns = lns0+lns1+lns2+lns3+lns5
278 #lns = lns1+lns3+lns5
279 #lns = lns0+lns3+lns5
280 lns = lns0+lns5
281 labs = [l.get_label() for l in lns]
282 # ax2.legend(lns, labs, loc='upper right', ncol=3,
283 # fancybox=True, framealpha=legend_alpha)
284 #ax1.set_ylim([-30, 200])
285 #ax2.set_ylim([-115, 125])
286
287 # ax1.text(bbox_x, bbox_y, 'A',
288 ax1.text(bbox_x, bbox_y, 'a',
289 horizontalalignment=horizontalalignment,
290 verticalalignment=verticalalignment,
291 fontweight=fontweight, bbox=bbox,
292 transform=ax1.transAxes)
293
294
295 # ax3: mu, ax4: unused
296 ax3 = plt.subplot(nplots, 1, 2, sharex=ax1)
297 #ax3.semilogy(t, mu, 'k', linewidth=linewidth)
298 alpha=1.0
299 if smooth_friction:
300 alpha=0.5
301 ax3.plot(t, mu, 'k', alpha=alpha, linewidth=linewidth)
302 if smooth_friction:
303 # smoothed
304 ax3.plot(t, smooth(mu, smooth_window), linewidth=2)
305 # label='', linewidth=1,
306 # alpha=alpha, color=color[c])
307 ax3.set_ylabel('Bulk friction $\\mu = \\tau\'/N\'$ [-]')
308 #ax3.set_ylabel('Bulk friction $\\mu = \\tau\'/N$ [-]')
309
310 # ax3.text(bbox_x, bbox_y, 'B',
311 ax3.text(bbox_x, bbox_y, 'b',
312 horizontalalignment=horizontalalignment,
313 verticalalignment=verticalalignment,
314 fontweight=fontweight, bbox=bbox,
315 transform=ax3.transAxes)
316
317
318 # ax7: n, ax8: unused
319 ax7 = plt.subplot(nplots, 1, 3, sharex=ax1)
320 if plotforcechains:
321 ax7.plot(t[:n.size], coordinationnumber, 'k', linewidth=linewidth)
322 ax7.set_ylabel('Coordination number $\\bar{n}$ [-]')
323 #ax7.semilogy(t, n - nkept, 'b', label='$\Delta n_\\text{heavy}$')
324 #ax7.set_ylim([1.0e1, 2.0e4])
325 #ax7.set_ylim([-0.2, 9.8])
326 ax7.set_ylim([-0.2, 5.2])
327
328 # ax7.text(bbox_x, bbox_y, 'D',
329 ax7.text(bbox_x, bbox_y, 'c',
330 horizontalalignment=horizontalalignment,
331 verticalalignment=verticalalignment,
332 fontweight=fontweight, bbox=bbox,
333 transform=ax7.transAxes)
334
335
336 # ax9: porosity or xvel, ax10: unused
337 #ax9 = plt.subplot(nplots, 1, 5, sharex=ax1)
338 ax9 = plt.subplot(nplots, 1, 4, sharex=ax1)
339 #poros_min = 0.375
340 #poros_max = 0.45
341 poros[:, 0] = poros[:, 2] # remove erroneous porosity increase
342 #cmap = matplotlib.cm.get_cmap('Blues_r')
343 cmap = matplotlib.cm.get_cmap('afmhot')
344 # im9 = ax9.pcolormesh(t, zpos_c, poros,
345 #zpos_c = zpos_c[:-1]
346
347 xvel = xvel[:-1]
348 # xvel[xvel < 0.0] = 0.0 # ignore negative velocities
349 # im9 = ax9.pcolormesh(t, zpos_c, poros,
350 im9 = ax9.pcolormesh(t, zpos_c, xvel,
351 cmap=cmap,
352 #vmin=poros_min, vmax=poros_max,
353 #norm=matplotlib.colors.LogNorm(vmin=1.0e-8, vmax=xvel.max()),
354 shading='goraud',
355 rasterized=rasterized)
356 ax9.set_ylim([zpos_c[0], sim.w_x[0]])
357 ax9.set_ylabel('Vertical position $z$ [m]')
358
359 cbaxes = fig.add_axes([0.32, 0.1, 0.4, 0.01]) # x,y,w,h
360
361 # ax9.add_patch(matplotlib.patches.Rectangle(
362 # (3.0, 0.04), # x,y
363 # 15., # dx
364 # .15, # dy
365 # fill=True,
366 # linewidth=1,
367 # facecolor='white'))
368 # ax9.add_patch(matplotlib.patches.Rectangle(
369 # (1.5, 0.04), # x,y
370 # 7., # dx
371 # .15, # dy
372 # fill=True,
373 # linewidth=1,
374 # facecolor='white',
375 # alpha=legend_alpha))
376
377 cb9 = plt.colorbar(im9, cax=cbaxes,
378 #ticks=[poros_min, poros_min + 0.5*(poros_max-poros_min), poros_max],
379 #ticks=[xvel.min(), xvel.min() + 0.5*(xvel.max()-xvel.min()), xvel.max()],
380 orientation='horizontal',
381 # extend='min',
382 cmap=cmap)
383 # cmap.set_under([8./255., 48./255., 107./255.]) # for poros
384 # cmap.set_under([1.0e-3, 1.0e-3, 1.0e-3]) # for xvel
385 # cb9.outline.set_color('w')
386 cb9.outline.set_edgecolor('w')
387
388 from matplotlib import ticker
389 tick_locator = ticker.MaxNLocator(nbins=4)
390 cb9.locator = tick_locator
391 cb9.update_ticks()
392
393 #cb9.set_label('Mean horizontal porosity [-]')
394 cb9.set_label('Norm. avg. horiz. vel. [-]', color='w')
395 '''
396 ax9.text(0.5, 0.4, 'Mean horizontal porosity [-]\\\\',
397 horizontalalignment='center',
398 verticalalignment='center',
399 bbox={'facecolor':'white', 'alpha':1.0, 'pad':3})
400 '''
401 cb9.solids.set_rasterized(rasterized)
402
403 # change text color of colorbar to white
404 #axes_obj = plt.getp(im9, 'axes')
405 #plt.setp(plt.getp(axes_obj, 'yticklabels'), color='w')
406 #plt.setp(plt.getp(axes_obj, 'xticklabels'), color='w')
407 #plt.setp(plt.getp(cb9.ax.axes, 'yticklabels'), color='w')
408 cb9.ax.yaxis.set_tick_params(color='w')
409 # cb9.yaxis.label.set_color(ax2color)
410 for tl in cb9.ax.get_xticklabels():
411 tl.set_color('w')
412 cb9.ax.yaxis.set_tick_params(color='w')
413
414 # ax9.text(bbox_x, bbox_y, 'E',
415 ax9.text(bbox_x, bbox_y, 'd',
416 horizontalalignment=horizontalalignment,
417 verticalalignment=verticalalignment,
418 fontweight=fontweight, bbox=bbox,
419 transform=ax9.transAxes)
420
421
422 plt.setp(ax1.get_xticklabels(), visible=False)
423 #plt.setp(ax2.get_xticklabels(), visible=False)
424 plt.setp(ax3.get_xticklabels(), visible=False)
425 #plt.setp(ax4.get_xticklabels(), visible=False)
426 #plt.setp(ax5.get_xticklabels(), visible=False)
427 #plt.setp(ax6.get_xticklabels(), visible=False)
428 plt.setp(ax7.get_xticklabels(), visible=False)
429 #plt.setp(ax8.get_xticklabels(), visible=False)
430
431 ax1.set_xlim([numpy.min(t), numpy.max(t)])
432 #ax2.set_ylim([1e-5, 1e-3])
433 ax3.set_ylim([-0.2, 1.2])
434
435 ax9.set_xlabel('Time [s]')
436 fig.tight_layout()
437 plt.subplots_adjust(hspace=0.05)
438
439 filename = sid + '-combined.' + outformat
440 plt.savefig(filename)
441 plt.close()
442 #shutil.copyfile(filename, '/home/adc/articles/own/3/graphics/' + filename)
443 print(filename)
444
445 # Loop through parameter values
446 for effective_stress in effective_stresses:
447 for velfac in velfacs:
448 for mu_s in mu_s_vals:
449 for mu_d in mu_s_vals:
450
451 jobname = jobname_prefix + '{}Pa-v={}-mu_s={}-mu_d={}'.format(
452 effective_stress,
453 velfac,
454 mu_s,
455 mu_d)
456
457 print(jobname)
458 sim = sphere.sim(jobname)
459 #sim.visualize('shear')
460 rateStatePlot(jobname)