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

Annotation of /trunk/finley/test/python/coalgas.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3501 - (hide annotations)
Wed Apr 13 04:07:53 2011 UTC (8 years, 3 months ago) by gross
File MIME type: text/x-python
File size: 25967 byte(s)


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

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26