/[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 3911 - (hide annotations)
Thu Jun 14 01:01:03 2012 UTC (7 years, 1 month ago) by jfenwick
File MIME type: text/x-python
File size: 25288 byte(s)
Copyright changes
1 gross 3494 ######################################################
2 gross 3458 #
3 jfenwick 3911 # Copyright (c) 2003-2012 by University of Queensland
4 gross 3458 # 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 jfenwick 3911 __copyright__="""Copyright (c) 2003-2012 by University of Queensland
16 gross 3458 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 jfenwick 3772 raise ValueError("length of interpolation nodes and value lists need to be identical.")
203 gross 3484 if len(x) < 1 :
204 jfenwick 3772 raise ValueError("length of interpolation nodes a list needs to at least one.")
205 gross 3478
206     x_ref=x[0]
207     for i in range(1,len(x)):
208     if x_ref >= x[i]:
209 jfenwick 3772 raise ValueError("interpolation nodes need to be given in increasing order.")
210 gross 3478 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 jfenwick 3772 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 jfenwick 3772 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 jfenwick 3772 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 jfenwick 3772 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 jfenwick 3772 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 gross 3527
346     r_el=0.11271331774384821*D[0]
347 gross 3484 self.__PI = 2 * math.pi * D[2] * sqrt(perm[1]*perm[0]) / (log(r_el/r) + s)
348     self.__r = r
349 gross 3458
350 gross 3484 self.__D=D
351     self.r_el=r_el
352 gross 3527
353     if self.domain.getDim() == 3:
354     self.geo_fac=1
355     else:
356     self.geo_fac=1./D[2]
357 gross 3494
358 gross 3527
359 gross 3494
360     def getProductivityIndex(self):
361 gross 3484 """
362     returns the productivity index of the well.
363     """
364     return self.__PI
365 gross 3494
366 gross 3484
367 gross 3478 class DualPorosity(object):
368     """
369     generic dual porosity model
370     """
371     def __init__(self, domain, phi_f=None, phi_m=None, phi=None,
372     S_fg=None, S_mg=None,
373     perm_f_0=None, perm_f_1=None, perm_f_2=None,
374     perm_m_0=None, perm_m_1=None, perm_m_2=None,
375     k_w =None, k_g=None, mu_w =None, mu_g =None,
376     rho_w =None, rho_g=None,
377 gross 3510 wells=[], g=9.81*U.m/U.sec**2):
378 gross 3478 """
379     set up
380    
381     :param domain: domain
382     :type domain: `esys.escript.Domain`
383     :param phi_f: porosity of the fractured rock as function of the gas pressure
384     :type phi_f: `MaterialPropertyWithDifferential`
385     :param phi_m: porosity of the coal matrix as function of the gas pressure, if None phi_m = phi-phi_f
386     :type phi_m: `MaterialPropertyWithDifferential` or None
387     :param phi: total porosity if None phi=1.
388     :type phi: scalar or None
389 gross 3484 :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=S_fg- p_f
390 gross 3478 :type S_fg: `MaterialPropertyWithDifferential`
391     :param S_mg: gas saturation in the coal matrix as function of the capillary pressure p_mc=p_mg-p_mw
392     :type S_mg: `MaterialPropertyWithDifferential`
393     :param perm_f_0: permeability in the x_0 direction in the fractured media
394     :type perm_f_0: scalar
395     :param perm_f_1: permeability in the x_1 direction in the fractured media
396     :type perm_f_1: scalar
397     :param perm_f_2: permeability in the x_2 direction in the fractured media
398     :type perm_f_2: scalar
399     :param perm_m_0: permeability in the x_0 direction in the coal matrix
400     :type perm_m_0: scalar
401     :param perm_m_1: permeability in the x_1 direction in the coal matrix
402     :type perm_m_1: scalar
403     :param perm_m_2: permeability in the x_2 direction in the coal matrix
404     :type perm_m_2: scalar
405     :param k_w: relative permeability of water as function of water saturation
406     :type k_w: `MaterialProperty`
407     :param k_g: relative permeability of gas as function of gas saturation
408     :type k_g: `MaterialProperty`
409     :param mu_w: viscosity of water as function of water pressure
410     :type mu_w: `MaterialProperty`
411     :param mu_g: viscosity of gas as function of gas pressure
412     :type mu_g: `MaterialProperty`
413     :param rho_w: density of water as function of water pressure
414     :type rho_w: `MaterialPropertyWithDifferential`
415     :param rho_g: density of gas as function of gas pressure
416     :type rho_g: `MaterialPropertyWithDifferential`
417     :param wells : list of wells
418     :type wells: list of `Well`
419     """
420 gross 3458
421 gross 3478 self.domain = domain
422     self.phi_f = phi_f
423     self.phi_m = phi_m
424     self.phi = phi
425     self.S_fg = S_fg
426     self.S_mg = S_mg
427     self.perm_f = [ perm_f_0, perm_f_1, perm_f_2 ]
428     self.perm_m = [ perm_m_0, perm_m_1, perm_m_2 ]
429     self.k_w = k_w
430     self.k_g = k_g
431     self.mu_w = mu_w
432     self.mu_g = mu_g
433     self.rho_w = rho_w
434     self.rho_g = rho_g
435     self.wells=wells
436 gross 3487 self.t =0
437 gross 3510 self.g=g
438 gross 3487
439     self.__iter_max=1
440     self.__rtol=1.e-4
441 gross 3494 self.verbose=False
442     self.XI=0.5
443     def setIterationControl(self, iter_max=None, rtol=None, verbose=None, xi=None):
444 gross 3487 """
445     sets parameters to control iteration process
446     """
447     if iter_max !=None: self.__iter_max=iter_max
448     if rtol !=None: self.__rtol = rtol
449 gross 3494 if verbose !=None: self.verbose=verbose
450     if xi !=None: self.XI=xi
451 gross 3478
452 gross 3487 def update(self, dt):
453 gross 3515 self.u, self.u_old =tuple([ v.copy() for v in self.u ] ), self.u
454 gross 3487 n=0
455 gross 3515 converged=False
456     while n < self.__iter_max and not converged:
457 gross 3487 u=self.solvePDE(dt)
458 jfenwick 3772 if self.verbose: print("iteration step %d:"%n)
459 gross 3515 converged=True
460     for i in range(len(u)):
461     if isinstance(u[i], Data):
462     norm_u=Lsup(u[i])
463     norm_e=Lsup(u[i]-self.u[i])
464     else:
465     norm_e=0.
466     norm_u=1.
467 gross 3494
468 gross 3515 if norm_u > 0:
469     rerr=norm_e/norm_u
470     else:
471     rerr=norm_e
472 gross 3526 if norm_e>self.__rtol * norm_u + 1.e-10: converged=False
473 jfenwick 3772 if self.verbose: print(" comp %i: change = %e (value = %e)"%(i, norm_e,norm_u))
474 gross 3487 n+=1
475     self.u=u
476 jfenwick 3772 print("iteration completed.")
477 gross 3487 self.t+=dt
478 gross 3458
479    
480 gross 3478 class PorosityOneHalfModel(DualPorosity):
481     """
482     Model for gas and water in double prosity model tracking water and gas
483     pressure in the fractured rock and gas concentration in the coal matrix.
484     This corresponds to the coal bed methan model in the ECLIPSE code.
485     """
486 gross 3458
487 gross 3484 def __init__(self, domain, phi_f=None, phi=None, L_g=None,
488 gross 3478 perm_f_0=None, perm_f_1=None, perm_f_2=None,
489     k_w =None, k_g=None, mu_w =None, mu_g =None,
490 gross 3494 rho_w =None, rho_g=None, sigma=0, A_mg=0, f_rg=1.,
491 gross 3510 wells=[], g=9.81*U.m/U.sec**2):
492 gross 3478 """
493     set up
494    
495     :param domain: domain
496     :type domain: `esys.escript.Domain`
497     :param phi_f: porosity of the fractured rock as function of the gas pressure
498     :type phi_f: `MaterialPropertyWithDifferential`
499     :param phi: total porosity if None phi=1.
500     :type phi: scalar or None
501     :param L_g: gas adsorption as function of gas pressure
502     :type L_g: `MaterialPropertyWithDifferential`
503 gross 3484 :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=S_fg- p_f
504 gross 3478 :type S_fg: `MaterialPropertyWithDifferential`
505     :param perm_f_0: permeability in the x_0 direction in the fractured media
506     :type perm_f_0: scalar
507     :param perm_f_1: permeability in the x_1 direction in the fractured media
508     :type perm_f_1: scalar
509     :param perm_f_2: permeability in the x_2 direction in the fractured media
510     :type perm_f_2: scalar
511     :param k_w: relative permeability of water as function of water saturation
512     :type k_w: `MaterialProperty`
513     :param k_g: relative permeability of gas as function of gas saturation
514     :type k_g: `MaterialProperty`
515     :param mu_w: viscosity of water as function of water pressure
516     :type mu_w: `MaterialProperty`
517     :param mu_g: viscosity of gas as function of gas pressure
518     :type mu_g: `MaterialProperty`
519     :param rho_w: density of water as function of water pressure
520     :type rho_w: `MaterialPropertyWithDifferential`
521     :param rho_g: density of gas as function of gas pressure
522     :type rho_g: `MaterialPropertyWithDifferential`
523     :param wells : list of wells
524     :type wells: list of `Well`
525 gross 3494 :param sigma: shape factor for gas matrix diffusion
526     :param A_mg: diffusion constant for gas adsorption
527     :param f_rg: gas re-adsorption factor
528 gross 3478 """
529    
530     DualPorosity.__init__(self, domain,
531     phi_f=phi_f, phi_m=None, phi=phi,
532 gross 3484 S_fg=None, S_mg=None,
533 gross 3478 perm_f_0=perm_f_0, perm_f_1=perm_f_1, perm_f_2=perm_f_2,
534     perm_m_0=None , perm_m_1=None, perm_m_2=None,
535     k_w =k_w, k_g=k_g, mu_w =mu_w, mu_g =mu_g,
536     rho_w =rho_w, rho_g=rho_g,
537 gross 3510 wells=wells, g=g)
538 gross 3478 self.L_g=L_g
539 gross 3494 self.sigma = sigma
540     self.A_mg = A_mg
541     self.f_rg = f_rg
542 gross 3478 self.__pde=LinearPDE(self.domain, numEquations=3, numSolutions =3)
543    
544 gross 3458
545 gross 3478 def getPDEOptions(self):
546     """
547     returns the `SolverOptions` of the underlying PDE
548     """
549     return self.__pde.getSolverOptions()
550    
551 gross 3486 def getState(self):
552 gross 3515 return self.u
553 gross 3486
554     def getOldState(self):
555 gross 3515 return self.u_old
556 gross 3486
557 gross 3510 def setInitialState(self, p_top=1.*U.atm, p_bottom= 1.*U.atm, S_fg=0, c_mg=None):
558 gross 3484 """
559     sets the initial state
560    
561     :param p: pressure
562     :param S_fg: gas saturation in fractured rock
563 gross 3510 :param c_mg: gas concentration in coal matrix at standart conditions. if not given it is calculated
564 gross 3484 using the gas adsorption curve.
565 gross 3485 """
566 gross 3515 if self.domain.getDim() == 2:
567     p_init=Scalar((p_top + p_bottom) /2, Solution(self.domain))
568 gross 3510 else:
569     z=self.u.getDomain().getX()[0]
570     z_top=sup(z)
571     z_bottom=inf(z)
572 gross 3515 p_init=(p_bottom-p_top)/(z_bottom-z_top)*(z-z_top) + p_top
573    
574     S_fg_init=Scalar(S_fg, Solution(self.domain))
575    
576 gross 3510 if c_mg == None:
577 gross 3515 c_mg_init=self.L_g(p_init)
578 gross 3484 else:
579 gross 3515 c_mg_init=Scalar(c_mg, Solution(self.domain))
580    
581 gross 3525 q_gas=Scalar(0., DiracDeltaFunctions(self.domain))
582     q_water=Scalar(0., DiracDeltaFunctions(self.domain))
583     BHP=interpolate(p_init, DiracDeltaFunctions(self.domain))
584    
585 gross 3526 self.u=(p_init,S_fg_init, c_mg_init, BHP, q_gas, q_water)
586 gross 3515
587 gross 3487 def solvePDE(self, dt):
588 gross 3478
589 gross 3526 p_f, S_fg, c_mg, BHP_old, q_gas_old, q_water_old =self.getState()
590     p_f_old, S_fg_old, c_mg_old, BHP, q_gas, q_water =self.getOldState()
591 gross 3458
592 gross 3487 S_fw=1-S_fg
593 gross 3458
594 gross 3494 if self.verbose:
595 jfenwick 3772 print("p_f range = ",inf(p_f),sup(p_f))
596     print("S_fg range = ",inf(S_fg),sup(S_fg))
597     print("S_fw range = ",inf(S_fw),sup(S_fw))
598     print("c_mg range = ",inf(c_mg),sup(c_mg))
599     print("BHP =",BHP)
600     print("q_gas =",q_gas)
601     print("q_water =",q_water)
602 gross 3494
603 gross 3487 k_fw=self.k_w(S_fw)
604 jfenwick 3772 if self.verbose: print("k_fw range = ",inf(k_fw),sup(k_fw))
605 gross 3494
606 gross 3515
607 gross 3487 k_fg=self.k_g(S_fg)
608 jfenwick 3772 if self.verbose: print("k_fg range = ",inf(k_fg),sup(k_fg))
609 gross 3494
610     mu_fw=self.mu_w(p_f)
611 jfenwick 3772 if self.verbose: print("mu_fw range = ",inf(mu_fw),sup(mu_fw))
612 gross 3494
613     mu_fg=self.mu_g(p_f)
614 jfenwick 3772 if self.verbose: print("mu_fg range = ",inf(mu_fg),sup(mu_fg))
615 gross 3487
616    
617 gross 3494 phi_f =self.phi_f.getValue(p_f)
618     dphi_fdp=self.phi_f.getValueDifferential(p_f)
619 jfenwick 3772 if self.verbose: print("phi_f range = ",inf(phi_f),sup(phi_f)," (slope %e,%e)"%(inf(dphi_fdp), sup(dphi_fdp)))
620 gross 3478
621 gross 3484 rho_fw = self.rho_w.getValue(p_f)
622 gross 3494 drho_fwdp = self.rho_w.getValueDifferential(p_f)
623 jfenwick 3772 if self.verbose: print("rho_fw range = ",inf(rho_fw),sup(rho_fw)," (slope %e,%e)"%(inf(drho_fwdp), sup(drho_fwdp)))
624 gross 3494
625     rho_fg = self.rho_g.getValue(p_f)
626 gross 3510 rho_g_surf = self.rho_g.rho_surf
627 gross 3484 drho_fgdp = self.rho_g.getValueDifferential(p_f)
628 gross 3510 if self.verbose:
629 jfenwick 3772 print("rho_fg range = ",inf(rho_fg),sup(rho_fg)," (slope %e,%e)"%(inf(drho_fgdp), sup(drho_fgdp)))
630     print("rho_fg surf = ",rho_g_surf)
631 gross 3502
632 gross 3510 L_g = self.L_g(p_f)
633     dL_gdp = self.L_g.getValueDifferential(p_f)
634 jfenwick 3772 if self.verbose: print("L_g range = ",inf(L_g),sup(L_g)," (slope %e,%e)"%(inf(dL_gdp), sup(dL_gdp)))
635 gross 3510
636 gross 3478 A_fw = rho_fw * k_fw/mu_fw
637     A_fg = rho_fg * k_fg/mu_fg
638    
639 gross 3510 m = whereNegative(L_g - c_mg)
640 gross 3501 H = self.sigma * self.A_mg * (m + (1-m) * self.f_rg * S_fg )
641 gross 3478
642 gross 3484 E_fpp= S_fw * ( rho_fw * dphi_fdp + phi_f * drho_fwdp )
643     E_fps= - phi_f * rho_fw
644 gross 3510 E_fsp= S_fg * ( rho_fg * dphi_fdp + phi_f * drho_fgdp )
645 gross 3484 E_fss= phi_f * rho_fg
646 gross 3524
647    
648 gross 3478
649 gross 3494 F_fw=0.
650     F_fg=0.
651     D_fw=0.
652     D_fg=0.
653 gross 3497
654 gross 3524 prod_index=Scalar(0., DiracDeltaFunctions(self.domain))
655 gross 3527 geo_fac=Scalar(0., DiracDeltaFunctions(self.domain))
656 gross 3515 for I in self.wells:
657 gross 3524 prod_index.setTaggedValue(I.name, I.getProductivityIndex() )
658 gross 3527 geo_fac.setTaggedValue(I.name, I.geo_fac)
659 gross 3494
660 gross 3527 F_fp_Y = A_fw * prod_index * BHP * geo_fac
661     F_fs_Y = A_fg * prod_index * BHP * geo_fac
662     D_fpp = A_fw * prod_index * geo_fac
663     D_fsp = A_fg * prod_index * geo_fac
664 gross 3494
665 gross 3478
666     if self.domain.getDim() == 3:
667 gross 3494 F_fp_X = ( A_fw * self.g * rho_fw * self.perm_f[2] ) * kronecker(self.domain)[2]
668     F_fs_X = ( A_fg * self.g * rho_fg * self.perm_f[2] ) * kronecker(self.domain)[2]
669 gross 3478 else:
670 gross 3484 F_fp_X = 0 * kronecker(self.domain)[1]
671     F_fs_X = 0 * kronecker(self.domain)[1]
672 gross 3478
673 gross 3510 F_mg_Y = H * L_g
674 gross 3458
675 gross 3494
676     D=self.__pde.createCoefficient("D")
677     A=self.__pde.createCoefficient("A")
678     Y=self.__pde.createCoefficient("Y")
679     X=self.__pde.createCoefficient("X")
680 gross 3524
681     d_dirac=self.__pde.createCoefficient("d_dirac")
682     y_dirac=self.__pde.createCoefficient("y_dirac")
683 gross 3494
684 gross 3524
685    
686     D[0,0]=E_fpp
687 gross 3484 D[0,1]=E_fps
688 gross 3524 d_dirac[0,0]=dt * D_fpp
689    
690     D[1,0]=E_fsp
691 gross 3494 D[1,1]=E_fss
692 gross 3524 D[1,2]= rho_g_surf
693     d_dirac[1,0]=dt * D_fsp
694    
695 gross 3510 D[0,2]= -dt * H * dL_gdp * 0
696     D[2,2]= 1 + dt * H
697 gross 3478
698    
699     for i in range(self.domain.getDim()):
700 gross 3510 A[0,i,0,i] = dt * A_fw * self.perm_f[i]
701     A[1,i,1,i] = dt * A_fg * self.perm_f[i]
702 gross 3494
703 gross 3510 X[0,:]= dt * F_fp_X
704     X[1,:]= dt * F_fs_X
705 gross 3494
706 gross 3524 Y[0] = E_fpp * p_f_old + E_fps * S_fg_old
707     Y[1] = E_fsp * p_f_old + E_fss * S_fg_old + rho_g_surf * c_mg_old
708 gross 3515 Y[2] = c_mg_old + dt * ( F_mg_Y - H * dL_gdp * p_f *0 )
709 gross 3478
710 gross 3524 y_dirac[0] =dt * F_fp_Y
711     y_dirac[1] =dt * F_fs_Y
712    
713 jfenwick 3772 print("HHH D[0,0] = ",D[0,0])
714     print("HHH D[0,1] = ",D[0,1])
715     print("HHH D[0,2] = ",D[0,2])
716     print("HHH D[1,0] = ",D[1,0])
717     print("HHH D[1,1] = ",D[1,1])
718     print("HHH D[1,2] = ",D[1,2])
719     print("HHH D[2,0] = ",D[2,0])
720     print("HHH D[2,1] = ",D[2,1])
721     print("HHH D[2,2] = ",D[2,2])
722     print("HHH A_fw = ",A_fw)
723     print("HHH A_fg = ",A_fg)
724     print("HHH A[0,0,0,0] = ",A[0,0,0,0])
725     print("HHH A[0,1,0,1] = ",A[0,1,0,1])
726     print("HHH A[1,0,1,0] = ",A[1,0,1,0])
727     print("HHH A[1,1,1,1] = ",A[1,1,1,1])
728     print("HHH Y[0] ",Y[0])
729     print("HHH Y[1] ",Y[1])
730     print("HHH Y[2] ",Y[2])
731     print("HHH F_fp_Y ",F_fp_Y)
732     print("HHH F_fs_Y ",F_fs_Y)
733     print("HHH F_mg_Y ",F_mg_Y)
734     print("HHH H = ",H)
735 gross 3494
736 gross 3524 self.__pde.setValue(A=A, D=D, X=X, Y=Y, d_dirac=d_dirac , y_dirac=y_dirac)
737 gross 3478
738     u2 = self.__pde.getSolution()
739 gross 3494 # this is deals with round-off errors to maintain physical meaningful values
740     # we should do this in a better way to detect values that are totally wrong
741 gross 3515 p_f=u2[0]
742     S_fg=clip(u2[1],minval=0, maxval=1)
743     c_mg=clip(u2[2],minval=0)
744    
745 gross 3510
746    
747 gross 3524 q=Scalar(0., DiracDeltaFunctions(self.domain))
748     BHP_limit=Scalar(0., DiracDeltaFunctions(self.domain))
749     water_well_mask=Scalar(0., DiracDeltaFunctions(self.domain))
750     for I in self.wells:
751     q.setTaggedValue(I.name, I.getFlowRate(self.t+dt))
752     BHP_limit.setTaggedValue(I.name, I.getBHPLimit())
753 gross 3515 if I.getPhase() == Well.WATER:
754 gross 3524 water_well_mask.setTaggedValue(I.name, 1)
755    
756     p_f_wells=interpolate(p_f, DiracDeltaFunctions(self.domain))
757     S_fg_wells=interpolate(S_fg, DiracDeltaFunctions(self.domain))
758     A_fw_wells= self.k_w(1-S_fg_wells)/self.mu_w(p_f_wells)*self.rho_w(p_f_wells)
759     A_fg_wells= self.k_g(S_fg_wells)/self.mu_g(p_f_wells)*self.rho_g(p_f_wells)
760    
761     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)
762 gross 3525 BHP=clip( p_f_wells - q/prod_index * self.rho_w.rho_surf/A_fw_wells, minval = BHP_limit)
763 gross 3527 q_gas = prod_index * A_fg_wells * (p_f_wells - BHP)/self.rho_g.rho_surf
764     q_water = prod_index * A_fw_wells * (p_f_wells - BHP)/self.rho_w.rho_surf
765    
766 gross 3526 return (p_f,S_fg, c_mg, BHP, q_gas, q_water )

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26