/[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 4657 - (hide annotations)
Thu Feb 6 06:12:20 2014 UTC (5 years, 5 months ago) by jfenwick
File MIME type: text/x-python
File size: 27193 byte(s)
I changed some files.
Updated copyright notices, added GeoComp.



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