thalfshear-darcy-strain-rate.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-strain-rate.py (7158B)
---
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
16 import seaborn as sns
17 #sns.set(style='ticks', palette='Set2')
18 #sns.set(style='ticks', palette='colorblind')
19 sns.set(style='ticks', palette='Set2')
20 sns.despine() # remove chartjunk
21
22 sigma0_list = [20000.0, 80000.0]
23 #cvals = ['dry', 1.0, 0.1, 0.01]
24 #cvals = ['dry', 3.5e-13, 3.5e-15]
25 #cvals = ['dry', 3.5e-13, 3.5e-14, 3.5e-15]
26 #cvals = ['dry', 3.5e-13, 3.5e-14, 3.5e-15]
27 #cvals = ['dry', 1.0]
28 c = 3.5e-15
29 muvals = ['dry', 1.797e-08, 1.797e-07, 1.797e-06]
30 #step = 1999
31
32 for sigma0 in sigma0_list:
33
34 sim = sphere.sim('halfshear-sigma0=' + str(sigma0) + '-shear')
35 sim.readfirst(verbose=False)
36
37
38 # particle z positions
39 zpos_p = [[], [], [], []]
40
41 # cell midpoint cell positions
42 zpos_c = [[], [], [], []]
43
44 # particle x displacements
45 xdisp = [[], [], [], []]
46 xdisp_mean = [[], [], [], []]
47
48 s = 0
49 for mu in muvals:
50
51 if mu == 'dry':
52 fluid = False
53 sid = 'halfshear-sigma0=' + str(sigma0) + '-shear'
54 else:
55 fluid = True
56 sid = 'halfshear-darcy-sigma0=' + str(sigma0) + '-k_c=' + str(c) + \
57 '-mu=' + str(mu) + '-velfac=1.0-shear'
58
59 sim = sphere.sim(sid, fluid=fluid)
60
61 if os.path.isfile('../output/' + sid + '.status.dat'):
62
63 sim.readlast(verbose=False)
64
65 zpos_c[s] = numpy.zeros(sim.num[2]*2)
66 dz = sim.L[2]/(sim.num[2]*2)
67 for i in numpy.arange(sim.num[2]*2):
68 zpos_c[s][i] = i*dz + 0.5*dz
69
70 xdisp[s] = numpy.zeros(sim.np)
71 xdisp_mean[s] = numpy.zeros(sim.num[2]*2)
72
73
74 zpos_p[s][:] = sim.x[:,2]
75
76 xdisp[s][:] = sim.xyzsum[:,0]
77
78 #shear_strain[s] += sim.shearStrain()/nsteps_avg
79
80 # calculate mean values of xdisp and f_pf
81 for iz in numpy.arange(sim.num[2]*2):
82 z_bot = iz*dz
83 z_top = (iz+1)*dz
84 I = numpy.nonzero((zpos_p[s][:] >= z_bot) & (zpos_p[s][:] < z_top))
85 if len(I) > 0:
86 xdisp_mean[s][iz] = numpy.mean(xdisp[s][I])
87
88 # normalize distance
89 max_dist = numpy.nanmax(xdisp_mean[s])
90 xdisp_mean[s] /= max_dist
91
92 else:
93 print(sid + ' not found')
94 s += 1
95
96
97 #fig = plt.figure(figsize=(8,4*(len(steps))+1))
98 #fig = plt.figure(figsize=(8,5*(len(steps))+1))
99 #fig = plt.figure(figsize=(8/2,6/2))
100 fig = plt.figure(figsize=(3.74,3.47)) # 3.14 inch = 80 mm, 3.74 = 95 mm
101 #fig = plt.figure(figsize=(8,6))
102
103 ax = []
104 #linetype = ['-', '--', '-.']
105 #linetype = ['-', '-', '-', '-']
106 linetype = ['-', '--', '-.', ':']
107 #color = ['b','g','c','y']
108 #color = ['k','g','c','y']
109 color = ['y','g','c','k']
110 #color = ['c','m','y','k']
111 for s in numpy.arange(len(muvals)):
112 #for s in numpy.arange(len(cvals)-1, -1, -1):
113
114 ax.append(plt.subplot(111))
115 #ax.append(plt.subplot(len(steps)*100 + 31 + s*3))
116 #ax.append(plt.subplot(len(steps)*100 + 32 + s*3, sharey=ax[s*4+0]))
117 #ax.append(plt.subplot(len(steps)*100 + 33 + s*3, sharey=ax[s*4+0]))
118 #ax.append(ax[s*4+2].twiny())
119
120 if muvals[s] == 'dry':
121 legend = 'dry'
122 elif muvals[s] == 1.797e-06:
123 legend = 'wet, ref.\\ shear vel.'
124 elif muvals[s] == 1.797e-07:
125 legend = 'wet, shear vel.\\ $\\times$ 0.1'
126 elif muvals[s] == 1.797e-08:
127 legend = 'wet, shear vel.\\ $\\times$ 0.01'
128 else:
129 legend = 'wet, $k_c$ = ' + str(muvals[s]) + ' m$^2$'
130
131 #ax[0].plot(xdisp[s], zpos_p[s], ',', color = '#888888')
132 #ax[0].plot(xdisp[s], zpos_p[s], ',', color=color[s], alpha=0.5)
133 ax[0].plot(xdisp_mean[s], zpos_c[s], linetype[s],
134 label=legend)#,
135 #color=color[s],
136 #linewidth=2.0)
137
138 ax[0].set_ylabel('Vertical position $z$ [m]')
139 #ax[0].set_xlabel('$\\boldsymbol{x}^x_\\text{p}$ [m]')
140 ax[0].set_xlabel('Normalized horizontal movement')
141
142 #ax[s*4+0].get_xaxis().set_major_locator(MaxNLocator(nbins=5))
143 #ax[s*4+1].get_xaxis().set_major_locator(MaxNLocator(nbins=5))
144 #ax[s*4+2].get_xaxis().set_major_locator(MaxNLocator(nbins=5))
145
146 #plt.setp(ax[s*4+0].xaxis.get_majorticklabels(), rotation=90)
147 #plt.setp(ax[s*4+1].xaxis.get_majorticklabels(), rotation=90)
148 #plt.setp(ax[s*4+2].xaxis.get_majorticklabels(), rotation=90)
149 #plt.setp(ax[s*4+3].xaxis.get_majorticklabels(), rotation=90)
150
151 #if s == 0:
152 #y = 0.95
153 #if s == 1:
154 #y = 0.55
155
156 #strain_str = 'Shear strain $\\gamma = %.3f$' % (shear_strain[s])
157 #fig.text(0.1, y, strain_str, horizontalalignment='left', fontsize=22)
158 #ax[s*4+0].annotate(strain_str, xytext=(0,1.1), textcoords='figure fraction',
159 #horizontalalignment='left', fontsize=22)
160 #plt.text(0.05, 1.06, strain_str, horizontalalignment='left', fontsize=22,
161 #transform=ax[s*4+0].transAxes)
162 #ax[s*4+0].set_title(strain_str)
163
164 #ax[s*4+0].grid()
165 #ax[s*4+1].grid()
166 #ax[s*4+2].grid()
167 #ax1.legend(loc='lower right', prop={'size':18})
168 #ax2.legend(loc='lower right', prop={'size':18})
169
170 # remove box at top and right
171 ax[0].spines['top'].set_visible(False)
172 ax[0].spines['right'].set_visible(False)
173 # remove ticks at top and right
174 ax[0].get_xaxis().tick_bottom()
175 ax[0].get_yaxis().tick_left()
176 ax[0].get_xaxis().grid(False) # horizontal grid lines
177 #ax[0].get_yaxis().grid(True, linestyle='--', linewidth=0.5) # vertical grid lines
178 ax[0].get_xaxis().grid(True, linestyle=':', linewidth=0.5) # vertical grid lines
179 ax[0].get_yaxis().grid(True, linestyle=':', linewidth=0.5) # vertical grid lines
180
181 # reverse legend order
182 handles, labels = ax[0].get_legend_handles_labels()
183 ax[0].legend(handles[::-1], labels[::-1], loc='best')
184
185 #legend_alpha=0.5
186 #ax[0].legend(loc='lower right', prop={'size':18}, fancybox=True, framealpha=legend_alpha)
187 #ax[0].legend(loc='best', prop={'size':18}, fancybox=True, framealpha=legend_alpha)
188 #ax[0].legend(loc='best')
189 #ax[0].grid()
190 #ax[0].set_xlim([-0.05, 1.01])
191 ax[0].set_xlim([-0.05, 1.05])
192 #ax[0].set_ylim([0.0, 0.47])
193 ax[0].set_ylim([0.20, 0.47])
194 plt.tight_layout()
195 plt.subplots_adjust(wspace = .05)
196 plt.MaxNLocator(nbins=4)
197
198 filename = 'halfshear-darcy-strain-rate.pdf'
199 if sigma0 == 80000.0:
200 filename = 'halfshear-darcy-strain-N80-rate.pdf'
201 plt.savefig(filename)
202 #shutil.copyfile(filename, '/Users/adc/articles/own/2/graphics/' + filename)
203 shutil.copyfile(filename, '/home/adc/articles/own/2/graphics/' + filename)
204 print(filename)