/[escript]/trunk/finley/test/python/runcoalgas.py
ViewVC logotype

Contents of /trunk/finley/test/python/runcoalgas.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (20 months, 1 week ago) by jfenwick
File MIME type: text/x-python
File size: 13476 byte(s)
Make everyone sad by touching all the files

Copyright dates update

1 #######################################################
2 #
3 # Copyright (c) 2003-2018 by The University of Queensland
4 # http://www.uq.edu.au
5 #
6 # Primary Business: Queensland, Australia
7 # Licensed under the Apache License, version 2.0
8 # http://www.apache.org/licenses/LICENSE-2.0
9 #
10 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 # Development 2012-2013 by School of Earth Sciences
12 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 #
14 ##############################################################################
15 """
16 Coal Seam gasL ECLIPSE test case
17 """
18
19 from __future__ import print_function, division
20
21 __copyright__="""Copyright (c) 2003-2018 by The University of Queensland
22 http://www.uq.edu.au
23 Primary Business: Queensland, Australia"""
24 __license__="""Licensed under the Apache License, version 2.0
25 http://www.apache.org/licenses/LICENSE-2.0"""
26 __url__="https://launchpad.net/escript-finley"
27
28
29 from esys.escript import *
30 from esys.weipa import saveVTK
31 from esys.escript import unitsSI as U
32 from coalgas import *
33 import time
34 from esys.finley import ReadMesh, Rectangle, Brick
35 from esys.escript.pdetools import Locator
36
37 SAVE_VTK=True and False
38 CONST_G = 9.81 * U.m/U.sec**2
39 P_0=1.*U.atm
40
41 CELL_X=2640*U.ft
42 CELL_Y=2640*U.ft
43 CELL_Z=33*U.ft
44
45 TOP=2310*U.ft
46
47 L_Z=CELL_Z
48
49 L_X=CELL_X*210
50 L_Y=CELL_Y*210
51 L_X=CELL_X*70
52 L_Y=CELL_Y*70
53
54 N_X=int(L_X/CELL_X)
55 N_Y=int(L_Y/CELL_Y)
56 N_Z=int(L_Z/CELL_Z)
57
58 OUTPUT_DIR="results"
59
60 PERM_F_X = 100 * U.mDarcy
61 PERM_F_Y = 100 * U.mDarcy
62 PERM_F_Z = 1e-4 * U.mDarcy
63
64 EQUIL = {
65 "DATUM_DEPTH" : 2310. * U.ft,
66 "DATUM_PRESS" : 1000. * U.psi,
67 "GWC_DEPTH" : -1000. * U.ft,
68 "GWC_PCOW" : 0. * U.psi
69 }
70
71 TOPS = 2310 * U.ft
72 PHI_F_0=0.01
73 SIGMA = 1. /U.ft**2
74
75 #DT=[0.1* U.day]*9+[1 * U.day,3* U.day,9* U.day, 17.5*U.day] + [ 30.5*U.day ] *20
76 DT=[1 * U.day,3* U.day,9* U.day, 17.5*U.day] + [ 30.5*U.day ] *20
77 DT=[.5 * U.day, .5 * U.day, 0.5*3* U.day, 0.5*3* U.day, 0.5*9* U.day, 0.5*9* U.day, 17.5 *.5 *U.day, 17.5*0.5*U.day] + [ 15.25*U.day ] *40
78 DT=[.25 * U.day, .25 * U.day, .25 * U.day, .25 * U.day,
79 0.25*3* U.day, 0.25*3* U.day, 0.25*3* U.day, 0.25*3* U.day,
80 0.25*9* U.day, 0.25*9* U.day, 0.25*9* U.day, 0.25*9* U.day,
81 17.5 *.25 *U.day, 17.5*0.25*U.day, 17.5 *.25 *U.day, 17.5*0.25*U.day] + \
82 [ 15.25*U.day * 0.5] *80
83 DT=[.125 * U.day, .125 * U.day, .125 * U.day, .125 * U.day, .125 * U.day, .125 * U.day, .125 * U.day, .125 * U.day,
84 0.125*3* U.day, 0.125*3* U.day, 0.125*3* U.day, 0.125*3* U.day, 0.125*3* U.day, 0.125*3* U.day, 0.125*3* U.day, 0.125*3* U.day,
85 0.125*9* U.day, 0.125*9* U.day, 0.125*9* U.day, 0.125*9* U.day, 0.125*9* U.day, 0.125*9* U.day, 0.125*9* U.day, 0.125*9* U.day,
86 17.5 *.125 *U.day, 17.5*0.125*U.day, 17.5 *.125 *U.day, 17.5*0.125*U.day, 17.5 *.125 *U.day, 17.5*0.125*U.day, 17.5 *.125 *U.day, 17.5*0.125*U.day] + \
87 [ 15.25*U.day * 0.25] *160
88 DT=[.125 * U.day, .125 * U.day, .125 * U.day, .125 * U.day, .125 * U.day, .125 * U.day, .125 * U.day, .125 * U.day]
89 # DT=[1 * U.day,2* U.day] + [4*U.day ] *200
90 DT=[.25 * U.day, .25 * U.day, .25 * U.day, .25 * U.day]
91 #DT=[.125 * U.day, .125 * U.day, .125 * U.day, .125 * U.day, .125 * U.day, .125 * U.day, .125 * U.day, .125 * U.day]
92 DT=[1./10. * U.day]*10 + [3./10. * U.day]*10 + [ 9./10.* U.day ] *10 + [ 17.5/10 *U.day] * 10 + [ 30.5/10*U.day ] *200
93
94 #[0.1 * U.day ] *20
95
96 PVTW={ "p_ref" : 1000 * U.psi,
97 "B_ref" : 0.997,
98 "C" : 3.084E-06/U.psi,
99 "mu_ref" : 0.68673 * U.cPoise,
100 "C_v" : 0/U.psi
101 }
102
103 GRAVITY = { "water" : 1.0,
104 "gas" : .553 }
105
106 ROCK = { "p_ref" : 1000 * U.psi,
107 "C" : 3.3E-4 * 1./U.psi }
108
109 DIFFCOAL = { "D" : 0.005 * U.ft**2/U.day,
110 "f_r": 1.}
111
112 LANGMUIR = [
113 [ 0 * U.psi , 0.00000000 * U.Mscf/U.ft**3],
114 [ 100 * U.psi , 0.00213886 * U.Mscf/U.ft**3],
115 [ 200 * U.psi , 0.00383259 * U.Mscf/U.ft**3],
116 [ 300 * U.psi , 0.00520706 * U.Mscf/U.ft**3],
117 [ 400 * U.psi , 0.00634474 * U.Mscf/U.ft**3],
118 [ 500 * U.psi , 0.00730199 * U.Mscf/U.ft**3],
119 [ 600 * U.psi , 0.00811857 * U.Mscf/U.ft**3],
120 [ 700 * U.psi , 0.00882336 * U.Mscf/U.ft**3],
121 [ 800 * U.psi , 0.00943786 * U.Mscf/U.ft**3],
122 [ 900 * U.psi , 0.00997836 * U.Mscf/U.ft**3],
123 [ 1000 * U.psi , 0.01045748 * U.Mscf/U.ft**3],
124 [ 1200 * U.psi , 0.01126912 * U.Mscf/U.ft**3] ]
125
126
127 PVDG = [
128 [ 14.70 * U.psi ,200.3800 * U.Barrel/U.Mscf , 0.012025 * U.cPoise ] , # psi, rb/Mscf,
129 [ 20.00 * U.psi ,146.0600 * U.Barrel/U.Mscf , 0.012030 * U.cPoise ] ,
130 [ 25.00 * U.psi ,116.1461 * U.Barrel/U.Mscf , 0.012034 * U.cPoise ] ,
131 [ 30.00 * U.psi ,96.3132 * U.Barrel/U.Mscf , 0.012038 * U.cPoise ] ,
132 [ 35.00 * U.psi ,82.2113 * U.Barrel/U.Mscf , 0.012043 * U.cPoise ] ,
133 [ 49.33 * U.psi ,57.7891 * U.Barrel/U.Mscf , 0.012055 * U.cPoise ] ,
134 [ 59.00 * U.psi ,48.0866 * U.Barrel/U.Mscf , 0.012064 * U.cPoise ] ,
135 [ 69.00 * U.psi ,40.9441 * U.Barrel/U.Mscf , 0.012073 * U.cPoise ] ,
136 [ 75.00 * U.psi ,37.5839 * U.Barrel/U.Mscf , 0.012078 * U.cPoise ] ,
137 [ 83.00 * U.psi ,33.8685 * U.Barrel/U.Mscf , 0.012085 * U.cPoise ] ,
138 [ 90.00 * U.psi ,31.1661 * U.Barrel/U.Mscf , 0.012092 * U.cPoise ] ,
139 [ 95.00 * U.psi ,29.4827 * U.Barrel/U.Mscf , 0.012097 * U.cPoise ] ,
140 [ 100.00 * U.psi ,27.9698 * U.Barrel/U.Mscf , 0.012101 * U.cPoise ] ,
141 [ 105.00 * U.psi ,26.6028 * U.Barrel/U.Mscf , 0.012106 * U.cPoise ] ,
142 [ 118.60 * U.psi ,23.4749 * U.Barrel/U.Mscf , 0.012119 * U.cPoise ] ,
143 [ 120.00 * U.psi ,23.1937 * U.Barrel/U.Mscf , 0.012120 * U.cPoise ] ,
144 [ 140.00 * U.psi ,19.7977 * U.Barrel/U.Mscf , 0.012140 * U.cPoise ] ,
145 [ 153.23 * U.psi ,18.0443 * U.Barrel/U.Mscf , 0.012153 * U.cPoise ] ,
146 [ 160.00 * U.psi ,17.2607 * U.Barrel/U.Mscf , 0.012159 * U.cPoise ] ,
147 [ 170.00 * U.psi ,16.2188 * U.Barrel/U.Mscf , 0.012169 * U.cPoise ] ,
148 [ 187.86 * U.psi ,14.6373 * U.Barrel/U.Mscf , 0.012188 * U.cPoise ] ,
149 [ 222.49 * U.psi ,12.3027 * U.Barrel/U.Mscf , 0.012224 * U.cPoise ] ,
150 [ 257.13 * U.psi ,10.6038 * U.Barrel/U.Mscf , 0.012262 * U.cPoise ] ,
151 [ 291.76 * U.psi ,9.3134 * U.Barrel/U.Mscf , 0.012301 * U.cPoise ] ,
152 [ 326.39 * U.psi ,8.3001 * U.Barrel/U.Mscf , 0.012341 * U.cPoise ] ,
153 [ 361.02 * U.psi ,7.4835 * U.Barrel/U.Mscf , 0.012383 * U.cPoise ] ,
154 [ 395.66 * U.psi ,6.8114 * U.Barrel/U.Mscf , 0.012425 * U.cPoise ] ,
155 [ 430.29 * U.psi ,6.2491 * U.Barrel/U.Mscf , 0.012470 * U.cPoise ] ,
156 [ 464.92 * U.psi ,5.7715 * U.Barrel/U.Mscf , 0.012515 * U.cPoise ] ,
157 [ 499.55 * U.psi ,5.3610 * U.Barrel/U.Mscf , 0.012562 * U.cPoise ] ,
158 [ 534.19 * U.psi ,5.0043* U.Barrel/U.Mscf , 0.012610 * U.cPoise ] ,
159 [ 568.82 * U.psi ,4.6917 * U.Barrel/U.Mscf , 0.012659 * U.cPoise ] ,
160 [ 603.45 * U.psi ,4.4154 * U.Barrel/U.Mscf , 0.012710 * U.cPoise ] ,
161 [ 638.08 * U.psi ,4.1695 * U.Barrel/U.Mscf , 0.012762 * U.cPoise ] ,
162 [ 672.72 * U.psi ,3.9491 * U.Barrel/U.Mscf , 0.012815 * U.cPoise ] ,
163 [ 707.35 * U.psi ,3.7507 * U.Barrel/U.Mscf , 0.012869 * U.cPoise ] ,
164 [ 741.98 * U.psi ,3.5711 * U.Barrel/U.Mscf , 0.012925 * U.cPoise ] ,
165 [ 776.61 * U.psi ,3.4076 * U.Barrel/U.Mscf , 0.012982 * U.cPoise ] ,
166 [ 811.25 * U.psi ,3.2583 * U.Barrel/U.Mscf , 0.013041 * U.cPoise ] ,
167 [ 845.88 * U.psi ,3.1214 * U.Barrel/U.Mscf , 0.013100 * U.cPoise ] ,
168 [ 880.51 * U.psi ,2.9953 * U.Barrel/U.Mscf , 0.013161 * U.cPoise ] ,
169 [ 915.14 * U.psi ,2.8790 * U.Barrel/U.Mscf , 0.013223 * U.cPoise ] ,
170 [ 949.78 * U.psi ,2.7712 * U.Barrel/U.Mscf , 0.013287 * U.cPoise ] ,
171 [ 984.41 * U.psi ,2.6711 * U.Barrel/U.Mscf , 0.013352 * U.cPoise ] ,
172 [ 1019.00 * U.psi ,2.5781 * U.Barrel/U.Mscf , 0.013418 * U.cPoise ] ,
173 [ 1053.70 * U.psi ,2.4909 * U.Barrel/U.Mscf , 0.013486 * U.cPoise ] ,
174 [ 1088.30 * U.psi ,2.4096 * U.Barrel/U.Mscf , 0.013554 * U.cPoise ] ,
175 [ 1122.90 * U.psi ,2.3334 * U.Barrel/U.Mscf , 0.013624 * U.cPoise ] ,
176 [ 1157.60 * U.psi ,2.2616 * U.Barrel/U.Mscf , 0.013696 * U.cPoise ] ,
177 [ 1192.20 * U.psi ,2.1942 * U.Barrel/U.Mscf , 0.013768 * U.cPoise ] ,
178 [ 1226.80 * U.psi ,2.1307 * U.Barrel/U.Mscf , 0.013842 * U.cPoise ] ,
179 [ 1261.50 * U.psi ,2.0705 * U.Barrel/U.Mscf , 0.013917 * U.cPoise ] ,
180 [ 1296.10 * U.psi ,2.0138 * U.Barrel/U.Mscf , 0.013994 * U.cPoise ] ,
181 [ 1330.70 * U.psi ,1.9600 * U.Barrel/U.Mscf , 0.014072 * U.cPoise ] ,
182 [ 1365.40 * U.psi ,1.9089 * U.Barrel/U.Mscf , 0.014151 * U.cPoise ] ]
183
184
185 SGFN = [
186 [ 0 , 0 , 0 * U.psi],
187 [ 0.05 , 0 , 0 * U.psi ],
188 [ 0.1333 , 0.00610 , 0 * U.psi ],
189 [ 0.2167 , 0.02990 , 0 * U.psi ],
190 [ 0.3 , 0.0759 , 0 * U.psi],
191 [ 0.3833 , 0.1471 , 0 * U.psi],
192 [ 0.46667 , 0.2458 , 0 * U.psi],
193 [ 0.55 , 0.3739 , 0 * U.psi],
194 [ 0.6333 , 0.53300 , 0 * U.psi],
195 [ 0.7167 , 0.7246 , 0 * U.psi],
196 [ 0.8 , 0.95 , 0 * U.psi] ]
197 SWFN = [
198 [ 0.20000 , 0.00000, 0 * U.psi],
199 [ 0.28330 , 0.03280, 0 * U.psi],
200 [ 0.36670 , 0.09270, 0 * U.psi],
201 [ 0.45000 , 0.17030, 0 * U.psi],
202 [ 0.53330 , 0.26220, 0 * U.psi],
203 [ 0.61670 , 0.36650, 0 * U.psi],
204 [ 0.70000 , 0.48170, 0 * U.psi],
205 [ 0.78330 , 0.60710, 0 * U.psi],
206 [ 0.86670 , 0.74170, 0 * U.psi],
207 [ 0.95000 , 0.88500, 0 * U.psi],
208 [ 1.00000 , 1.00000, 0 * U.psi] ]
209
210
211 wellspecs = {
212 'P1' : { "X0" :[ (N_X/2+0.5)*CELL_X, (N_Y/2+0.5)*CELL_Y, 0.5*CELL_Z],
213 "r" : 0.8333 * U.ft,
214 "s" : 0,
215 "Q" : [0., 2000*U.Barrel/U.day ],
216 "BHP" : 75*U.psi,
217 "schedule" : [0.*U.yr, 2*U.yr]
218 }
219 }
220
221
222 # print input
223 print(("<%s> Execution started."%time.asctime()))
224 DIM=2
225 domain=Rectangle(N_X, N_Y, l0=L_X, l1=L_Y)
226
227 N=1000
228 for I in wellspecs:
229 domain.setTagMap(I, N)
230 N+=1
231 domain.addDiracPoint(wellspecs[I]["X0"][:DIM], I)
232 print(("<%s> Well %s introduced to domain."%(time.asctime(), I)))
233
234 #domain=Brick(N_X, N_Y,N_Z,l0=L_X, l1=L_Y,l2=L_Z)
235
236 print(("<%s> Domain has been generated."%time.asctime()))
237
238 print("length x-direction = %f km"%(sup(domain.getX()[0])/U.km))
239 print("cell size in x direction = %f m"%(CELL_X/U.m))
240 print("length y-direction = %f km"%(sup(domain.getX()[1])/U.km))
241 print("cell size in y direction = %f m"%(CELL_Y/U.m))
242 print("fracture permeability in x direction= %f mD"%(PERM_F_X/(U.mDarcy)))
243 print("fracture permeability in y direction= %f mD"%(PERM_F_Y/(U.mDarcy)))
244 print("fracture permeability in z direction= %f mD"%(PERM_F_Z/(U.mDarcy)))
245
246
247 mkDir(OUTPUT_DIR)
248
249 print(("<%s> Mesh set up completed."%time.asctime()))
250 well_P1=VerticalPeacemanWell('P1', domain, BHP_limit=wellspecs['P1' ]["BHP"],
251 Q=wellspecs['P1']["Q"],
252 r=wellspecs['P1']["r"],
253 X0=[ wellspecs['P1' ]["X0"][0], wellspecs['P1']["X0"][1], wellspecs['P1']["X0"][2]] ,
254 D=[CELL_X, CELL_Y, CELL_Z],
255 perm=[PERM_F_X, PERM_F_Y, PERM_F_Z],
256 schedule=wellspecs['P1']["schedule"],
257 s=wellspecs['P1']["s"])
258 rho_w = WaterDensity(B_ref=PVTW["B_ref"], p_ref = PVTW["p_ref"], C=PVTW["C"], gravity=GRAVITY["water"])
259 p_top = EQUIL["DATUM_PRESS"] + P_0
260 p_bottom=p_top + CONST_G * CELL_Z * rho_w(p_top)
261
262 model = PorosityOneHalfModel(domain,
263 phi_f=Porosity(phi_0=PHI_F_0, p_0=(p_bottom +p_top)/2., p_ref=ROCK["p_ref"], C = ROCK["C"]),
264 L_g=InterpolationTable([ l[0] for l in LANGMUIR ], [ l[1] for l in LANGMUIR ] ),
265 perm_f_0=PERM_F_X,
266 perm_f_1=PERM_F_Y,
267 perm_f_2=PERM_F_Z,
268 k_w =InterpolationTable([ l[0] for l in SWFN ], [ l[1] for l in SWFN ], obeyBounds=False ),
269 k_g= InterpolationTable([ l[0] for l in SGFN ], [ l[1] for l in SGFN ], obeyBounds=False ),
270 mu_w = WaterViscosity(mu_ref = PVTW["mu_ref"], p_ref=PVTW["p_ref"], C=PVTW["C_v"]),
271 mu_g = InterpolationTable([ l[0] for l in PVDG ], [ l[2] for l in PVDG ] ),
272 rho_w = rho_w,
273 rho_g=GasDensity( p = [ l[0] for l in PVDG ], B = [ l[1] for l in PVDG ], gravity=GRAVITY["gas"]),
274 sigma=SIGMA,
275 A_mg=DIFFCOAL["D"],
276 f_rg=DIFFCOAL["f_r"],
277 wells=[ well_P1, ], g= CONST_G)
278 # this needs to be revised:.
279 model.setInitialState(S_fg=0, c_mg=None, p_top=p_top, p_bottom=p_bottom)
280 model.getPDEOptions().setVerbosityOn()
281 model.getPDEOptions().setSolverMethod(model.getPDEOptions().DIRECT)
282 model.setIterationControl(iter_max=10, rtol=1.e-4, verbose=True)
283 print("<%s> Problem set up completed."%time.asctime())
284 t=0
285 n_t = 0
286
287 p, S_fg, c_mg, BHP, q_gas,q_water =model.getState()
288
289 if SAVE_VTK:
290 FN=os.path.join(OUTPUT_DIR, "state.%d.vtu"%n_t)
291 saveVTK(FN,p=p, S_fg=S_fg, c_mg=c_mg)
292 print("<%s> Initial state saved to file %s."%(time.asctime(),FN))
293 print(t/U.day, well_P1.locator(p)/U.psi, well_P1.locator(S_fg), well_P1.locator(c_mg)/U.Mscf*U.ft**3)
294 print(t/U.day, well_P1.locator(BHP)/U.psi, well_P1.locator(q_gas)/U.Mcf*U.day, well_P1.locator(q_water)/U.Barrel*U.day)
295
296
297 for dt in DT:
298 print("<%s>Time step %d, time = %e days started:"%(time.asctime(), n_t+1, (t+dt)/U.day))
299
300 model.update(dt)
301
302 p, S_fg, c_mg, BHP, q_gas,q_water = model.getState()
303
304 if SAVE_VTK:
305 FN=os.path.join(OUTPUT_DIR, "state.%d.vtu"%(n_t+1))
306 saveVTK(FN,p=p, S_fg=S_fg, c_mg=c_mg)
307 print("<%s>State %s saved to file %s."%(time.asctime(),n_t+1,FN))
308 print((t+dt)/U.day, well_P1.locator(p)/U.psi, well_P1.locator(S_fg), well_P1.locator(c_mg)/U.Mscf*U.ft**3)
309 print((t+dt)/U.day, well_P1.locator(BHP)/U.psi, well_P1.locator(q_gas)/U.Mcf*U.day, well_P1.locator(q_water)/U.Barrel*U.day)
310
311 n_t += 1
312 t += dt
313

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26