thalfshear-darcy-creep-dynamics.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-creep-dynamics.py (13651B)
---
1 #!/usr/bin/env python
2 import matplotlib
3 matplotlib.use('Agg')
4 import shutil
5
6 import os
7 import numpy
8 import sphere
9 import matplotlib.pyplot as plt
10 import matplotlib.patches
11 import matplotlib.colors
12 import matplotlib.ticker
13 import matplotlib.cm
14
15 matplotlib.rcParams.update({'font.size': 7, 'font.family': 'sans-serif'})
16 matplotlib.rc('text', usetex=True)
17 matplotlib.rcParams['text.latex.preamble'] = [r"\usepackage{amsmath}"]
18 matplotlib.rc('grid', linestyle=':', linewidth=0.2)
19
20 # import seaborn as sns
21 # sns.set(style='ticks', palette='Set2')
22 # sns.set(style='ticks', palette='colorblind')
23 # sns.set(style='ticks', palette='muted')
24 # sns.set(style='ticks', palette='pastel')
25 # sns.set(style='ticks', palette='pastel')
26 # sns.set(style='ticks')
27 # sns.despine() # remove right and top spines
28
29 outformat = 'pdf'
30
31 scatter = False
32 plotContacts = False
33 # plotContacts = False
34 plotForceChains = True
35 # plotForceChains = False
36 plotHistograms = False
37 plotCombinedHistogram = True
38
39 sid = 'halfshear-darcy-sigma0=80000.0' +\
40 '-k_c=3.5e-13-mu=1.04e-07-ss=10000.0-A=70000.0-f=0.2'
41 timescaling = 5.787e-5
42
43 n = 6 # wave number
44 period = 500 # 1 period = 500 steps
45 steps = numpy.arange(n*period + 0.25*period, n*period + 0.75*period,
46 dtype=numpy.int)
47 ''' # 6.5, 6.6, 6.7 d
48 plotsteps = numpy.array([
49 n*period + 0.50*period,
50 n*period + 0.60*period,
51 n*period + 0.65*period],
52 dtype=numpy.int) # slow creep, fast creep, slip
53 '''
54 plotsteps = numpy.array([
55 n*period + 0.50*period,
56 n*period + 0.60*period,
57 n*period + 0.65*period],
58 dtype=numpy.int) # slow creep, fast creep, slip
59 # steps[-1] - 10], dtype=numpy.int) # slow creep, fast creep, slip
60 # plotsteps = numpy.array([1670])
61 contactfigs = []
62 contactidx = []
63
64 datalists = [[], [], []]
65 strikelists = [[], [], []]
66 diplists = [[], [], []]
67 forcemagnitudes = [[], [], []]
68 alphas = [[], [], []]
69 f_n_maxs = [[], [], []]
70 taus = [[], [], []]
71 ts = [[], [], []]
72 Ns = [[], [], []]
73 vs = [[], [], []]
74
75 # f_min = 1.0
76 # f_max = 1.0e16
77 # lower_limit = 0.3
78 lower_limit = 0.08
79 # upper_limit = 0.5
80 upper_limit = 1.0
81 f_n_max = 500 # for force chain plots
82
83 N = numpy.zeros_like(steps, dtype=numpy.float64)
84 t = numpy.zeros_like(steps, dtype=numpy.float64)
85 v = numpy.zeros_like(steps, dtype=numpy.float64)
86
87 vel_avg = numpy.zeros_like(N)
88 angvel_avg = numpy.zeros_like(N)
89
90 # insert plot positions
91 Lx = [.17, .37, .65]
92 Ly = [.3, .7, .4]
93 dx = .23
94 dy = .23
95 xytext_x = [Lx[0]+.15*dx, Lx[1]+.6*dx, Lx[2]+1.15*dx]
96 xytext_y = [Ly[0]+dy+.07, Ly[1]+.15, .2]
97
98 fc_x = [Lx[0]-.007, Lx[1]+dx+.08, Lx[2]-.007]
99 fc_y = [Ly[0]-.07*dy, Ly[1]+.035, Ly[2]-.7*dy]
100 fc_dx = [dx, dx, dx]
101 fc_dy = [dy*.7, dy*.7, dy*.7]
102
103
104 sim = sphere.sim(sid, fluid=True)
105 # print(sim.status())
106 t_DEM_to_t_real = timescaling
107
108 i = 0
109 i_scatter = 0
110 for step in steps:
111
112 sim.readstep(step, verbose=False)
113 if i == 0:
114 L = sim.L
115
116 N[i] = sim.currentNormalStress('defined')
117 t[i] = sim.currentTime()
118 vel_avg[i] = numpy.average(sim.vel[:, 0])
119 angvel_avg[i] = numpy.average(sim.angvel[:, 1])
120 v[i] = sim.shearVelocity()
121
122 if plotContacts:
123 outfolder = '../img_out/' + sim.id() + '/'
124 if not os.path.exists(outfolder):
125 os.makedirs(outfolder)
126 sim.plotContacts(outfolder=outfolder,
127 alpha=0.9,
128 lower_limit=lower_limit,
129 upper_limit=upper_limit)
130
131 if (step == plotsteps).any():
132 datalists[i_scatter], strikelists[i_scatter], diplists[i_scatter],\
133 forcemagnitudes[i_scatter], alphas[i_scatter], \
134 f_n_maxs[i_scatter] = sim.plotContacts(return_data=True,
135 alpha=0.9,
136 lower_limit=lower_limit,
137 upper_limit=upper_limit)
138 # contactfigs.append(
139 # sim.plotContacts(return_fig=True,
140 # f_min=f_min,
141 # f_max=f_max))
142 # datalists.append(data)
143 # strikelists.append(strikelist)
144 # diplists.append(diplist)
145 # forcemagnitudes.append(forcemagnitude)
146 # alphas.append(alpha)
147 # f_n_maxs.append(f_n_max)
148 ts[i_scatter] = t[i]
149 Ns[i_scatter] = N[i]
150 vs[i_scatter] = sim.shearVelocity()
151 # taus.append(sim.shearStress('defined'))
152 taus[i_scatter] = sim.shearStress('defined')
153
154 # contactidx.append(step)
155 i_scatter += 1
156
157 i += 1
158
159 # PLOTTING ######################################################
160
161 # Time in days
162 scalingfactor = 1./t_DEM_to_t_real / (24.*3600.)
163 t_scaled = t*scalingfactor
164
165 fig = plt.figure(figsize=[3.5, 3.5])
166
167 #plt.figtext(0.05, 0.95, 'A', horizontalalignment='left', weight='bold')
168 #plt.figtext(0.05, 0.35, 'B', horizontalalignment='left', weight='bold')
169 plt.figtext(0.05, 0.95, 'a', horizontalalignment='left', weight='bold')
170 plt.figtext(0.05, 0.35, 'b', horizontalalignment='left', weight='bold')
171
172 # ax1 = plt.subplot(1, 1, 1)
173 ax1 = plt.subplot2grid((2, 3), (0, 0), colspan=3)
174
175 ax1.add_patch(matplotlib.patches.Rectangle((6.61, -100.), .2, 500., alpha=0.7,
176 facecolor='orange',
177 edgecolor='none'))
178 # shade stick periods
179 # collection = matplotlib.collections.BrokenBarHCollection.span_where(
180 # t_scaled, ymin=-20.0, ymax=200.0,
181 # where=numpy.isclose(v, 0.0),
182 # facecolor='black', alpha=0.2,
183 # linewidth=0)
184 # ax1.add_collection(collection)
185
186 lns0 = ax1.plot(t_scaled, N/1000., '-k', label='$N$', clip_on=False)
187 # ax1.plot([t_scaled[int(len(t_scaled)*0.60)], 10],
188 # [taus[0]/.3/1000., taus[0]/.3/1000.], '--', color='gray')
189 ax1.set_xlabel('Time $t$ [d]')
190 ax1.set_ylabel('Effective normal stress $N$ [kPa]')
191
192 ax2 = ax1.twinx()
193 lns1 = ax2.semilogy(t_scaled, numpy.abs(vel_avg)*timescaling, '-b',
194 label='$|\\bar{\\boldsymbol{v}}_x|$',
195 clip_on=False)
196 lns2 = ax2.semilogy(t_scaled, numpy.abs(angvel_avg)*timescaling, '-r',
197 label='$|\\bar{\\boldsymbol{\\omega}}_y|$',
198 clip_on=False)
199 ax2.set_ylabel('Linear $\\bar{\\boldsymbol{v}}_x$ [m/s] or '
200 + 'angular velocity $\\bar{\\boldsymbol{\\omega}}_y$ [rad/s]')
201
202 ax1.set_ylim([-20, 150])
203 ax1.text(6.42, -5., 'Creep', horizontalalignment='center')
204 ax1.text(6.68, -5., 'Slip', horizontalalignment='center')
205
206 ax1.spines['right'].set_visible(False)
207 ax1.spines['top'].set_visible(False)
208 # Only show ticks on the left and bottom spines
209 ax1.yaxis.set_ticks_position('left')
210 ax1.xaxis.set_ticks_position('bottom')
211
212 lns = lns0+lns1+lns2
213 labs = [l.get_label() for l in lns]
214 # bbox_to_anchor=(0.5, 1.12) for legend centered above box
215 ax2.legend(lns, labs, loc='upper center', ncol=3, bbox_to_anchor=(0.5, 1.25),
216 fancybox=True, framealpha=1.0)
217
218 ax1.set_xlim([numpy.min(t_scaled), numpy.max(t_scaled)])
219 # ax1.locator_params(axis='x', nbins=5)
220 ax1.locator_params(axis='y', nbins=4)
221
222 # Contact scatter plots
223
224 # Scatter plot 1
225 for sc in range(len(Lx)):
226 xy = (ts[sc]*scalingfactor, Ns[sc]/1000.)
227 # print xytext
228 # print xy
229 '''
230 ax1.annotate('',
231 xytext=(xytext_x[sc], xytext_y[sc]),
232 textcoords='axes fraction',
233 xy=xy, xycoords='data',
234 arrowprops=dict(arrowstyle='->', connectionstyle='arc3'))
235 '''
236
237 if scatter:
238 axsc1 = fig.add_axes([Lx[sc], Ly[sc], dx, dy], polar=True)
239 cs = axsc1.scatter(strikelists[sc], 90. - diplists[sc], marker='o',
240 c=forcemagnitudes[sc],
241 s=forcemagnitudes[sc]/f_n_maxs[sc]*5.,
242 alpha=alphas[sc],
243 edgecolors='none',
244 vmin=f_n_maxs[sc]*lower_limit,
245 vmax=f_n_maxs[sc]*upper_limit,
246 cmap=matplotlib.cm.get_cmap('afmhot_r'))
247 # norm=matplotlib.colors.LogNorm())
248 # tick locations
249 # thetaticks = numpy.arange(0,360,90)
250 # set ticklabels location at 1.3 times the axes' radius
251 # ax.set_thetagrids(thetaticks, frac=1.3)
252 axsc1.set_xticklabels([])
253 axsc1.set_yticklabels([])
254
255 if sc == 0:
256 title = '1. Slow creep'
257 elif sc == 1:
258 title = '2. Fast creep'
259 elif sc == 2:
260 title = '3. Slip'
261 else:
262 title = 'Unknown'
263 axsc1.set_title('\\textbf{' + title + '}', fontsize=7)
264 if upper_limit < 1.0:
265 cbar = plt.colorbar(cs, extend='max', fraction=0.035, pad=0.04)
266 else:
267 cbar = plt.colorbar(cs, fraction=0.045, pad=0.04)
268 # cbar.set_label('$||\\boldsymbol{f}_n||$')
269 cbar.set_label('$\\boldsymbol{f}_\\text{n}$ [N]')
270 cbar.locator = matplotlib.ticker.MaxNLocator(nbins=4)
271 cbar.update_ticks()
272
273 # plot defined max compressive stress from tau/N ratio
274 axsc1.scatter(0., numpy.degrees(numpy.arctan(taus[sc]/Ns[sc])),
275 marker='+', c='none', edgecolor='red', s=100)
276 axsc1.scatter(0., numpy.degrees(numpy.arctan(taus[sc]/Ns[sc])),
277 marker='o', c='none', edgecolor='red', s=100)
278 '''
279 ax.sc1scatter(0., # actual stress
280 numpy.degrees(numpy.arctan(
281 self.shearStress('effective')/
282 self.currentNormalStress('effective'))),
283 marker='+', color='blue', s=300)
284 '''
285
286 axsc1.set_rmax(90)
287 axsc1.set_rticks([])
288 # axsc1.grid(False)
289
290 if plotHistograms:
291 axhistload1 = fig.add_axes([Lx[sc], Ly[sc], dx, dy*.5])
292 axhistload1.hist(datalists[sc][:, 6], alpha=0.75, facecolor='gray',
293 log=True)
294 # plt.yscale('log', nonposy='clip')
295 axhistload1.set_xlim([0, 50.])
296 axhistload1.set_xlabel('Contact load [N]')
297
298 axdipload1 = fig.add_axes([Lx[sc], Ly[sc]-.8*dy, dx, dy*.5])
299 axdipload1.hist(diplists[sc], alpha=0.75, facecolor='gray',
300 log=False)
301 axdipload1.set_xlabel('Contact angle')
302 # axdipload1.set_ylabel('Count')
303
304 # force chain plot
305
306 if plotForceChains:
307 axfc1 = plt.subplot2grid((2, 3), (1, sc))
308 # axfc1 = fig.add_axes([fc_x[sc], fc_y[sc], fc_dx[sc], fc_dy[sc]])
309
310 data = datalists[sc]
311
312 # find the max. value of the normal force
313 # f_n_max = numpy.amax(data[:,6])
314
315 # specify the lower limit of force chains to do statistics on
316 f_n_lim = lower_limit * f_n_max
317
318 # find the indexes of these contacts
319 I = numpy.nonzero(data[:, 6] >= f_n_lim)
320
321 # color = matplotlib.cm.spectral(data[:,6]/f_n_max)
322 for i in I[0]:
323
324 x1 = data[i, 0]
325 # y1 = data[i, 1]
326 z1 = data[i, 2]
327 x2 = data[i, 3]
328 # y2 = data[i, 4]
329 z2 = data[i, 5]
330 f_n = data[i, 6]
331
332 lw_max = 1.2
333 if f_n >= f_n_max:
334 lw = lw_max
335 else:
336 lw = (f_n - f_n_lim)/(f_n_max - f_n_lim)*lw_max
337
338 # print lw
339 axfc1.plot([x1, x2], [z1, z2], '-k', linewidth=lw)
340 # axfc1.plot([x1,x2], [z1,z2], '-', linewidth=lw, color=color)
341
342 if sc == 0:
343 title = 'Slow creep'
344 elif sc == 1:
345 title = 'Fast creep'
346 elif sc == 2:
347 title = 'Slip'
348 else:
349 title = 'Unknown'
350 axfc1.set_title(title, fontsize=7)
351 if sc == 0:
352 axfc1.set_ylabel('Stress network', fontsize=7)
353 axfc1.set_xlabel('$t$ = {:.3} d\n'.format(ts[sc]*scalingfactor) +
354 '$v$ = {:.3} m/s'.format(vs[sc]*scalingfactor))
355
356 axfc1.spines['right'].set_visible(False)
357 axfc1.spines['left'].set_visible(False)
358 # Only show ticks on the left and bottom spines
359 axfc1.xaxis.set_ticks_position('none')
360 axfc1.yaxis.set_ticks_position('none')
361 axfc1.set_xticklabels([])
362 axfc1.set_yticklabels([])
363 axfc1.set_xlim([numpy.min(data[I[0], 0]), numpy.max(data[I[0], 0])])
364 axfc1.set_ylim([numpy.min(data[I[0], 2]), numpy.max(data[I[0], 2])])
365 axfc1.set_aspect('equal')
366
367 fig.tight_layout()
368 # plt.subplots_adjust(hspace=0.05)
369 # plt.subplots_adjust(right=0.82)
370
371 filename = sid + '-creep-dynamics.' + outformat
372 plt.savefig(filename)
373 plt.close()
374 shutil.copyfile(filename, '/home/adc/articles/own/3/graphics/' + filename)
375 print(filename)
376
377 if plotCombinedHistogram:
378 fig = plt.figure(figsize=[3.5, 3.5])
379
380 ax1 = plt.subplot(1, 1, 1)
381
382 f_min = 0.0
383 f_max = 0.0
384 for sc in range(len(Lx)):
385 if f_max < numpy.max(datalists[sc][:, 6]):
386 f_max = numpy.max(datalists[sc][:, 6])
387 f_max = numpy.ceil(f_max/10.)*10.
388 print(f_max)
389
390 for sc in range(len(Lx)):
391 # axhistload1.hist(datalists[sc][:,6], alpha=0.75, facecolor='gray',
392 # log=True)
393 hist, bin_edges = numpy.histogram(datalists[sc][:, 6],
394 range=(f_min, f_max))
395
396 ax1.semilogy((bin_edges[1:] - bin_edges[:-1])/2 + bin_edges[:-1],
397 hist, 'o-', label='$N$ = {:3.1f} kPa'.format(Ns[sc]/1000.))
398
399 # plt.yscale('log', nonposy='clip')
400 ax1.set_xlabel('Contact load [N]')
401 ax1.set_ylabel('Number of contacts')
402 ax1.legend(loc='upper right', fancybox=True, framealpha=1.0)
403
404 filename = sid + '-creep-dynamics-hist.' + outformat
405 plt.savefig(filename)
406 plt.close()
407 shutil.copyfile(filename, '/home/adc/articles/own/3/graphics/' + filename)
408 print(filename)