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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

Legend:
Removed from v.3478  
changed lines
  Added in v.4576

  ViewVC Help
Powered by ViewVC 1.1.26