tpermeability-results.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
---
tpermeability-results.py (5010B)
---
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
15 #dp_list = numpy.array([1.0e3, 2.0e3, 4.0e3, 10.0e3, 20.0e3, 40.0e3])
16 dp_list = numpy.array([1.0e3, 2.0e3, 4.0e3, 10.0e3])
17 cvals = [1.0, 0.1, 0.01]
18 c_phi = 1.0
19
20 K = [[], [], []]
21 dpdz = [[], [], []]
22 Q = [[], [], []]
23 phi_bar = [[], [], []]
24 Re = [[], [], []]
25 fp_fsum = [[], [], []]
26
27
28 c = 0
29 for c_grad_p in cvals:
30
31 sids = []
32 for dp in dp_list:
33 if c_grad_p == 1.0:
34 sids.append('permeability-dp=' + str(dp))
35 else:
36 sids.append('permeability-dp=' + str(dp) + '-c_phi=' + \
37 str(c_phi) + '-c_grad_p=' + str(c_grad_p))
38
39 K[c] = numpy.zeros(len(sids))
40 dpdz[c] = numpy.zeros_like(K[c])
41 Q[c] = numpy.zeros_like(K[c])
42 phi_bar[c] = numpy.zeros_like(K[c])
43 Re[c] = numpy.zeros_like(K[c])
44 fp_fsum[c] = numpy.zeros_like(K[c])
45 i = 0
46
47 for sid in sids:
48 if os.path.isfile('../output/' + sid + '.status.dat'):
49 pc = PermeabilityCalc(sid, plot_evolution=False,
50 print_results=False, verbose=False)
51 K[c][i] = pc.conductivity()
52 pc.findPressureGradient()
53 pc.findCrossSectionalFlux()
54 dpdz[c][i] = pc.dPdL[2]
55 Q[c][i] = pc.Q[2]
56 pc.findMeanPorosity()
57 #pc.plotEvolution()
58 phi_bar[c][i] = pc.phi_bar
59
60 sim = sphere.sim(sid, fluid=True)
61 sim.readlast(verbose=False)
62 Re[c][i] = numpy.mean(sim.ReynoldsNumber())
63
64 #sim.writeVTKall()
65
66 # find magnitude of fluid pressure force and total interaction force
67 '''
68 fp_magn = numpy.empty(sim.np)
69 fsum_magn = numpy.empty(sim.np)
70 for i in numpy.arange(sim.np):
71 fp_magn[i] = sim.f_p[i,:].dot(sim.f_p[i,:])
72 fsum_magn[i] = sim.f_sum[i,:].dot(sim.f_sum[i,:])
73
74 fp_fsum[c][i] = numpy.mean(fp_magn/fsum_magn)
75 # interaction forces not written in these old output files!
76 '''
77
78
79 else:
80 print(sid + ' not found')
81
82 i += 1
83
84 # produce VTK files
85 #for sid in sids:
86 #sim = sphere.sim(sid, fluid=True)
87 #sim.writeVTKall()
88 c += 1
89
90 fig = plt.figure(figsize=(8,12))
91
92 ax1 = plt.subplot(3,1,1)
93 ax2 = plt.subplot(3,1,2, sharex=ax1)
94 ax3 = plt.subplot(3,1,3, sharex=ax1)
95 #ax4 = plt.subplot(4,1,4, sharex=ax1)
96 colors = ['g', 'r', 'c', 'y']
97 #lines = ['-', '--', '-.', ':']
98 lines = ['-', '-', '-', '-']
99 #markers = ['o', 'x', '^', '+']
100 markers = ['x', 'x', 'x', 'x']
101 for c in range(len(cvals)):
102 dpdz[c] /= 1000.0
103 #plt.plot(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
104 #plt.semilogx(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
105 #plt.semilogy(dpdz[c], K[c], 'o-', label='$c$ = %.2f' % (cvals[c]))
106 ax1.loglog(dpdz[c], K[c], label='$c$ = %.2f' % (cvals[c]),
107 linestyle=lines[c], marker=markers[c], color=colors[c])
108 ax2.semilogx(dpdz[c], phi_bar[c], label='$c$ = %.2f' % (cvals[c]),
109 linestyle=lines[c], marker=markers[c], color=colors[c])
110 ax3.loglog(dpdz[c], Re[c], label='$c$ = %.2f' % (cvals[c]),
111 linestyle=lines[c], marker=markers[c], color=colors[c])
112 #ax4.loglog(dpdz[c], fp_fsum[c], label='$c$ = %.2f' % (cvals[c]),
113 #linestyle=lines[c], marker=markers[c], color='black')
114
115 ax1.set_ylabel('Hydraulic conductivity $K$ [ms$^{-1}$]')
116 #ax1.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
117
118 ax2.set_ylabel('Mean porosity $\\bar{\\phi}$ [-]')
119
120 ax3.set_ylabel('Mean Reynolds number $\\bar{Re}$ [-]')
121
122 #ax4.set_ylabel('$\\bar{\\boldsymbol{f}_{\\Delta p}/\\bar{\\boldsymbol{f}_\\text{pf}}$ [-]')
123
124 ax3.set_xlabel('Pressure gradient $\\Delta p/\\Delta z$ [kPa m$^{-1}$]')
125
126 plt.setp(ax1.get_xticklabels(), visible=False)
127 plt.setp(ax2.get_xticklabels(), visible=False)
128
129 ax1.grid()
130 ax2.grid()
131 ax3.grid()
132 #ax4.grid()
133
134 #plt.subplot(3,1,2)
135 #plt.xlabel('Pressure gradient $\\Delta p/\\Delta z$ [Pa m$^{-1}$]')
136 #plt.ylabel('Hydraulic flux $Q$ [m$^3$s$^{-1}$]')
137 #plt.plot(dpdz, Q, '+')
138 #plt.grid()
139
140 #plt.subplot(3,1,3)
141 #plt.xlabel('Pressure gradient $\\Delta p/\\Delta z$ [Pa m$^{-1}$]')
142 #plt.ylabel('Mean porosity $\\bar{\\phi}$ [-]')
143 #plt.plot(dpdz, phi_bar, '+')
144 #plt.grid()
145
146 ax1.legend(loc='best', prop={'size':18}, fancybox=True, framealpha=0.5)
147 ax2.legend(loc='best', prop={'size':18}, fancybox=True, framealpha=0.5)
148 ax3.legend(loc='best', prop={'size':18}, fancybox=True, framealpha=0.5)
149
150 plt.tight_layout()
151 plt.subplots_adjust(hspace = .12)
152 filename = 'permeability-dpdz-vs-K-vs-c.pdf'
153 #print(os.getcwd() + '/' + filename)
154 plt.savefig(filename)
155 shutil.copyfile(filename, '/home/adc/articles/own/2-org/' + filename)
156 print(filename)