/[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 4154 by jfenwick, Tue Jan 22 09:30:23 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]
       z=(y[i]-y[i-1])/(x[i]-x[i-1]) * (x-x[i-1]) + y[i-1]  
231        out = out * m0 + z * (1-m0)        out = out * m0 + z * (1-m0)
232        m0=whereNegative(x-x[i])        m0=whereNonPositive(x-X[i])
233              
234          if self.__obeyBounds:
235            if inf(m0) < 1 :
236               raise ValueError("interpolation argument out of range [%e, %e]"%(X[0],X[-1]))
237          else:
238            out = out * m0 + X[-1] * (1-m0)
239          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:        if self.__obeyBounds:
257          if inf(m0) < 1:          if inf(m0) < 1:
258             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]))
259        else:        else:
260          out = out * m0 + y[-1] * (1-m0)          out = out * m0
261        return out        return out  
262          
263    
264  class Well(object):  class Well(object):
265     """     """
# 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            
# Line 263  class PorosityOneHalfModel(DualPorosity) Line 501  class PorosityOneHalfModel(DualPorosity)
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
# Line 285  class PorosityOneHalfModel(DualPorosity) Line 523  class PorosityOneHalfModel(DualPorosity)
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.sigma = sigma
541         self.A_mg = A_mg
542         self.f_rg  = f_rg
543       self.__pde=LinearPDE(self.domain, numEquations=3, numSolutions =3)       self.__pde=LinearPDE(self.domain, numEquations=3, numSolutions =3)
544            
545            
# Line 305  class PorosityOneHalfModel(DualPorosity) Line 549  class PorosityOneHalfModel(DualPorosity)
549       """       """
550       return self.__pde.getSolverOptions()       return self.__pde.getSolverOptions()
551            
552        def solvePDE(self):        def getState(self):
553             return self.u
      p_fw=self.u[0]  
      p_fg=self.u[1]  
      p_fc=p_fg-p_fw  
      C_mg=self.u[3]  
         
         
      p_fw_old=self.u_old[0]  
      p_fg_old=self.u_old[1]  
      C_mg_old=self.u_old[3]  
