/[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 3525 - (hide annotations)
Tue May 24 04:04:58 2011 UTC (8 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 24960 byte(s)
and more typos fixed.
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 gross 3515 from esys.escript import *
27 gross 3494 from esys.weipa import saveVTK
28 gross 3484 import math
29 gross 3458
30 gross 3510 USE_NODAL_WELL = False or True
31    
32 gross 3478 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 gross 3458
57 gross 3478 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 gross 3484
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 gross 3478
88 gross 3484 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 gross 3510 rho = rho_surf (1 + X + X*X/2) with X= C * ( p - p_ref )
108 gross 3484
109 gross 3510 with rho_surf =rho_s/B_ref * gravity
110 gross 3484 """
111 gross 3510 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 gross 3484 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 gross 3510 return self.__rho_0 * (1+ X * (1+ X/2) )
123 gross 3484
124     def getValueDifferential(self, p):
125     """
126     """
127     X= self.C * ( p - self.p_ref )
128 gross 3510 return self.__rho_0 * self.C * (1+ X)
129 gross 3497
130     def getFormationFactor(self, p):
131 gross 3510 return self.rho_surf/self.getValue(p)
132 gross 3484
133 gross 3497
134    
135 gross 3484 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 gross 3510 self.rho_surf =rho_air_s * gravity
168 gross 3484 self.tab=InterpolationTable(x=p, y=B)
169    
170     def getValue(self, p):
171     """
172     returns the density for given pressure p
173     """
174 gross 3510 return self.rho_surf/self.getFormationFactor(p)
175 gross 3484
176     def getValueDifferential(self, p):
177     """
178     """
179 gross 3497 B = self.getFormationFactor(p)
180     dBdp = self.getFormationFactorDifferential(p)
181 gross 3510 return -self.rho_surf * dBdp/(B * B)
182 gross 3497
183     def getFormationFactor(self, p):
184     return self.tab.getValue(p)
185    
186     def getFormationFactorDifferential(self, p):
187     return self.tab.getValueDifferential(p)
188 gross 3484
189 gross 3497
190 gross 3484 class InterpolationTable(MaterialPropertyWithDifferential):
191 gross 3478 """
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 gross 3484 MaterialPropertyWithDifferential.__init__(self)
201 gross 3478 if len(x) != len(y):
202     raise ValueError,"length of interpolation nodes and value lists need to be identical."
203 gross 3484 if len(x) < 1 :
204 gross 3478 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 gross 3484
215 gross 3478 def getValue(self, x):
216     """
217     returns the interpolated values of x
218     """
219 gross 3486 X=self.__x
220     Y=self.__y
221    
222     x0=X[0]
223 gross 3478 m0=whereNegative(x-x0)
224     if self.__obeyBounds:
225     if sup(m0) > 0:
226 gross 3486 raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
227 gross 3478 out=self.__x[0]
228 gross 3486 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 gross 3478 out = out * m0 + z * (1-m0)
231 gross 3494 m0=whereNonPositive(x-X[i])
232    
233 gross 3478 if self.__obeyBounds:
234 gross 3494 if inf(m0) < 1 :
235 gross 3486 raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
236 gross 3478 else:
237 gross 3502 out = out * m0 + X[-1] * (1-m0)
238 gross 3478 return out
239 gross 3458
240 gross 3484 def getValueDifferential(self, x):
241 gross 3486 X=self.__x
242     Y=self.__y
243    
244     x0=X[0]
245 gross 3484 m0=whereNegative(x-x0)
246     if self.__obeyBounds:
247     if sup(m0) > 0:
248 gross 3486 raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
249 gross 3484 out=0.
250 gross 3486 for i in range(1,len(X)):
251     z=(Y[i]-Y[i-1])/(X[i]-X[i-1])
252 gross 3484 out = out * m0 + z * (1-m0)
253 gross 3494 m0=whereNonPositive(x-X[i])
254 gross 3484
255     if self.__obeyBounds:
256     if inf(m0) < 1:
257 gross 3486 raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
258 gross 3484 else:
259     out = out * m0
260     return out
261    
262    
263 gross 3478 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 gross 3524 def __init__(self, name, domain, Q=[0.], schedule=[0.], phase="Water", BHP_limit=1.*U.atm, X0=[0.,0.,0.], *args, **kwargs):
273 gross 3478 """
274     set up well
275     """
276 gross 3510 if not len(schedule) == len(Q):
277     raise ValueError,"length of schedule and Q must match."
278 gross 3484 self.__schedule=schedule
279 gross 3510 self.__Q = Q
280 gross 3484 self.__phase=phase
281 gross 3494 self.__BHP_limit=BHP_limit
282     self.name=name
283     self.domain=domain
284 gross 3525 self.locator=Locator(DiracDeltaFunctions(self.domain),X0[:self.domain.getDim()])
285 gross 3524 self.X0=self.locator.getX()
286 gross 3484
287 gross 3524 def getLocation(self):
288     return self.X0
289 gross 3458
290 gross 3478 def getProductivityIndex(self):
291     """
292     returns the productivity index of the well.
293     typically a Gaussian profile around the location of the well.
294    
295     :note: needs to be overwritten
296     """
297     raise NotImplementedError
298    
299 gross 3510 def getFlowRate(self,t):
300 gross 3478 """
301     returns the flow rate
302     """
303 gross 3515 out=0.
304     for i in range(len(self.__schedule)):
305     if t <= self.__schedule[i]:
306     out=self.__Q[i]
307     break
308 gross 3510 return out
309 gross 3494
310     def getBHPLimit(self):
311     """
312     return bottom-hole pressure limit
313    
314     :note: needs to be overwritten
315     """
316     return self.__BHP_limit
317    
318 gross 3478 def getPhase(self):
319     """
320     returns the pahse the well works on
321     """
322 gross 3484 return self.__phase
323    
324     class VerticalPeacemanWell(Well):
325 gross 3478 """
326     defines a well using Peaceman's formula
327     """
328 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],
329 gross 3510 perm=[1.*U.cPoise,1.*U.cPoise, 1.*U.cPoise], phase=Well.WATER, s=0, decay_factor = 5):
330 gross 3502 # reset_factor=0.1):
331 gross 3484 """
332     set up well
333    
334     :param BHP: ottom-hole pressure
335     :param Q: flow rate
336     :param r: well radius
337     :param X: location
338     :param D: dimension of well block
339     :param perm: permeabilities
340     :param s: skin factor
341     """
342 gross 3524 Well.__init__(self, name, domain, Q=Q, schedule=schedule, phase=phase,BHP_limit=BHP_limit, X0=X0)
343 gross 3494 r_el=0.28 * sqrt( sqrt(perm[1]/perm[0]) * D[0]**2 + sqrt(perm[0]/perm[1]) * D[1]**2 )\
344 gross 3484 / ( (perm[1]/perm[0])**0.25 + (perm[1]/perm[0])**0.25 )
345     self.__PI = 2 * math.pi * D[2] * sqrt(perm[1]*perm[0]) / (log(r_el/r) + s)
346     self.__r = r
347 gross 3458
348 gross 3484 self.__D=D
349     self.r_el=r_el
350 gross 3494
351    
352     def getProductivityIndex(self):
353 gross 3484 """
354     returns the productivity index of the well.
355     """
356     return self.__PI
357 gross 3494
358 gross 3484
359 gross 3478 class DualPorosity(object):
360     """
361     generic dual porosity model
362     """
363     def __init__(self, domain, phi_f=None, phi_m=None, phi=None,
364     S_fg=None, S_mg=None,
365     perm_f_0=None, perm_f_1=None, perm_f_2=None,
366     perm_m_0=None, perm_m_1=None, perm_m_2=None,
367     k_w =None, k_g=None, mu_w =None, mu_g =None,
368     rho_w =None, rho_g=None,
369 gross 3510 wells=[], g=9.81*U.m/U.sec**2):
370 gross 3478 """
371     set up
372    
373     :param domain: domain
374     :type domain: `esys.escript.Domain`
375     :param phi_f: porosity of the fractured rock as function of the gas pressure
376     :type phi_f: `MaterialPropertyWithDifferential`
377     :param phi_m: porosity of the coal matrix as function of the gas pressure, if None phi_m = phi-phi_f
378     :type phi_m: `MaterialPropertyWithDifferential` or None
379     :param phi: total porosity if None phi=1.
380     :type phi: scalar or None
381 gross 3484 :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=S_fg- p_f
382 gross 3478 :type S_fg: `MaterialPropertyWithDifferential`
383     :param S_mg: gas saturation in the coal matrix as function of the capillary pressure p_mc=p_mg-p_mw
384     :type S_mg: `MaterialPropertyWithDifferential`
385     :param perm_f_0: permeability in the x_0 direction in the fractured media
386     :type perm_f_0: scalar
387     :param perm_f_1: permeability in the x_1 direction in the fractured media
388     :type perm_f_1: scalar
389     :param perm_f_2: permeability in the x_2 direction in the fractured media
390     :type perm_f_2: scalar
391     :param perm_m_0: permeability in the x_0 direction in the coal matrix
392     :type perm_m_0: scalar
393     :param perm_m_1: permeability in the x_1 direction in the coal matrix
394     :type perm_m_1: scalar
395     :param perm_m_2: permeability in the x_2 direction in the coal matrix
396     :type perm_m_2: scalar
397     :param k_w: relative permeability of water as function of water saturation
398     :type k_w: `MaterialProperty`
399     :param k_g: relative permeability of gas as function of gas saturation
400     :type k_g: `MaterialProperty`
401     :param mu_w: viscosity of water as function of water pressure
402     :type mu_w: `MaterialProperty`
403     :param mu_g: viscosity of gas as function of gas pressure
404     :type mu_g: `MaterialProperty`
405     :param rho_w: density of water as function of water pressure
406     :type rho_w: `MaterialPropertyWithDifferential`
407     :param rho_g: density of gas as function of gas pressure
408     :type rho_g: `MaterialPropertyWithDifferential`
409     :param wells : list of wells
410     :type wells: list of `Well`
411     """
412 gross 3458
413 gross 3478 self.domain = domain
414     self.phi_f = phi_f
415     self.phi_m = phi_m
416     self.phi = phi
417     self.S_fg = S_fg
418     self.S_mg = S_mg
419     self.perm_f = [ perm_f_0, perm_f_1, perm_f_2 ]
420     self.perm_m = [ perm_m_0, perm_m_1, perm_m_2 ]
421     self.k_w = k_w
422     self.k_g = k_g
423     self.mu_w = mu_w
424     self.mu_g = mu_g
425     self.rho_w = rho_w
426     self.rho_g = rho_g
427     self.wells=wells
428 gross 3487 self.t =0
429 gross 3510 self.g=g
430 gross 3487
431     self.__iter_max=1
432     self.__rtol=1.e-4
433 gross 3494 self.verbose=False
434     self.XI=0.5
435     def setIterationControl(self, iter_max=None, rtol=None, verbose=None, xi=None):
436 gross 3487 """
437     sets parameters to control iteration process
438     """
439     if iter_max !=None: self.__iter_max=iter_max
440     if rtol !=None: self.__rtol = rtol
441 gross 3494 if verbose !=None: self.verbose=verbose
442     if xi !=None: self.XI=xi
443 gross 3478
444 gross 3487 def update(self, dt):
445 gross 3515 self.u, self.u_old =tuple([ v.copy() for v in self.u ] ), self.u
446 gross 3487 n=0
447 gross 3515 converged=False
448     while n < self.__iter_max and not converged:
449 gross 3487 u=self.solvePDE(dt)
450 gross 3515 if self.verbose: print "iteration step %d:"
451     converged=True
452     for i in range(len(u)):
453     if isinstance(u[i], Data):
454     norm_u=Lsup(u[i])
455     norm_e=Lsup(u[i]-self.u[i])
456     else:
457     norm_e=0.
458     norm_u=1.
459 gross 3494
460 gross 3515 if norm_u > 0:
461     rerr=norm_e/norm_u
462     else:
463     rerr=norm_e
464     if rerr>self.__rtol: converged=False
465     if self.verbose: print " comp %i: relative change = %e"%(i, rerr)
466 gross 3487 n+=1
467     self.u=u
468     print "iteration completed."
469     self.t+=dt
470 gross 3458
471    
472 gross 3478 class PorosityOneHalfModel(DualPorosity):
473     """
474     Model for gas and water in double prosity model tracking water and gas
475     pressure in the fractured rock and gas concentration in the coal matrix.
476     This corresponds to the coal bed methan model in the ECLIPSE code.
477     """
478 gross 3458
479 gross 3484 def __init__(self, domain, phi_f=None, phi=None, L_g=None,
480 gross 3478 perm_f_0=None, perm_f_1=None, perm_f_2=None,
481     k_w =None, k_g=None, mu_w =None, mu_g =None,
482 gross 3494 rho_w =None, rho_g=None, sigma=0, A_mg=0, f_rg=1.,
483 gross 3510 wells=[], g=9.81*U.m/U.sec**2):
484 gross 3478 """
485     set up
486    
487     :param domain: domain
488     :type domain: `esys.escript.Domain`
489     :param phi_f: porosity of the fractured rock as function of the gas pressure
490     :type phi_f: `MaterialPropertyWithDifferential`
491     :param phi: total porosity if None phi=1.
492     :type phi: scalar or None
493     :param L_g: gas adsorption as function of gas pressure
494     :type L_g: `MaterialPropertyWithDifferential`
495 gross 3484 :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=S_fg- p_f
496 gross 3478 :type S_fg: `MaterialPropertyWithDifferential`
497     :param perm_f_0: permeability in the x_0 direction in the fractured media
498     :type perm_f_0: scalar
499     :param perm_f_1: permeability in the x_1 direction in the fractured media
500     :type perm_f_1: scalar
501     :param perm_f_2: permeability in the x_2 direction in the fractured media
502     :type perm_f_2: scalar
503     :param k_w: relative permeability of water as function of water saturation
504     :type k_w: `MaterialProperty`
505     :param k_g: relative permeability of gas as function of gas saturation
506     :type k_g: `MaterialProperty`
507     :param mu_w: viscosity of water as function of water pressure
508     :type mu_w: `MaterialProperty`
509     :param mu_g: viscosity of gas as function of gas pressure
510     :type mu_g: `MaterialProperty`
511     :param rho_w: density of water as function of water pressure
512     :type rho_w: `MaterialPropertyWithDifferential`
513     :param rho_g: density of gas as function of gas pressure
514     :type rho_g: `MaterialPropertyWithDifferential`
515     :param wells : list of wells
516     :type wells: list of `Well`
517 gross 3494 :param sigma: shape factor for gas matrix diffusion
518     :param A_mg: diffusion constant for gas adsorption
519     :param f_rg: gas re-adsorption factor
520 gross 3478 """
521    
522     DualPorosity.__init__(self, domain,
523     phi_f=phi_f, phi_m=None, phi=phi,
524 gross 3484 S_fg=None, S_mg=None,
525 gross 3478 perm_f_0=perm_f_0, perm_f_1=perm_f_1, perm_f_2=perm_f_2,
526     perm_m_0=None , perm_m_1=None, perm_m_2=None,
527     k_w =k_w, k_g=k_g, mu_w =mu_w, mu_g =mu_g,
528     rho_w =rho_w, rho_g=rho_g,
529 gross 3510 wells=wells, g=g)
530 gross 3478 self.L_g=L_g
531 gross 3494 self.sigma = sigma
532     self.A_mg = A_mg
533     self.f_rg = f_rg
534 gross 3478 self.__pde=LinearPDE(self.domain, numEquations=3, numSolutions =3)
535    
536 gross 3458
537 gross 3478 def getPDEOptions(self):
538     """
539     returns the `SolverOptions` of the underlying PDE
540     """
541     return self.__pde.getSolverOptions()
542    
543 gross 3486 def getState(self):
544 gross 3515 return self.u
545 gross 3486
546     def getOldState(self):
547 gross 3515 return self.u_old
548 gross 3486
549 gross 3510 def setInitialState(self, p_top=1.*U.atm, p_bottom= 1.*U.atm, S_fg=0, c_mg=None):
550 gross 3484 """
551     sets the initial state
552    
553     :param p: pressure
554     :param S_fg: gas saturation in fractured rock
555 gross 3510 :param c_mg: gas concentration in coal matrix at standart conditions. if not given it is calculated
556 gross 3484 using the gas adsorption curve.
557 gross 3485 """
558 gross 3515 if self.domain.getDim() == 2:
559     p_init=Scalar((p_top + p_bottom) /2, Solution(self.domain))
560 gross 3510 else:
561     z=self.u.getDomain().getX()[0]
562     z_top=sup(z)
563     z_bottom=inf(z)
564 gross 3515 p_init=(p_bottom-p_top)/(z_bottom-z_top)*(z-z_top) + p_top
565    
566     S_fg_init=Scalar(S_fg, Solution(self.domain))
567    
568 gross 3510 if c_mg == None:
569 gross 3515 c_mg_init=self.L_g(p_init)
570 gross 3484 else:
571 gross 3515 c_mg_init=Scalar(c_mg, Solution(self.domain))
572    
573 gross 3525 q_gas=Scalar(0., DiracDeltaFunctions(self.domain))
574     q_water=Scalar(0., DiracDeltaFunctions(self.domain))
575     BHP=interpolate(p_init, DiracDeltaFunctions(self.domain))
576     wells_init={ "bhp" : BHP, "q_gas": q_gas, "q_water": q_water}
577    
578 gross 3515 self.u=(p_init,S_fg_init, c_mg_init, wells_init)
579    
580 gross 3487 def solvePDE(self, dt):
581 gross 3478
582 gross 3515 p_f, S_fg, c_mg, wells_state =self.getState()
583     p_f_old, S_fg_old, c_mg_old, wells_sate_old=self.getOldState()
584 gross 3458
585 gross 3487 S_fw=1-S_fg
586 gross 3458
587 gross 3494 if self.verbose:
588     print "p_f range = ",inf(p_f),sup(p_f)
589     print "S_fg range = ",inf(S_fg),sup(S_fg)
590     print "S_fw range = ",inf(S_fw),sup(S_fw)
591 gross 3510 print "c_mg range = ",inf(c_mg),sup(c_mg)
592 gross 3515 print "wells state =",wells_state
593 gross 3494
594 gross 3487 k_fw=self.k_w(S_fw)
595 gross 3494 if self.verbose: print "k_fw range = ",inf(k_fw),sup(k_fw)
596    
597 gross 3515
598 gross 3487 k_fg=self.k_g(S_fg)
599 gross 3494 if self.verbose: print "k_fg range = ",inf(k_fg),sup(k_fg)
600    
601     mu_fw=self.mu_w(p_f)
602     if self.verbose: print "mu_fw range = ",inf(mu_fw),sup(mu_fw)
603    
604     mu_fg=self.mu_g(p_f)
605     if self.verbose: print "mu_fg range = ",inf(mu_fg),sup(mu_fg)
606 gross 3487
607    
608 gross 3494 phi_f =self.phi_f.getValue(p_f)
609     dphi_fdp=self.phi_f.getValueDifferential(p_f)
610     if self.verbose: print "phi_f range = ",inf(phi_f),sup(phi_f)," (slope %e,%e)"%(inf(dphi_fdp), sup(dphi_fdp))
611 gross 3478
612 gross 3484 rho_fw = self.rho_w.getValue(p_f)
613 gross 3494 drho_fwdp = self.rho_w.getValueDifferential(p_f)
614     if self.verbose: print "rho_fw range = ",inf(rho_fw),sup(rho_fw)," (slope %e,%e)"%(inf(drho_fwdp), sup(drho_fwdp))
615    
616     rho_fg = self.rho_g.getValue(p_f)
617 gross 3510 rho_g_surf = self.rho_g.rho_surf
618 gross 3484 drho_fgdp = self.rho_g.getValueDifferential(p_f)
619 gross 3510 if self.verbose:
620     print "rho_fg range = ",inf(rho_fg),sup(rho_fg)," (slope %e,%e)"%(inf(drho_fgdp), sup(drho_fgdp))
621     print "rho_fg surf = ",rho_g_surf
622 gross 3502
623 gross 3510 L_g = self.L_g(p_f)
624     dL_gdp = self.L_g.getValueDifferential(p_f)
625     if self.verbose: print "L_g range = ",inf(L_g),sup(L_g)," (slope %e,%e)"%(inf(dL_gdp), sup(dL_gdp))
626    
627 gross 3478 A_fw = rho_fw * k_fw/mu_fw
628     A_fg = rho_fg * k_fg/mu_fg
629    
630 gross 3510 m = whereNegative(L_g - c_mg)
631 gross 3501 H = self.sigma * self.A_mg * (m + (1-m) * self.f_rg * S_fg )
632 gross 3478
633 gross 3484 E_fpp= S_fw * ( rho_fw * dphi_fdp + phi_f * drho_fwdp )
634     E_fps= - phi_f * rho_fw
635 gross 3510 E_fsp= S_fg * ( rho_fg * dphi_fdp + phi_f * drho_fgdp )
636 gross 3484 E_fss= phi_f * rho_fg
637 gross 3524
638    
639 gross 3478
640 gross 3494 F_fw=0.
641     F_fg=0.
642     D_fw=0.
643     D_fg=0.
644 gross 3497
645 gross 3524 prod_index=Scalar(0., DiracDeltaFunctions(self.domain))
646 gross 3515 for I in self.wells:
647 gross 3524 prod_index.setTaggedValue(I.name, I.getProductivityIndex() )
648 gross 3494
649 gross 3524 F_fp_Y = A_fw * prod_index * wells_state["bhp"]
650     F_fs_Y = A_fg * prod_index * wells_state["bhp"]
651     D_fpp = A_fw * prod_index
652     D_fsp = A_fg * prod_index
653 gross 3494
654 gross 3478
655     if self.domain.getDim() == 3:
656 gross 3494 F_fp_X = ( A_fw * self.g * rho_fw * self.perm_f[2] ) * kronecker(self.domain)[2]
657     F_fs_X = ( A_fg * self.g * rho_fg * self.perm_f[2] ) * kronecker(self.domain)[2]
658 gross 3478 else:
659 gross 3484 F_fp_X = 0 * kronecker(self.domain)[1]
660     F_fs_X = 0 * kronecker(self.domain)[1]
661 gross 3478
662 gross 3510 F_mg_Y = H * L_g
663 gross 3458
664 gross 3494
665     D=self.__pde.createCoefficient("D")
666     A=self.__pde.createCoefficient("A")
667     Y=self.__pde.createCoefficient("Y")
668     X=self.__pde.createCoefficient("X")
669 gross 3524
670     d_dirac=self.__pde.createCoefficient("d_dirac")
671     y_dirac=self.__pde.createCoefficient("y_dirac")
672 gross 3494
673 gross 3524
674    
675     D[0,0]=E_fpp
676 gross 3484 D[0,1]=E_fps
677 gross 3524 d_dirac[0,0]=dt * D_fpp
678    
679     D[1,0]=E_fsp
680 gross 3494 D[1,1]=E_fss
681 gross 3524 D[1,2]= rho_g_surf
682     d_dirac[1,0]=dt * D_fsp
683    
684 gross 3510 D[0,2]= -dt * H * dL_gdp * 0
685     D[2,2]= 1 + dt * H
686 gross 3478
687    
688     for i in range(self.domain.getDim()):
689 gross 3510 A[0,i,0,i] = dt * A_fw * self.perm_f[i]
690     A[1,i,1,i] = dt * A_fg * self.perm_f[i]
691 gross 3494
692 gross 3510 X[0,:]= dt * F_fp_X
693     X[1,:]= dt * F_fs_X
694 gross 3494
695 gross 3524 Y[0] = E_fpp * p_f_old + E_fps * S_fg_old
696     Y[1] = E_fsp * p_f_old + E_fss * S_fg_old + rho_g_surf * c_mg_old
697 gross 3515 Y[2] = c_mg_old + dt * ( F_mg_Y - H * dL_gdp * p_f *0 )
698 gross 3478
699 gross 3524 y_dirac[0] =dt * F_fp_Y
700     y_dirac[1] =dt * F_fs_Y
701    
702 gross 3494 print "HHH D[0,0] = ",D[0,0]
703     print "HHH D[0,1] = ",D[0,1]
704     print "HHH D[0,2] = ",D[0,2]
705 gross 3510 print "HHH D[1,0] = ",D[1,0]
706     print "HHH D[1,1] = ",D[1,1]
707     print "HHH D[1,2] = ",D[1,2]
708     print "HHH D[2,0] = ",D[2,0]
709     print "HHH D[2,1] = ",D[2,1]
710     print "HHH D[2,2] = ",D[2,2]
711     print "HHH A_fw = ",A_fw
712     print "HHH A_fg = ",A_fg
713     print "HHH A[0,0,0,0] = ",A[0,0,0,0]
714     print "HHH A[0,1,0,1] = ",A[0,1,0,1]
715     print "HHH A[1,0,1,0] = ",A[1,0,1,0]
716     print "HHH A[1,1,1,1] = ",A[1,1,1,1]
717     print "HHH Y[0] ",Y[0]
718     print "HHH Y[1] ",Y[1]
719     print "HHH Y[2] ",Y[2]
720     print "HHH F_fp_Y ",F_fp_Y
721     print "HHH F_fs_Y ",F_fs_Y
722     print "HHH F_mg_Y ",F_mg_Y
723     print "HHH H = ",H
724 gross 3494
725 gross 3524 self.__pde.setValue(A=A, D=D, X=X, Y=Y, d_dirac=d_dirac , y_dirac=y_dirac)
726 gross 3478
727     u2 = self.__pde.getSolution()
728 gross 3494 # this is deals with round-off errors to maintain physical meaningful values
729     # we should do this in a better way to detect values that are totally wrong
730 gross 3515 p_f=u2[0]
731     S_fg=clip(u2[1],minval=0, maxval=1)
732     c_mg=clip(u2[2],minval=0)
733    
734 gross 3510
735    
736 gross 3524 q=Scalar(0., DiracDeltaFunctions(self.domain))
737     BHP_limit=Scalar(0., DiracDeltaFunctions(self.domain))
738     water_well_mask=Scalar(0., DiracDeltaFunctions(self.domain))
739     for I in self.wells:
740     q.setTaggedValue(I.name, I.getFlowRate(self.t+dt))
741     BHP_limit.setTaggedValue(I.name, I.getBHPLimit())
742 gross 3515 if I.getPhase() == Well.WATER:
743 gross 3524 water_well_mask.setTaggedValue(I.name, 1)
744    
745     p_f_wells=interpolate(p_f, DiracDeltaFunctions(self.domain))
746     S_fg_wells=interpolate(S_fg, DiracDeltaFunctions(self.domain))
747     A_fw_wells= self.k_w(1-S_fg_wells)/self.mu_w(p_f_wells)*self.rho_w(p_f_wells)
748     A_fg_wells= self.k_g(S_fg_wells)/self.mu_g(p_f_wells)*self.rho_g(p_f_wells)
749    
750     BHP=clip( p_f_wells - q/prod_index * ( m * self.rho_w.rho_surf/A_fw_wells + (1-m)* self.rho_g.rho_surf/A_fg_wells ), minval = BHP_limit)
751 gross 3525 BHP=clip( p_f_wells - q/prod_index * self.rho_w.rho_surf/A_fw_wells, minval = BHP_limit)
752    
753     print "AAA",BHP
754 gross 3524 wells_state_new={ "bhp" : BHP, "q_gas": A_fg_wells * prod_index*(p_f_wells-BHP)/self.rho_g.rho_surf, "q_water": A_fw_wells * prod_index*(p_f_wells-BHP)/self.rho_w.rho_surf }
755 gross 3525
756 gross 3515 return (p_f,S_fg, c_mg, wells_state_new)

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26