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

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

  ViewVC Help
Powered by ViewVC 1.1.26