554    
555          def getOldState(self):
556         return self.u_old
557    
558       phi_f   =self.phi_f.getValue(p_fg)        def setInitialState(self, p_top=1.*U.atm, p_bottom= 1.*U.atm, S_fg=0,  c_mg=None):
559       dphi_fdp=self.phi_f.getValueDifferential(p_fg)          """
560                sets the initial state
561       S_fg=  self.S_fg.getValue(p_fc)          
562       dS_fgdp=   self.S_fg.getValueDifferential(p_fc)          :param p: pressure
563       S_fw=  1-S_fg          :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       rho_fw     = self.rho_w.getValue(p_fw)              using the  gas adsorption curve.
566       drho_fwdp  = self.rho_w.getValueDifferential(p_fw)          """    
567       rho_fg = self.rho_g.getValue(p_fg)          if self.domain.getDim() == 2:
568       drho_fgdp  = self.rho_g.getValueDifferential(p_fg)             p_init=Scalar((p_top + p_bottom) /2, Solution(self.domain))
569                else:
570       L_g       = self.getValue(p_fg)             z=self.u.getDomain().getX()[0]
571       dL_gdp_fg =  self.rho_w.getValueDifferential(p_fg)             z_top=sup(z)
572                   z_bottom=inf(z)
573               p_init=(p_bottom-p_top)/(z_bottom-z_top)*(z-z_top) + p_top
574    
575                S_fg_init=Scalar(S_fg, Solution(self.domain))
576    
577            if c_mg == None:
578                  c_mg_init=self.L_g(p_init)
579            else:
580                  c_mg_init=Scalar(c_mg, Solution(self.domain))
581    
582                q_gas=Scalar(0., DiracDeltaFunctions(self.domain))
583                q_water=Scalar(0., DiracDeltaFunctions(self.domain))
584                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    
588          def solvePDE(self, dt):
589        
590         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    
595         if self.verbose:
596              print("p_f range = ",inf(p_f),sup(p_f))
597              print("S_fg range = ",inf(S_fg),sup(S_fg))
598              print("S_fw range = ",inf(S_fw),sup(S_fw))
599              print("c_mg range = ",inf(c_mg),sup(c_mg))
600                  print("BHP =",BHP)
601                  print("q_gas =",q_gas)
602                  print("q_water =",q_water)
603    
604             k_fw=self.k_w(S_fw)
605             if self.verbose: print("k_fw range = ",inf(k_fw),sup(k_fw))
606    
607    
608             k_fg=self.k_g(S_fg)
609             if self.verbose: print("k_fg range = ",inf(k_fg),sup(k_fg))
610    
611         mu_fw=self.mu_w(p_f)
612             if self.verbose: print("mu_fw range = ",inf(mu_fw),sup(mu_fw))
613    
614         mu_fg=self.mu_g(p_f)
615             if self.verbose: print("mu_fg range = ",inf(mu_fg),sup(mu_fg))
616            
617    
618         phi_f   =self.phi_f.getValue(p_f)
619         dphi_fdp=self.phi_f.getValueDifferential(p_f)
620             if self.verbose: print("phi_f range = ",inf(phi_f),sup(phi_f)," (slope %e,%e)"%(inf(dphi_fdp), sup(dphi_fdp)))
621        
622         rho_fw     = self.rho_w.getValue(p_f)
623         drho_fwdp  = self.rho_w.getValueDifferential(p_f)
624             if self.verbose: print("rho_fw range = ",inf(rho_fw),sup(rho_fw)," (slope %e,%e)"%(inf(drho_fwdp), sup(drho_fwdp)))
625    
626             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             if self.verbose:
630                    print("rho_fg range = ",inf(rho_fg),sup(rho_fg)," (slope %e,%e)"%(inf(drho_fgdp), sup(drho_fgdp)))
631                    print("rho_fg surf = ",rho_g_surf)
632              
633         L_g = self.L_g(p_f)
634         dL_gdp = self.L_g.getValueDifferential(p_f)
635             if self.verbose: print("L_g range = ",inf(L_g),sup(L_g)," (slope %e,%e)"%(inf(dL_gdp), sup(dL_gdp)))
636            
637       A_fw = rho_fw * k_fw/mu_fw       A_fw = rho_fw * k_fw/mu_fw
638       A_fg = rho_fg * k_fg/mu_fg       A_fg = rho_fg * k_fg/mu_fg
639            
640       m = whereNegative(L_g - C_mg)         m = whereNegative(L_g - c_mg)
641       B = self.sigma * self.a_mg * (m + (1-m) * self.f_rg * S_fg )       H = self.sigma * self.A_mg * (m + (1-m) * self.f_rg * S_fg )
642            
643             E_fpp= S_fw * (  rho_fw * dphi_fdp + phi_f  * drho_fwdp )
644       E_fww= phi_f * ( rho_fw * dS_fgdp + S_fw * drho_fwdp )       E_fps=  -  phi_f * rho_fw
645       E_fwg= rho_fw * ( S_fw * dphi_fdp - phi_f * dS_fgdp )       E_fsp= S_fg * ( rho_fg * dphi_fdp + phi_f * drho_fgdp )
646       E_fgw= - phi_f * rho_fg * dS_fgdp       E_fss= phi_f * rho_fg
647       E_fgg= S_fg * rho_fg  * dphi_fdp + phi_f * rho_fg * dS_fgdp + phi_f  * S_fg * drho_fgdp      
648        
649        
650         F_fw=0.
651         F_fg=0.
652         D_fw=0.
653         D_fg=0.
654    
655             prod_index=Scalar(0., DiracDeltaFunctions(self.domain))
656             geo_fac=Scalar(0., DiracDeltaFunctions(self.domain))
657             for I in self.wells:
658             prod_index.setTaggedValue(I.name, I.getProductivityIndex() )
659             geo_fac.setTaggedValue(I.name, I.geo_fac)
660    
661         F_fp_Y = A_fw * prod_index * BHP  * geo_fac
662         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            
      Q_w=0.  
      Q_g=0.  
      R_w=0.  
      R_g=0.  
       
      for I in self.wells:  
         chi_I= I.getGeometry()  
         Pi_I = I.getProductivityIndex()  
         P_I = I.getBHP()  
         Q_I = I.getFlowRate()  
           
         if self.getPhase == Well.WATER:  
            Q_w = Q_w + Q_I        * chi_I  
            Q_g = Q_g + Pi_I * P_I * chi_I  
            R_g = R_g + Pi_I       * chi_I  
         else:  
            Q_g = Q_g + Q_I        * chi_I  
            Q_w = Q_w + Pi_I * P_I * chi_I  
            R_w = R_w + Pi_I       * chi_I  
             
      F_fw_Y = A_fw * Q_w  
      F_fg_Y = A_fg * Q_g  
      D_fgg = - A_fg * R_g  
      D_fww = - A_fw * R_w  
       
