thalfshear-darcy-rate-strength.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-rate-strength.py (7497B)
---
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 import scipy.optimize
16
17 import seaborn as sns
18 #sns.set(style='ticks', palette='Set2')
19 sns.set(style='ticks', palette='colorblind')
20 #sns.set(style='ticks', palette='muted')
21 #sns.set(style='ticks', palette='pastel')
22 sns.despine() # remove right and top spines
23
24 pressures = True
25 zflow = False
26 contact_forces = False
27 smooth_friction = True
28
29
30 sids = [
31 'halfshear-darcy-sigma0=20000.0-k_c=3.5e-13-mu=1.797e-06-velfac=1.0-shear',
32 'halfshear-darcy-sigma0=20000.0-k_c=3.5e-14-mu=1.797e-06-velfac=1.0-shear',
33 'halfshear-darcy-sigma0=20000.0-k_c=3.5e-15-mu=1.797e-06-velfac=1.0-shear',
34
35 #'halfshear-darcy-sigma0=20000.0-k_c=3.5e-15-mu=3.594e-07-velfac=1.0-shear',
36 #'halfshear-darcy-sigma0=20000.0-k_c=3.5e-15-mu=6.11e-07-velfac=1.0-shear',
37
38 'halfshear-darcy-sigma0=20000.0-k_c=3.5e-15-mu=1.797e-07-velfac=1.0-shear',
39 'halfshear-darcy-sigma0=20000.0-k_c=3.5e-15-mu=1.797e-08-velfac=1.0-shear',
40
41
42 #'halfshear-sigma0=20000.0-shear'
43
44 # from halfshear-darcy-perm.sh
45 #'halfshear-darcy-sigma0=20000.0-k_c=3.5e-13-mu=1.797e-07-velfac=1.0-shear',
46 #'halfshear-darcy-sigma0=20000.0-k_c=3.5e-13-mu=1.797e-08-velfac=1.0-shear',
47 #'halfshear-darcy-sigma0=20000.0-k_c=3.5e-13-mu=1.797e-09-velfac=1.0-shear',
48 ]
49 fluids = [
50 True,
51 True,
52 True,
53 True,
54 True,
55 #False,
56 True,
57 True,
58 True,
59 ]
60
61 def smooth(x, window_len=10, window='hanning'):
62 """smooth the data using a window with requested size.
63
64 This method is based on the convolution of a scaled window with the signal.
65 The signal is prepared by introducing reflected copies of the signal
66 (with the window size) in both ends so that transient parts are minimized
67 in the begining and end part of the output signal.
68
69 input:
70 x: the input signal
71 window_len: the dimension of the smoothing window
72 window: the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'
73 flat window will produce a moving average smoothing.
74
75 output:
76 the smoothed signal
77
78 example:
79
80 import numpy as np
81 t = np.linspace(-2,2,0.1)
82 x = np.sin(t)+np.random.randn(len(t))*0.1
83 y = smooth(x)
84
85 see also:
86
87 numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve
88 scipy.signal.lfilter
89
90 TODO: the window parameter could be the window itself if an array instead of a string
91 """
92
93 if x.ndim != 1:
94 raise ValueError, "smooth only accepts 1 dimension arrays."
95
96 if x.size < window_len:
97 raise ValueError, "Input vector needs to be bigger than window size."
98
99 if window_len < 3:
100 return x
101
102 if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
103 raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"
104
105 s=numpy.r_[2*x[0]-x[window_len:1:-1], x, 2*x[-1]-x[-1:-window_len:-1]]
106 #print(len(s))
107
108 if window == 'flat': #moving average
109 w = numpy.ones(window_len,'d')
110 else:
111 w = getattr(numpy, window)(window_len)
112 y = numpy.convolve(w/w.sum(), s, mode='same')
113 return y[window_len-1:-window_len+1]
114
115
116 max_step = 400
117 friction = numpy.zeros(max_step)
118
119 velratios = []
120 peakfrictions = []
121
122 sim = sphere.sim(sids[0], fluid=True)
123 sim.readfirst(verbose=False)
124 fig = plt.figure(figsize=(3.5, 3))
125
126 if False:
127 it=0
128 for sid in sids:
129
130 print '\n' + sid
131 sim.id(sid)
132 sim.fluid=fluids[it]
133 it += 1
134
135 sim.readfirst(verbose=False)
136
137 velratio = 1.0
138 if sim.fluid:
139 if sim.k_c[0] > 3.5e-15 and sim.mu[0] < 1.797e-06:
140 print 'Large permeability and low viscosity'
141 velratio = 3.5e-15/sim.k_c[0] \
142 * sim.mu[0]/1.797e-06
143 elif sim.k_c[0] > 3.5e-15:
144 print 'Large permeability'
145 velratio = 3.5e-15/sim.k_c[0]
146 elif sim.mu < 1.797e-06:
147 print 'Low viscosity'
148 velratio = sim.mu[0]/1.797e-06
149 elif numpy.isclose(sim.k_c[0], 3.5e-15) and \
150 numpy.isclose(sim.mu[0], 1.797e-06):
151 print 'Normal permeability and viscosity'
152 velratio = 1.
153 else:
154 raise Exception('Could not determine velratio')
155 else:
156 velratio = 0.001
157
158 print 'velratio = ' + str(velratio)
159
160 for i in numpy.arange(max_step):
161 sim.readstep(i+1, verbose=False)
162
163 friction[i] = sim.shearStress('effective')/\
164 sim.currentNormalStress('defined')
165
166 smoothfriction = smooth(friction, 20, window='hanning')
167 peakfriction = numpy.amax(smoothfriction[14:])
168
169 plt.plot(numpy.arange(max_step), friction, alpha=0.3, color='b')
170 plt.plot(numpy.arange(max_step), smoothfriction, color='b')
171 plt.title(str(peakfriction))
172 plt.savefig(sid + '-friction.pdf')
173 plt.clf()
174
175 velratios.append(velratio)
176 peakfrictions.append(peakfriction)
177 else:
178 velratios =\
179 [0.01,
180 0.099999999999999992,
181 1.0,
182 0.10000000000000001,
183 #0.010000000000000002,
184 #0.001,
185 #0.0001,
186 #1.0000000000000001e-05
187 ]
188 peakfrictions = \
189 [0.619354807290315,
190 0.6536161052814875,
191 0.70810354077280957,
192 0.64649301787774571,
193 #0.65265739261434697,
194 #0.85878138368962764,
195 #1.6263903846405066,
196 #0.94451353171977692
197 ]
198
199
200 #def frictionfunction(velratio, a, b=numpy.min(peakfrictions)):
201 def frictionfunction(x, a, b):
202 #return numpy.exp(a*velratio) + numpy.min(peakfrictions)
203 #return -numpy.log(a*velratio) + numpy.min(peakfrictions)
204 #return a*numpy.log(velratio) + b
205 return a*10**(x) + b
206
207 popt, pvoc = scipy.optimize.curve_fit(
208 frictionfunction,
209 peakfrictions,
210 velratios)
211 print popt
212 print pvoc
213 a=popt[0]
214 b=popt[1]
215
216 #a = 0.025
217 #b = numpy.min(peakfrictions)
218 #b = numpy.max(peakfrictions)
219
220 shearvel = sim.shearVel()/1000.*numpy.asarray(velratios) * (3600.*24.*365.25)
221
222 '''
223 plt.semilogx(
224 #[1, numpy.min(shearvel)],
225 [1, 6],
226 #[1, 1000],
227 #[numpy.min(peakfrictions), numpy.min(peakfrictions)],
228 [0.615, 0.615],
229 '--', color='gray', linewidth=2)
230 '''
231
232 #plt.semilogx(velratios, peakfrictions, 'o')
233 #plt.semilogx(shearvel, peakfrictions, 'o')
234 plt.semilogx(shearvel, peakfrictions, 'o')
235 #plt.plot(velratios, peakfrictions, 'o')
236
237 xfit = numpy.linspace(numpy.min(velratios), numpy.max(velratios))
238 #yfit = frictionfunction(xfit, popt[0], popt[1])
239 yfit = frictionfunction(xfit, a, b)
240 #plt.semilogx(xfit, yfit)
241 #plt.plot(xfit, yfit)
242
243 plt.xlabel('Shear velocity [m a$^{-1}$]')
244 plt.ylabel('Peak shear friction, $\\max(\\tau/\\sigma_0)$ [-]')
245 #plt.xlim([0.8, 110])
246 #plt.xlim([0.008, 1.1])
247 plt.xlim([1., 1000,])
248 plt.ylim([0.6, 0.72])
249 plt.tight_layout()
250 sns.despine() # remove right and top spines
251 filename = 'halfshear-darcy-rate-strength.pdf'
252 plt.savefig(filename)
253 plt.close()
254 print(filename)