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

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26