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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3509 by gross, Thu Apr 28 05:06:24 2011 UTC revision 3510 by gross, Fri May 13 06:09:49 2011 UTC
# Line 26  from esys.weipa import saveVTK Line 26  from esys.weipa import saveVTK
26  from esys.escript import unitsSI as U  from esys.escript import unitsSI as U
27  from coalgas import *  from coalgas import *
28  import time  import time
29  from esys.finley import ReadMesh, Rectangle  from esys.finley import ReadMesh, Rectangle, Brick
30  from esys.escript.pdetools import Locator  from esys.escript.pdetools import Locator
31    CONST_G = 9.81  * U.m/U.sec**2
32    P_0=1.*U.atm
33    
 CELL_X=800*U.m  
 CELL_Y=800*U.m  
 CELL_Z=10*U.m  
34    
35    CELL_X=2640*U.ft
36    CELL_Y=2640*U.ft
37    CELL_Z=33*U.ft
38    
39  L_Z=10*U.m  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    
 L_X=168*U.km/3  
 L_Y=168*U.km/3  
47  N_X=int(L_X/CELL_X)  N_X=int(L_X/CELL_X)
48  N_Y=int(L_Y/CELL_Y)  N_Y=int(L_Y/CELL_Y)
49    N_Z=int(L_Z/CELL_Z)
50    
51    
52  OUTPUT_DIR="results"  OUTPUT_DIR="results"
53    
54    
55  PERM_F_X = 10 * U.mDarcy  PERM_F_X = 100 * U.mDarcy
56  PERM_F_Y = 10 * U.mDarcy  PERM_F_Y = 100 * U.mDarcy
57  PERM_F_Z = 1e-4 * U.mDarcy  PERM_F_Z = 1e-4 * U.mDarcy
58    
59  EQUIL = {  EQUIL = {
# Line 56  EQUIL = { Line 64  EQUIL = {
64      "GWC_PCOW" : 0. * U.psi      "GWC_PCOW" : 0. * U.psi
65  }  }
66    
67    TOPS = 2310 * U.ft
68  PHI_F_0=0.01  PHI_F_0=0.01
69  SIGMA = 1. /U.ft**2  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  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  #[0.1 * U.day ] *20
76    
# Line 94  LANGMUIR = [ Line 104  LANGMUIR = [
104  [ 1000  * U.psi , 0.01045748 * U.Mscf/U.ft**3],  [ 1000  * U.psi , 0.01045748 * U.Mscf/U.ft**3],
105  [ 1200  * U.psi , 0.01126912 * U.Mscf/U.ft**3] ]  [ 1200  * U.psi , 0.01126912 * U.Mscf/U.ft**3] ]
106    
107    
108  PVDG = [  PVDG = [
109  [ 14.70 * U.psi ,200.3800 * U.Barrel/U.Mscf , 0.012025 * U.cPoise ] , # psi, rb/Mscf,    [ 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 ] ,  [ 20.00 * U.psi ,146.0600 * U.Barrel/U.Mscf , 0.012030 * U.cPoise ] ,
# Line 179  SWFN = [ Line 190  SWFN = [
190    
191    
192  wellspecs = {  wellspecs = {
193    'P1' : { "X0" :[ (N_X/2+0.5)*CELL_X,  (N_Y/2+0.5)*CELL_Y],    '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,             "r"  : 0.8333 * U.ft,
195             "s"  : 0,             "s"  : 0,
196             "Q"  : 2000*U.Barrel/U.day,             "Q"  : [0., 2000*U.Barrel/U.day, 2000*U.Barrel/U.day, 0.],
197             "BHP" : 75*U.psi,             "BHP" : 75*U.psi,
198             "schedule" : [0.*U.yr, 2*U.yr]             "schedule" : [0.*U.yr,  1.*U.day, 2.*U.yr-1*U.day, 2*U.yr]
199           }           }
200  }  }
201    
# Line 193  wellspecs = { Line 204  wellspecs = {
204  # print input  # print input
205  print("<%s> Execution started."%time.asctime())  print("<%s> Execution started."%time.asctime())
206  domain=Rectangle(N_X, N_Y, l0=L_X, l1=L_Y)  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)  #domain=ReadMesh("reservoir.fly",2)
210  print("<%s> Domain has been generated."%time.asctime())  print("<%s> Domain has been generated."%time.asctime())
211    
# Line 203  print "cell size in y direction = %f m"% Line 216  print "cell size in y direction = %f m"%
216  print "fracture permeability in x direction= %f mD"%(PERM_F_X/(U.mDarcy))  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))  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))  print "fracture permeability in z direction= %f mD"%(PERM_F_Z/(U.mDarcy))
 print "initial porosity in fractured rock= %f"%EQUIL["DATUM_PRESS"]  
