thalfshear-darcy-fft.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-fft.py (4245B)
---
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.fftpack
16
17 #sigma0_list = numpy.array([1.0e3, 2.0e3, 4.0e3, 10.0e3, 20.0e3, 40.0e3])
18 #sigma0 = 10.0e3
19 sigma0 = 20000.0
20 k_c_vals = [3.5e-13, 3.5e-15]
21 mu_f = 1.797e-06
22 velfac = 1.0
23 #cvals = [1.0, 0.1, 0.01]
24 #cvals = [1.0]
25
26
27 shear_strain = [[], [], [], []]
28 friction = [[], [], [], []]
29 dilation = [[], [], [], []]
30 p_min = [[], [], [], []]
31 p_mean = [[], [], [], []]
32 p_max = [[], [], [], []]
33 f_n_mean = [[], [], [], []]
34 f_n_max = [[], [], [], []]
35 v_f_z_mean = [[], [], [], []]
36 t_total = []
37
38 fluid=True
39
40 # dry shear
41 #sid = 'shear-sigma0=' + sys.argv[1] + '-hw'
42 # halfshear-darcy-sigma0=20000.0-k_c=3.5e-13-mu=1.797e-06-velfac=1.0-shear
43 sid = 'halfshear-sigma0=' + str(sigma0) + '-shear'
44 sim = sphere.sim(sid)
45 sim.readlast(verbose=False)
46 sim.visualize('shear')
47 shear_strain[0] = sim.shear_strain
48 #shear_strain[0] = numpy.arange(sim.status()+1)
49 friction[0] = sim.tau/sim.sigma_eff
50 dilation[0] = sim.dilation
51 t_total.append(sim.time_total[0])
52
53
54 # wet shear
55 c = 1
56 for c in numpy.arange(1,len(k_c_vals)+1):
57 k_c = k_c_vals[c-1]
58
59 # halfshear-darcy-sigma0=20000.0-k_c=3.5e-13-mu=1.797e-06-velfac=1.0-shear
60 sid = 'halfshear-darcy-sigma0=' + str(sigma0) + '-k_c=' + str(k_c) + \
61 '-mu=' + str(mu_f) + '-velfac=' + str(velfac) + '-shear'
62 #sid = 'halfshear-sigma0=' + str(sigma0) + '-c_v=' + str(c_v) +\
63 #'-c_a=0.0-velfac=1.0-shear'
64 if os.path.isfile('../output/' + sid + '.status.dat'):
65
66 sim = sphere.sim(sid, fluid=fluid)
67 shear_strain[c] = numpy.zeros(sim.status())
68 friction[c] = numpy.zeros_like(shear_strain[c])
69 dilation[c] = numpy.zeros_like(shear_strain[c])
70
71 sim.readlast(verbose=False)
72 sim.visualize('shear')
73 t_total.append(sim.time_total[0])
74 shear_strain[c] = sim.shear_strain
75 #shear_strain[c] = numpy.arange(sim.status()+1)
76 friction[c] = sim.tau/sim.sigma_eff
77 dilation[c] = sim.dilation
78
79 else:
80 print(sid + ' not found')
81
82 # produce VTK files
83 #for sid in sids:
84 #sim = sphere.sim(sid, fluid=True)
85 #sim.writeVTKall()
86 c += 1
87
88 fig = plt.figure(figsize=(8,8)) # (w,h)
89 #fig.subplots_adjust(hspace=0.0)
90
91 #ax1 = plt.subplot(111)
92 ax1 = plt.subplot(211)
93 ax2 = plt.subplot(212, sharex=ax1)
94 alpha = 1.0
95 #ax1.plot(shear_strain[0], friction[0], label='dry', linewidth=1, alpha=alpha)
96
97 color = ['b','g','r','c']
98 for c in numpy.arange(0,len(k_c_vals)+1):
99
100 if c == 0:
101 label = 'dry'
102 elif c == 1:
103 label = 'wet, relatively permeable'
104 elif c == 2:
105 label = 'wet, relatively impermeable'
106 else:
107 label = '$k_c$ = %.1e m$^2$' % (k_c_vals[c-1])
108
109 str_arr = shear_strain[c][200:1999]
110 dil_arr = dilation[c][200:1999]
111 t = numpy.linspace(0.0, sim.time_total[0], shear_strain[c].size)
112
113 freqs = scipy.fftpack.fftfreq(str_arr.size, t[1]-t[0])
114 str_yf = numpy.abs(scipy.fftpack.fft(str_arr))
115 dil_yf = numpy.abs(scipy.fftpack.fft(dil_arr))
116
117 ax1.plot(freqs, str_yf, label=label, linewidth=1, alpha=alpha)
118 ax2.plot(freqs, dil_yf, label=label, linewidth=1, alpha=alpha)
119
120 ax2.set_xlabel('Frequency [s$^{-1}$]')
121
122 #ax1.set_ylabel('Shear friction $\\tau/\\sigma\'$ [-]')
123 #ax2.set_ylabel('Dilation $\\Delta h/(2r)$ [-]')
124
125 ax1.set_xlim([0.0,12.0])
126 plt.setp(ax1.get_xticklabels(), visible=False)
127
128 ax1.grid()
129 ax2.grid()
130
131 legend_alpha=0.5
132 ax1.legend(loc='best', prop={'size':18}, fancybox=True,
133 framealpha=legend_alpha)
134 #ax2.legend(loc='lower right', prop={'size':18}, fancybox=True,
135 #framealpha=legend_alpha)
136
137 #ax1.set_ylim([-0.1, 1.9])
138
139 plt.tight_layout()
140 plt.subplots_adjust(hspace=0.05)
141 #filename = 'shear-' + str(int(sigma0/1000.0)) + 'kPa-stress-dilation.pdf'
142 filename = 'halfshear-darcy-fft.pdf'
143 #print(os.getcwd() + '/' + filename)
144 plt.savefig(filename)
145 shutil.copyfile(filename, '/home/adc/articles/own/2/graphics/' + filename)
146 print(filename)