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

Contents of /trunk/finley/test/python/coalgas.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: 27430 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 Gas in Coal Seam (fully coupled version)
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 from esys.escript.linearPDEs import LinearPDE
24 from esys.escript import unitsSI as U
25 from esys.escript.pdetools import Locator
26 from esys.escript import sqrt, log, whereNegative, sup, inf, whereNonPositive, integrate, Function, kronecker, grad, Lsup, clip, maximum, exp, ContinuousFunction
27 from esys.weipa import saveVTK
28 import math
29
30 USE_NODAL_WELL = False or True
31
32 class MaterialProperty(object):
33 """
34 generic class for material properties depending on one (or several) status
35 variable(s)
36 """
37 def __init__(self, *args, **kwargs):
38 """
39 set up
40 """
41 pass
42 def __call__(self, *args, **kwargs):
43 """
44 return value of material property for given status
45 """
46 return self.getValue( *args, **kwargs)
47
48
49 def getValue(self, *args, **kwargs):
50 """
51 return value of material property for given status
52
53 :remark: needs to be overwritten
54 """
55 raise NotImplementedError
56
57 class MaterialPropertyWithDifferential(MaterialProperty):
58 """
59 generic class for material properties depending on one (or several) status
60 variable(s) where the derivative of the property with respect to the
61 status variables is available
62 """
63 def getValueDifferential(self, *args, **kwargs):
64 """
65 returns the value of the derivative of material property for given status variable
66
67 :remark: needs to be overwritten
68 """
69 raise NotImplementedError
70
71 class Porosity(MaterialPropertyWithDifferential):
72 """
73 defines porosity phi as function of pressure
74
75 phi = phi_p_ref /(1 + X + X**2/2 ) with X= C * (p-p_ref)
76
77 phi_p_ref is claculted from the initial porosity phi_0 at pressure p_0
78 """
79 def __init__(self, phi_0, p_0, p_ref=1.*U.atm , C=4.0e-5/U.bar):
80 """
81 """
82 self.phi_p_ref=1 # will be overwritten later
83 self.p_ref=p_ref
84 self.C=C
85 # reset phi_p_ref to get phi(p_0)=phi_0
86 self.phi_p_ref=phi_0/self.getValue(p_0)
87
88 def getValue(self, p):
89 """
90 returns the porosity for given pressure p
91 """
92 X= self.C * ( p - self.p_ref )
93 return self.phi_p_ref * (1. + X * (1. + X/2 ) )
94
95 def getValueDifferential(self, p):
96 """
97 returns the porosity for given pressure p
98 """
99 X= self.C * ( p - self.p_ref )
100 return self.phi_p_ref * self.C * (1. + X )
101
102
103 class WaterDensity(MaterialPropertyWithDifferential):
104 """
105 set water density as
106
107 rho = rho_surf (1 + X + X*X/2) with X= C * ( p - p_ref )
108
109 with rho_surf =rho_s/B_ref * gravity
110 """
111 def __init__(self, B_ref=1., p_ref = 1.*U.atm, gravity=1., C = 0./U.bar, rho_s= 999.014 * U.kg/U.m**3):
112 self.rho_surf = rho_s * gravity
113 self.__rho_0 = self.rho_surf/B_ref
114 self.C=C
115 self.p_ref=p_ref
116
117 def getValue(self, p):
118 """
119 returns the density for given pressure p
120 """
121 X= self.C * ( p - self.p_ref )
122 return self.__rho_0 * (1+ X * (1+ X/2) )
123
124 def getValueDifferential(self, p):
125 """
126 """
127 X= self.C * ( p - self.p_ref )
128 return self.__rho_0 * self.C * (1+ X)
129
130 def getFormationFactor(self, p):
131 return self.rho_surf/self.getValue(p)
132
133
134
135 class WaterViscosity(MaterialProperty):
136 """
137 defines viscosity of water
138
139 mu=mu_ref /(1 + X + X*X/2)
140
141 with X=-C*(p-p_ref)
142 """
143
144 def __init__(self, mu_ref = 0.3*U.cPoise, p_ref = 1.*U.atm , C = 0./U.bar):
145 """
146 set up
147 """
148 self.mu_ref = mu_ref
149 self.p_ref = p_ref
150 self.C=C
151 def getValue(self, p):
152 """
153 returns the viscosity for given pressure p
154 """
155 X= -self.C * ( p - self.p_ref )
156 return self.mu_ref/(1+ X * (1+ X/2) )
157
158 class GasDensity(MaterialPropertyWithDifferential):
159 """
160 set water density as
161
162 rho = gravity * rho_air_s /B(p)
163
164 where B is given by an interpolation table
165 """
166 def __init__(self, p, B, gravity=1., rho_air_s=1.2041*U.kg/U.m**3):
167 self.rho_surf =rho_air_s * gravity
168 self.tab=InterpolationTable(x=p, y=B)
169
170 def getValue(self, p):
171 """
172 returns the density for given pressure p
173 """
174 return self.rho_surf/self.getFormationFactor(p)
175
176 def getValueDifferential(self, p):
177 """
178 """
179 B = self.getFormationFactor(p)
180 dBdp = self.getFormationFactorDifferential(p)
181 return -self.rho_surf * dBdp/(B * B)
182
183 def getFormationFactor(self, p):
184 return self.tab.getValue(p)
185
186 def getFormationFactorDifferential(self, p):
187 return self.tab.getValueDifferential(p)
188
189
190 class InterpolationTable(MaterialPropertyWithDifferential):
191 """
192 a simple 1D interpolation table for escript Data object with non-equidistant nodes
193 """
194 def __init__(self,x=[], y=[], obeyBounds=True):
195 """
196 set up interpolation table. obeyBounds is set an exception is thrown if
197 the interpolation argument is below min(x) or larger than max(x). Otherwise
198 the value for x is set to y[0] for
199 """
200 MaterialPropertyWithDifferential.__init__(self)
201 if len(x) != len(y):
202 raise ValueError,"length of interpolation nodes and value lists need to be identical."
203 if len(x) < 1 :
204 raise ValueError,"length of interpolation nodes a list needs to at least one."
205
206 x_ref=x[0]
207 for i in range(1,len(x)):
208 if x_ref >= x[i]:
209 raise ValueError,"interpolation nodes need to be given in increasing order."
210 x_ref=x[i]
211 self.__x = x
212 self.__y = y
213 self.__obeyBounds = obeyBounds
214
215 def getValue(self, x):
216 """
217 returns the interpolated values of x
218 """
219 X=self.__x
220 Y=self.__y
221
222 x0=X[0]
223 m0=whereNegative(x-x0)
224 if self.__obeyBounds:
225 if sup(m0) > 0:
226 raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
227 out=self.__x[0]
228 for i in range(1,len(X)):
229 z=(Y[i]-Y[i-1])/(X[i]-X[i-1]) * (x-X[i-1]) + Y[i-1]
230 out = out * m0 + z * (1-m0)
231 m0=whereNonPositive(x-X[i])
232
233 if self.__obeyBounds:
234 if inf(m0) < 1 :
235 raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
236 else:
237 out = out * m0 + X[-1] * (1-m0)
238 return out
239
240 def getValueDifferential(self, x):
241 X=self.__x
242 Y=self.__y
243
244 x0=X[0]
245 m0=whereNegative(x-x0)
246 if self.__obeyBounds:
247 if sup(m0) > 0:
248 raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
249 out=0.
250 for i in range(1,len(X)):
251 z=(Y[i]-Y[i-1])/(X[i]-X[i-1])
252 out = out * m0 + z * (1-m0)
253 m0=whereNonPositive(x-X[i])
254
255 if self.__obeyBounds:
256 if inf(m0) < 1:
257 raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
258 else:
259 out = out * m0
260 return out
261
262
263 class Well(object):
264 """
265 generic well
266
267 :var WATER: phase identifier for water well
268 :var GAS: phase identifier for gas well
269 """
270 WATER="Water"
271 GAS="Gas"
272 def __init__(self, name, domain, Q=[0.], schedule=[0.], phase="Water", BHP_limit=1.*U.atm, *args, **kwargs):
273 """
274 set up well
275 """
276 if not len(schedule) == len(Q):
277 raise ValueError,"length of schedule and Q must match."
278 self.__schedule=schedule
279 self.__Q = Q
280 self.__phase=phase
281 self.__BHP_limit=BHP_limit
282 self.name=name
283 self.domain=domain
284
285 def getGeometry(self):
286 """
287 returns a function representing the geometry of the well.
288 typically a Gaussian profile around the location of the well.
289
290 :note: needs to be overwritten
291 """
292 raise NotImplementedError
293
294 def getProductivityIndex(self):
295 """
296 returns the productivity index of the well.
297 typically a Gaussian profile around the location of the well.
298
299 :note: needs to be overwritten
300 """
301 raise NotImplementedError
302
303 def getFlowRate(self,t):
304 """
305 returns the flow rate
306 """
307 if t <= self.__schedule[0]:
308 out = self.__Q[0]
309 elif t>= self.__schedule[-1]:
310 out=self.__Q[-1]
311 else:
312 for i in range(1, len(self.__schedule)) :
313 if t < self.__schedule[i]:
314 out = ( self.__Q[i]-self.__Q[i-1] )/(self.__schedule[i] - self.__schedule[i-1]) * ( t - self.__schedule[i-1] ) + self.__Q[i-1]
315 break
316 return out
317
318 def getBHPLimit(self):
319 """
320 return bottom-hole pressure limit
321
322 :note: needs to be overwritten
323 """
324 return self.__BHP_limit
325
326 def getBHP(self):
327 """
328 return bottom-hole pressure
329
330 :note: needs to be overwritten
331 """
332 return self.__BHP
333
334 def setBHP(self, BHP):
335 """
336 sets bottom-hole pressure
337 """
338 self.__BHP= BHP
339
340 def getPhase(self):
341 """
342 returns the pahse the well works on
343 """
344 return self.__phase
345
346 def setWaterProduction(self, v):
347 self.water_production=v
348 def getWaterProduction(self):
349 return self.water_production
350 def setGasProduction(self, v):
351 self.gas_production=v
352 def getGasProduction(self):
353 return self.gas_production
354
355 class VerticalPeacemanWell(Well):
356 """
357 defines a well using Peaceman's formula
358 """
359 def __init__(self,name, domain, schedule = [0.], BHP_limit=1.*U.atm, Q=0, r=10*U.cm, X0=[0.,0.,0.], D=[1.*U.km,1.*U.km, 1.*U.m],
360 perm=[1.*U.cPoise,1.*U.cPoise, 1.*U.cPoise], phase=Well.WATER, s=0, decay_factor = 5):
361 # reset_factor=0.1):
362 """
363 set up well
364
365 :param BHP: ottom-hole pressure
366 :param Q: flow rate
367 :param r: well radius
368 :param X: location
369 :param D: dimension of well block
370 :param perm: permeabilities
371 :param s: skin factor
372 """
373 Well.__init__(self, name, domain, Q=Q, schedule=schedule, phase=phase,BHP_limit=BHP_limit)
374 r_el=0.28 * sqrt( sqrt(perm[1]/perm[0]) * D[0]**2 + sqrt(perm[0]/perm[1]) * D[1]**2 )\
375 / ( (perm[1]/perm[0])**0.25 + (perm[1]/perm[0])**0.25 )
376 self.__PI = 2 * math.pi * D[2] * sqrt(perm[1]*perm[0]) / (log(r_el/r) + s)
377 self.__r = r
378
379 self.__D=D
380 self.r_el=r_el
381 if USE_NODAL_WELL:
382 self.X0=Locator(ContinuousFunction(self.domain),X0[:self.domain.getDim()]).getX()
383 else:
384 self.X0=X0[:self.domain.getDim()]
385
386 if USE_NODAL_WELL:
387 x=domain.getX()
388 self.chi = 1.
389 for i in range(domain.getDim()):
390 self.chi = self.chi * exp(-((x[i]-self.X0[i])/(0.1*D[i]))**2)
391 #saveVTK("t.vtu", chi=self.chi)
392 #1/0
393 else:
394 x=Function(domain).getX()
395 self.chi = 1.
396 if decay_factor > 0:
397 for i in range(domain.getDim()):
398 r=decay_factor * D[i]
399 self.chi = self.chi * clip( (D[i]/2.+r-abs(X0[i]-x[i]))/r, maxval=1., minval=0.)
400 else:
401 for i in range(domain.getDim()):
402 r=D[i]
403 self.chi = self.chi * whereNonPositive(abs(X0[i]-x[i])-D[i]/2.)
404
405 l=integrate(self.chi)
406 self.chi*=1./l
407
408
409 #self.chi=0.00000001
410 def getLocation(self):
411 return self.X0
412
413 def getGeometry(self):
414 """
415 returns a function representing the geometry of the well.
416 typically a Gaussian profile around the location of the well.
417
418 :note: needs to be overwritten
419 """
420 if self.domain.getDim() <3:
421 return self.chi/self.__D[2]
422 else:
423 return self.chi
424
425 def getProductivityIndex(self):
426 """
427 returns the productivity index of the well.
428 """
429 return self.__PI
430
431
432 class DualPorosity(object):
433 """
434 generic dual porosity model
435 """
436 def __init__(self, domain, phi_f=None, phi_m=None, phi=None,
437 S_fg=None, S_mg=None,
438 perm_f_0=None, perm_f_1=None, perm_f_2=None,
439 perm_m_0=None, perm_m_1=None, perm_m_2=None,
440 k_w =None, k_g=None, mu_w =None, mu_g =None,
441 rho_w =None, rho_g=None,
442 wells=[], g=9.81*U.m/U.sec**2):
443 """
444 set up
445
446 :param domain: domain
447 :type domain: `esys.escript.Domain`
448 :param phi_f: porosity of the fractured rock as function of the gas pressure
449 :type phi_f: `MaterialPropertyWithDifferential`
450 :param phi_m: porosity of the coal matrix as function of the gas pressure, if None phi_m = phi-phi_f
451 :type phi_m: `MaterialPropertyWithDifferential` or None
452 :param phi: total porosity if None phi=1.
453 :type phi: scalar or None
454 :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=S_fg- p_f
455 :type S_fg: `MaterialPropertyWithDifferential`
456 :param S_mg: gas saturation in the coal matrix as function of the capillary pressure p_mc=p_mg-p_mw
457 :type S_mg: `MaterialPropertyWithDifferential`
458 :param perm_f_0: permeability in the x_0 direction in the fractured media
459 :type perm_f_0: scalar
460 :param perm_f_1: permeability in the x_1 direction in the fractured media
461 :type perm_f_1: scalar
462 :param perm_f_2: permeability in the x_2 direction in the fractured media
463 :type perm_f_2: scalar
464 :param perm_m_0: permeability in the x_0 direction in the coal matrix
465 :type perm_m_0: scalar
466 :param perm_m_1: permeability in the x_1 direction in the coal matrix
467 :type perm_m_1: scalar
468 :param perm_m_2: permeability in the x_2 direction in the coal matrix
469 :type perm_m_2: scalar
470 :param k_w: relative permeability of water as function of water saturation
471 :type k_w: `MaterialProperty`
472 :param k_g: relative permeability of gas as function of gas saturation
473 :type k_g: `MaterialProperty`
474 :param mu_w: viscosity of water as function of water pressure
475 :type mu_w: `MaterialProperty`
476 :param mu_g: viscosity of gas as function of gas pressure
477 :type mu_g: `MaterialProperty`
478 :param rho_w: density of water as function of water pressure
479 :type rho_w: `MaterialPropertyWithDifferential`
480 :param rho_g: density of gas as function of gas pressure
481 :type rho_g: `MaterialPropertyWithDifferential`
482 :param wells : list of wells
483 :type wells: list of `Well`
484 """
485
486 self.domain = domain
487 self.phi_f = phi_f
488 self.phi_m = phi_m
489 self.phi = phi
490 self.S_fg = S_fg
491 self.S_mg = S_mg
492 self.perm_f = [ perm_f_0, perm_f_1, perm_f_2 ]
493 self.perm_m = [ perm_m_0, perm_m_1, perm_m_2 ]
494 self.k_w = k_w
495 self.k_g = k_g
496 self.mu_w = mu_w
497 self.mu_g = mu_g
498 self.rho_w = rho_w
499 self.rho_g = rho_g
500 self.wells=wells
501 self.t =0
502 self.g=g
503
504 self.__iter_max=1
505 self.__rtol=1.e-4
506 self.verbose=False
507 self.XI=0.5
508 def setIterationControl(self, iter_max=None, rtol=None, verbose=None, xi=None):
509 """
510 sets parameters to control iteration process
511 """
512 if iter_max !=None: self.__iter_max=iter_max
513 if rtol !=None: self.__rtol = rtol
514 if verbose !=None: self.verbose=verbose
515 if xi !=None: self.XI=xi
516
517 def update(self, dt):
518 self.u, self.u_old = self.u.copy(), self.u
519 n=0
520 rerr=1.
521 while n < self.__iter_max and rerr > self.__rtol:
522 u=self.solvePDE(dt)
523 norm_u=Lsup(u)
524 norm_e=Lsup(u-self.u)
525
526 if norm_u > 0:
527 rerr=norm_e/norm_u
528 else:
529 rerr=norm_e
530 if self.verbose: print "iteration step %d: relative change = %e"%(n, rerr)
531 n+=1
532 self.u=u
533 print "iteration completed."
534 self.t+=dt
535
536
537 class PorosityOneHalfModel(DualPorosity):
538 """
539 Model for gas and water in double prosity model tracking water and gas
540 pressure in the fractured rock and gas concentration in the coal matrix.
541 This corresponds to the coal bed methan model in the ECLIPSE code.
542 """
543
544 def __init__(self, domain, phi_f=None, phi=None, L_g=None,
545 perm_f_0=None, perm_f_1=None, perm_f_2=None,
546 k_w =None, k_g=None, mu_w =None, mu_g =None,
547 rho_w =None, rho_g=None, sigma=0, A_mg=0, f_rg=1.,
548 wells=[], g=9.81*U.m/U.sec**2):
549 """
550 set up
551
552 :param domain: domain
553 :type domain: `esys.escript.Domain`
554 :param phi_f: porosity of the fractured rock as function of the gas pressure
555 :type phi_f: `MaterialPropertyWithDifferential`
556 :param phi: total porosity if None phi=1.
557 :type phi: scalar or None
558 :param L_g: gas adsorption as function of gas pressure
559 :type L_g: `MaterialPropertyWithDifferential`
560 :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=S_fg- p_f
561 :type S_fg: `MaterialPropertyWithDifferential`
562 :param perm_f_0: permeability in the x_0 direction in the fractured media
563 :type perm_f_0: scalar
564 :param perm_f_1: permeability in the x_1 direction in the fractured media
565 :type perm_f_1: scalar
566 :param perm_f_2: permeability in the x_2 direction in the fractured media
567 :type perm_f_2: scalar
568 :param k_w: relative permeability of water as function of water saturation
569 :type k_w: `MaterialProperty`
570 :param k_g: relative permeability of gas as function of gas saturation
571 :type k_g: `MaterialProperty`
572 :param mu_w: viscosity of water as function of water pressure
573 :type mu_w: `MaterialProperty`
574 :param mu_g: viscosity of gas as function of gas pressure
575 :type mu_g: `MaterialProperty`
576 :param rho_w: density of water as function of water pressure
577 :type rho_w: `MaterialPropertyWithDifferential`
578 :param rho_g: density of gas as function of gas pressure
579 :type rho_g: `MaterialPropertyWithDifferential`
580 :param wells : list of wells
581 :type wells: list of `Well`
582 :param sigma: shape factor for gas matrix diffusion
583 :param A_mg: diffusion constant for gas adsorption
584 :param f_rg: gas re-adsorption factor
585 """
586
587 DualPorosity.__init__(self, domain,
588 phi_f=phi_f, phi_m=None, phi=phi,
589 S_fg=None, S_mg=None,
590 perm_f_0=perm_f_0, perm_f_1=perm_f_1, perm_f_2=perm_f_2,
591 perm_m_0=None , perm_m_1=None, perm_m_2=None,
592 k_w =k_w, k_g=k_g, mu_w =mu_w, mu_g =mu_g,
593 rho_w =rho_w, rho_g=rho_g,
594 wells=wells, g=g)
595 self.L_g=L_g
596 self.sigma = sigma
597 self.A_mg = A_mg
598 self.f_rg = f_rg
599 self.__pde=LinearPDE(self.domain, numEquations=3, numSolutions =3)
600
601
602 def getPDEOptions(self):
603 """
604 returns the `SolverOptions` of the underlying PDE
605 """
606 return self.__pde.getSolverOptions()
607
608 def getState(self):
609 return self.u[0], self.u[1], self.u[2]
610
611 def getOldState(self):
612 return self.u_old[0], self.u_old[1], self.u_old[2]
613
614 def setInitialState(self, p_top=1.*U.atm, p_bottom= 1.*U.atm, S_fg=0, c_mg=None):
615 """
616 sets the initial state
617
618 :param p: pressure
619 :param S_fg: gas saturation in fractured rock
620 :param c_mg: gas concentration in coal matrix at standart conditions. if not given it is calculated
621 using the gas adsorption curve.
622 """
623 self.u=self.__pde.createSolution()
624 if self.u.getDomain().getDim() == 2:
625 self.u[0]=(p_top + p_bottom) /2
626 else:
627 z=self.u.getDomain().getX()[0]
628 z_top=sup(z)
629 z_bottom=inf(z)
630 self.u[0]=(p_bottom-p_top)/(z_bottom-z_top)*(z-z_top) + p_top
631 self.u[1]=S_fg
632 if c_mg == None:
633 self.u[2]= self.L_g(self.u[0])
634 else:
635 self.u[2]=c_mg
636
637 def solvePDE(self, dt):
638
639 C_couple=1
640
641 p_f, S_fg, c_mg=self.getState()
642 p_f_old, S_fg_old, c_mg_old=self.getOldState()
643
644 S_fw=1-S_fg
645
646 if self.verbose:
647 print "p_f range = ",inf(p_f),sup(p_f)
648 print "S_fg range = ",inf(S_fg),sup(S_fg)
649 print "S_fw range = ",inf(S_fw),sup(S_fw)
650 print "c_mg range = ",inf(c_mg),sup(c_mg)
651
652 k_fw=self.k_w(S_fw)
653 if self.verbose: print "k_fw range = ",inf(k_fw),sup(k_fw)
654
655 k_fg=self.k_g(S_fg)
656 if self.verbose: print "k_fg range = ",inf(k_fg),sup(k_fg)
657
658 mu_fw=self.mu_w(p_f)
659 if self.verbose: print "mu_fw range = ",inf(mu_fw),sup(mu_fw)
660
661 mu_fg=self.mu_g(p_f)
662 if self.verbose: print "mu_fg range = ",inf(mu_fg),sup(mu_fg)
663
664
665 phi_f =self.phi_f.getValue(p_f)
666 dphi_fdp=self.phi_f.getValueDifferential(p_f)
667 if self.verbose: print "phi_f range = ",inf(phi_f),sup(phi_f)," (slope %e,%e)"%(inf(dphi_fdp), sup(dphi_fdp))
668
669 rho_fw = self.rho_w.getValue(p_f)
670 drho_fwdp = self.rho_w.getValueDifferential(p_f)
671 if self.verbose: print "rho_fw range = ",inf(rho_fw),sup(rho_fw)," (slope %e,%e)"%(inf(drho_fwdp), sup(drho_fwdp))
672
673 rho_fg = self.rho_g.getValue(p_f)
674 rho_g_surf = self.rho_g.rho_surf
675 drho_fgdp = self.rho_g.getValueDifferential(p_f)
676 if self.verbose:
677 print "rho_fg range = ",inf(rho_fg),sup(rho_fg)," (slope %e,%e)"%(inf(drho_fgdp), sup(drho_fgdp))
678 print "rho_fg surf = ",rho_g_surf
679
680 L_g = self.L_g(p_f)
681 dL_gdp = self.L_g.getValueDifferential(p_f)
682 if self.verbose: print "L_g range = ",inf(L_g),sup(L_g)," (slope %e,%e)"%(inf(dL_gdp), sup(dL_gdp))
683
684 A_fw = rho_fw * k_fw/mu_fw
685 A_fg = rho_fg * k_fg/mu_fg
686
687 m = whereNegative(L_g - c_mg)
688 H = self.sigma * self.A_mg * (m + (1-m) * self.f_rg * S_fg )
689
690 E_fpp= S_fw * ( rho_fw * dphi_fdp + phi_f * drho_fwdp )
691 E_fps= - phi_f * rho_fw
692 E_fsp= S_fg * ( rho_fg * dphi_fdp + phi_f * drho_fgdp )
693 E_fss= phi_f * rho_fg
694
695 F_fw=0.
696 F_fg=0.
697 D_fw=0.
698 D_fg=0.
699
700 for I in self.wells:
701 chi_I= I.getGeometry()
702 Pi_I = I.getProductivityIndex()
703 BHP_limit_I=I.getBHPLimit()
704 p_f_I=integrate(p_f*I.chi)
705 S_fg_I=integrate(S_fg*I.chi)
706 if not self.t > 0:
707 print "AAA", self.t/U.day, p_f_I/U.psi, inf(p_f)/U.psi, S_fg_I, integrate(c_mg*I.chi)/U.Mscf*U.ft**3
708 A_fw_I= self.k_w(1-S_fg_I)/self.mu_w(p_f_I)
709 A_fg_I= self.k_g(S_fg_I)/self.mu_g(p_f_I)
710
711 if self.verbose: print "well %s is open."%I.name
712 if I.getPhase() == Well.WATER:
713 q_I = I.getFlowRate(self.t+dt)
714 p_I_ref=p_f_I-q_I/(A_fw_I * Pi_I)
715 else:
716 q_I = I.getFlowRate(self.t+dt)
717 p_I_ref=p_f_I-q_I/(A_fg_I * Pi_I)
718
719 if BHP_limit_I > p_I_ref:
720 BHP_I=BHP_limit_I
721 D_fw = D_fw + Pi_I * A_fw_I * chi_I
722 F_fw = F_fw - Pi_I * A_fw_I * BHP_I *chi_I
723 D_fg = D_fg + Pi_I * A_fg_I * chi_I
724 F_fg = F_fg - Pi_I * A_fg_I * BHP_I *chi_I
725 else:
726 BHP_I=p_I_ref
727 if I.getPhase() == Well.WATER:
728 F_fw = F_fw + q_I * chi_I
729 F_fg = F_fg + A_fg_I/A_fw_I * q_I * chi_I
730 else:
731 F_fg = F_fg + q_I * chi_I
732 F_fw = F_fw + A_fw_I/A_fg_I * q_I *chi_I
733
734
735 if self.verbose: print "well %s:BHP = %e"%(I.name, BHP_I)
736 I.setBHP(BHP_I)
737
738 F_fp_Y = - rho_fw * F_fw
739 F_fs_Y = - rho_fg * F_fg
740 D_fpp = rho_fw * D_fw
741 D_fsp = rho_fg * D_fg
742
743
744 if self.domain.getDim() == 3:
745 F_fp_X = ( A_fw * self.g * rho_fw * self.perm_f[2] ) * kronecker(self.domain)[2]
746 F_fs_X = ( A_fg * self.g * rho_fg * self.perm_f[2] ) * kronecker(self.domain)[2]
747 else:
748 F_fp_X = 0 * kronecker(self.domain)[1]
749 F_fs_X = 0 * kronecker(self.domain)[1]
750
751 F_mg_Y = H * L_g
752
753
754 D=self.__pde.createCoefficient("D")
755 A=self.__pde.createCoefficient("A")
756 Y=self.__pde.createCoefficient("Y")
757 X=self.__pde.createCoefficient("X")
758
759
760 D[0,0]=E_fpp + dt * D_fpp
761 D[0,1]=E_fps
762 D[1,0]=E_fsp + dt * D_fsp
763 D[1,1]=E_fss
764 #D[1,2]=rho_g_surf * C_couple
765 D[1,2]= rho_g_surf # rho_g_surf
766 D[0,2]= -dt * H * dL_gdp * 0
767 D[2,2]= 1 + dt * H
768
769
770 for i in range(self.domain.getDim()):
771 A[0,i,0,i] = dt * A_fw * self.perm_f[i]
772 A[1,i,1,i] = dt * A_fg * self.perm_f[i]
773
774 X[0,:]= dt * F_fp_X
775 X[1,:]= dt * F_fs_X
776
777 Y[0] = E_fpp * p_f_old + E_fps * S_fg_old + dt * F_fp_Y
778 Y[1] = E_fsp * p_f_old + E_fss * S_fg_old + C_couple * rho_g_surf * c_mg_old + dt * F_fs_Y
779 Y[2] = c_mg_old + dt * ( F_mg_Y - H * dL_gdp * p_f *0 )
780
781 print "HHH D[0,0] = ",D[0,0]
782 print "HHH D[0,1] = ",D[0,1]
783 print "HHH D[0,2] = ",D[0,2]
784 print "HHH D[1,0] = ",D[1,0]
785 print "HHH D[1,1] = ",D[1,1]
786 print "HHH D[1,2] = ",D[1,2]
787 print "HHH D[2,0] = ",D[2,0]
788 print "HHH D[2,1] = ",D[2,1]
789 print "HHH D[2,2] = ",D[2,2]
790 print "HHH A_fw = ",A_fw
791 print "HHH A_fg = ",A_fg
792 print "HHH A[0,0,0,0] = ",A[0,0,0,0]
793 print "HHH A[0,1,0,1] = ",A[0,1,0,1]
794 print "HHH A[1,0,1,0] = ",A[1,0,1,0]
795 print "HHH A[1,1,1,1] = ",A[1,1,1,1]
796 print "HHH Y[0] ",Y[0]
797 print "HHH Y[1] ",Y[1]
798 print "HHH Y[2] ",Y[2]
799 print "HHH F_fp_Y ",F_fp_Y
800 print "HHH F_fs_Y ",F_fs_Y
801 print "HHH F_mg_Y ",F_mg_Y
802 print "HHH H = ",H
803
804 self.__pde.setValue(A=A, D=D, X=X, Y=Y)
805
806 u2 = self.__pde.getSolution()
807 # this is deals with round-off errors to maintain physical meaningful values
808 # we should do this in a better way to detect values that are totally wrong
809 u2[1]=clip(u2[1],minval=0, maxval=1)
810 u2[2]=clip(u2[2],minval=0)
811
812 #saveVTK("t.vtu", chi = u2[0] )
813
814 print "XXX dc_mg_,t = ", (u2[2]-c_mg_old)/dt
815 print "XXX d_Sfg_,t = ", (u2[1]- S_fg_old)/dt
816
817 # update well production
818 p_f=u2[0]
819 S_fg=u2[1]
820 for I in self.wells:
821
822 Pi_I = I.getProductivityIndex()
823 q_I = I.getFlowRate(self.t+dt)
824 p_f_I=integrate(p_f*I.chi)
825 S_fg_I=integrate(S_fg*I.chi)
826 A_fw_I= self.k_w(1-S_fg_I)/self.mu_w(p_f_I)
827 A_fg_I= self.k_g(S_fg_I)/self.mu_g(p_f_I)
828 print "AAA", (self.t+dt)/U.day, p_f_I/U.psi, inf(p_f)/U.psi, S_fg_I, integrate(u2[2]*I.chi)/U.Mscf*U.ft**3
829
830 BHP_limit_I=I.getBHPLimit()
831 BHP_I=I.getBHP()
832
833 if self.verbose: print "reservoir pressure for well %s = %e gas)."%(I.name, p_f_I)
834
835 if BHP_limit_I < BHP_I:
836 if I.getPhase() == Well.WATER:
837 q_I_w = q_I
838 q_I_g = Pi_I * A_fg_I * (p_f_I - BHP_I)
839 else:
840 q_I_g = q_I
841 q_I_w = Pi_I * A_fw_I * (p_f_I - BHP_I)
842 else:
843 q_I_w = Pi_I * A_fw_I * (p_f_I - BHP_I)
844 q_I_g = Pi_I * A_fg_I * (p_f_I - BHP_I)
845 I.setWaterProduction(q_I_w)
846 I.setGasProduction(q_I_g)
847 return u2

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26