219    
220    
221  mkDir(OUTPUT_DIR)  mkDir(OUTPUT_DIR)
# Line 214  print("<%s> Mesh set up completed."%time Line 226  print("<%s> Mesh set up completed."%time
226  well_P1=VerticalPeacemanWell('P1', domain, BHP_limit=wellspecs['P1' ]["BHP"],  well_P1=VerticalPeacemanWell('P1', domain, BHP_limit=wellspecs['P1' ]["BHP"],
227                                  Q=wellspecs['P1']["Q"],                                  Q=wellspecs['P1']["Q"],
228                                  r=wellspecs['P1']["r"],                                  r=wellspecs['P1']["r"],
229                                  X0=[ wellspecs['P1' ]["X0"][0], wellspecs['P1']["X0"][1]] ,                                  X0=[ wellspecs['P1' ]["X0"][0], wellspecs['P1']["X0"][1], wellspecs['P1']["X0"][2]] ,
230                                  D=[CELL_X, CELL_Y, L_Z],                                  D=[CELL_X, CELL_Y, CELL_Z],
231                                  perm=[PERM_F_X, PERM_F_Y, PERM_F_Z],                                  perm=[PERM_F_X, PERM_F_Y, PERM_F_Z],
232                                  schedule=wellspecs['P1']["schedule"],                                  schedule=wellspecs['P1']["schedule"],
233                                  s=wellspecs['P1']["s"], decay_factor = 1)                                  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,  model = PorosityOneHalfModel(domain,
239                               phi_f=Porosity(phi_0=PHI_F_0, p_0=EQUIL["DATUM_PRESS"], p_ref=ROCK["p_ref"], C = ROCK["C"]),                               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 ] ),                               L_g=InterpolationTable([ l[0] for l in LANGMUIR ], [ l[1] for l in LANGMUIR ] ),
241                   perm_f_0=PERM_F_X,                   perm_f_0=PERM_F_X,
242                   perm_f_1=PERM_F_Y,                   perm_f_1=PERM_F_Y,
# Line 230  model = PorosityOneHalfModel(domain, Line 245  model = PorosityOneHalfModel(domain,
245                   k_g= InterpolationTable([ l[0] for l in SGFN ], [ l[1] for l in SGFN ], obeyBounds=False ),                     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"]),                         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 ] ),                   mu_g = InterpolationTable([ l[0] for l in PVDG ], [ l[2] for l in PVDG ] ),
248                   rho_w = WaterDensity(B_ref=PVTW["B_ref"], p_ref = PVTW["p_ref"], C=PVTW["C"], gravity=GRAVITY["water"]),                   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"]),                   rho_g=GasDensity( p = [ l[0] for l in PVDG ], B = [ l[1] for l in PVDG ], gravity=GRAVITY["gas"]),
250                   sigma=SIGMA,                   sigma=SIGMA,
251                   A_mg=DIFFCOAL["D"],                   A_mg=DIFFCOAL["D"],
252                   f_rg=DIFFCOAL["f_r"],                   f_rg=DIFFCOAL["f_r"],
253                   wells=[ well_P1, ] )                   wells=[ well_P1, ], g= CONST_G)
254        # this needs to be revised:.
255  model.setInitialState(p=EQUIL["DATUM_PRESS"], S_fg=0,  C_mg=None)  model.setInitialState(S_fg=0,  c_mg=None, p_top=p_top, p_bottom=p_bottom)
256  model.getPDEOptions().setVerbosityOn()  model.getPDEOptions().setVerbosityOn()
257  model.getPDEOptions().setSolverMethod(model.getPDEOptions().DIRECT)  model.getPDEOptions().setSolverMethod(model.getPDEOptions().DIRECT)
258  model.setIterationControl(iter_max=1, rtol=1.e-4, verbose=True, xi=1.0)  model.setIterationControl(iter_max=22, rtol=1.e-4, verbose=True)
259  print "<%s> Problem set up completed."%time.asctime()  print "<%s> Problem set up completed."%time.asctime()
260    
261  p, S_fg, C_mg=model.getState()  p, S_fg, c_mg=model.getState()
262    
263  FN=os.path.join(OUTPUT_DIR, "state.%d.vtu"%0)  FN=os.path.join(OUTPUT_DIR, "state.%d.vtu"%0)
264  saveVTK(FN,p=p, S_fg=S_fg, C_mg=C_mg)  saveVTK(FN,p=p, S_fg=S_fg, c_mg=c_mg)
265  print "<%s> Initial state saved to file %s."%(time.asctime(),FN)  print "<%s> Initial state saved to file %s."%(time.asctime(),FN)
266    
267    
# Line 258  for dt in DT: Line 273  for dt in DT:
273        
274    model.update(dt)    model.update(dt)
275    
276    p, S_fg, C_mg=model.getState()    p, S_fg, c_mg=model.getState()
277        
278    FN=os.path.join(OUTPUT_DIR, "state.%d.vtu"%(n_t+1))    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)    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    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 )    print "<%s>State %s saved to file %s."%(time.asctime(),n_t+1,FN )
282    
283    n_t+=1    n_t+=1
284    t+=dt    t+=dt
   
 #DDD 1.0 620.718083467 0.0 2000.0  
 #DDD 4.0 616.356950541 0.0 2000.0  
 #DDD 13.0 602.731301066 0.0 2000.0  
 #DDD 30.5 559.432880289 0.0210138086777 2000.0  
 #DDD 61.0 476.375599867 0.242467255289 2000.0  
 #DDD 91.5 292.277380173 2.30397035606 2000.0  
 #DDD 122.0 75.0 7942.07091041 1748.12930007  
285      if n_t > 2:
286        1/0

Legend:
Removed from v.3509  
changed lines
  Added in v.3510

  ViewVC Help
Powered by ViewVC 1.1.26