667       if self.domain.getDim() == 3:       if self.domain.getDim() == 3:
668          F_fw_X = A_fw * self.g * rho_fw * self.perm_f[2] * kronecker(self.domain)[2]          F_fp_X = ( A_fw * self.g * rho_fw * self.perm_f[2] ) * kronecker(self.domain)[2]
669          F_fg_X = A_fg * self.g * rho_fg * self.perm_f[2] * kronecker(self.domain)[2]          F_fs_X = ( A_fg * self.g * rho_fg * self.perm_f[2] ) * kronecker(self.domain)[2]
670       else:       else:
671          F_fw_X = 0 * kronecker(self.domain)[1]          F_fp_X = 0 * kronecker(self.domain)[1]
672          F_fg_X = 0 * kronecker(self.domain)[1]          F_fs_X = 0 * kronecker(self.domain)[1]
673                    
674       F_mg = B * ( L_g - dL_gdp_fg * p_fg )           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=self.__pde.getNewCoefficient("D")       d_dirac=self.__pde.createCoefficient("d_dirac")
683       A=self.__pde.getNewCoefficient("A")       y_dirac=self.__pde.createCoefficient("y_dirac")
684       Y=self.__pde.getNewCoefficient("Y")  
685       X=self.__pde.getNewCoefficient("X")  
686                
687       dtXI = dt*self.XI       D[0,0]=E_fpp
688       dtcXI = dt*(1-self.XI)       D[0,1]=E_fps
689             d_dirac[0,0]=dt * D_fpp
      D[0,0]=E_fww + dtXI * D_fww  
      D[0,1]=E_fwg  
      D[1,0]=E_fgw  
      D[1,1]=E_fgg + dtXI * D_fgg  
      D[1,2]=rho_fg  
      D[2,1]= dtXI * B * dL_gdp_fg  
      D[2,2]= 1 + dtXI * B  
690            
691       H_fw = dtcXI * A_fw * grad(p_fw_old)       D[1,0]=E_fsp
692       H_fg = dtcXI * A_fg * grad(p_fg_old)       D[1,1]=E_fss
693         D[1,2]= rho_g_surf
694         d_dirac[1,0]=dt * D_fsp
695            
696       for i in range(self.domain.getDim()):       D[0,2]= -dt * H * dL_gdp * 0
697          A[0,i,0,i] = dtXI * A_fw * self.perm_f[i]       D[2,2]= 1 + dt * H
         A[1,i,1,i] = dtXI * A_fg * self.perm_f[i]  
           
         X[0,i]= dt * F_fw_X - H_fw * self.perm_f[i]  
         X[1,i]= dt * F_fg_X - H_fg * self.perm_f[i]  
           
      Y[0] = E_fww * p_fw_old + E_fwg * p_fg_old +                     dt * F_fw_Y - dtcXI * D_fww * p_fw_old  
      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  
      Y[2] = C_mg_old                                                + dt * F_mg_Y - dtcXI * ( dL_gdp_fg * p_fg_old + C_mg_old)  
698            
699       self.__pde.setValue(A=A, D=D, X=X, Y=Y)      
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()       u2 = self.__pde.getSolution()
740       return u2       # 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.4154

  ViewVC Help
Powered by ViewVC 1.1.26