/[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 3509 by gross, Thu Apr 28 05:06:24 2011 UTC revision 3510 by gross, Fri May 13 06:09:49 2011 UTC
# Line 23  __url__="https://launchpad.net/escript-f Line 23  __url__="https://launchpad.net/escript-f
23  from esys.escript.linearPDEs import LinearPDE  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  from esys.escript.pdetools import Locator
26  from esys.escript import sqrt, log, whereNegative, sup, inf, whereNonPositive, integrate, Function, kronecker, grad, Lsup, clip, maximum, exp  from esys.escript import sqrt, log, whereNegative, sup, inf, whereNonPositive, integrate, Function, kronecker, grad, Lsup, clip, maximum, exp, ContinuousFunction
27  from esys.weipa import saveVTK  from esys.weipa import saveVTK
28  import math  import math
29    
30    USE_NODAL_WELL = False or True
31    
32  class MaterialProperty(object):  class MaterialProperty(object):
33     """     """
34     generic class for material properties depending on one (or several) status     generic class for material properties depending on one (or several) status
# Line 102  class WaterDensity(MaterialPropertyWithD Line 104  class WaterDensity(MaterialPropertyWithD
104      """      """
105      set water density as      set water density as
106                
107            rho = rho_ref (1 + X + X*X/2) with X= C * ( p - p_ref )            rho = rho_surf (1 + X + X*X/2) with X= C * ( p - p_ref )
108                        
109      with rho_ref =rho_s/B_ref * gravity      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= 998.2 * U.kg/U.m**3):      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_0 = rho_s * gravity        self.rho_surf = rho_s * gravity
113        self.rho_ref = self.rho_0/B_ref        self.__rho_0 = self.rho_surf/B_ref
114        self.C=C        self.C=C
115        self.p_ref=p_ref        self.p_ref=p_ref
116            
# Line 117  class WaterDensity(MaterialPropertyWithD Line 119  class WaterDensity(MaterialPropertyWithD
119        returns the density for given pressure p        returns the density for given pressure p
120        """        """
121        X= self.C * ( p - self.p_ref )        X= self.C * ( p - self.p_ref )
122        return self.rho_ref * (1+ X * (1+ X/2) )          return self.__rho_0 * (1+ X * (1+ X/2) )  
123                
124      def getValueDifferential(self,  p):      def getValueDifferential(self,  p):
125        """        """
126        """        """
127        X= self.C * ( p - self.p_ref )        X= self.C * ( p - self.p_ref )
128        return self.rho_ref * self.C * (1+ X)        return self.__rho_0 * self.C * (1+ X)
129                
130      def getFormationFactor(self, p):      def getFormationFactor(self, p):
131        return self.rho_0/self.getValue(p)        return self.rho_surf/self.getValue(p)
132    
133    
134    
# Line 162  class GasDensity(MaterialPropertyWithDif Line 164  class GasDensity(MaterialPropertyWithDif
164      where B is given by an interpolation table      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):      def __init__(self, p, B, gravity=1., rho_air_s=1.2041*U.kg/U.m**3):
167        self.rho_ref =rho_air_s * gravity        self.rho_surf =rho_air_s * gravity
       self.rho_0 =rho_air_s * gravity  
168        self.tab=InterpolationTable(x=p, y=B)        self.tab=InterpolationTable(x=p, y=B)
169            
170      def getValue(self, p):      def getValue(self, p):
171        """        """
172        returns the density for given pressure p        returns the density for given pressure p
173        """        """
174        return self.rho_ref/self.getFormationFactor(p)        return self.rho_surf/self.getFormationFactor(p)
175                
176      def getValueDifferential(self,  p):      def getValueDifferential(self,  p):
177        """        """
178        """        """
179        B    = self.getFormationFactor(p)        B    = self.getFormationFactor(p)
180        dBdp = self.getFormationFactorDifferential(p)        dBdp = self.getFormationFactorDifferential(p)
181        return  -self.rho_ref * dBdp/(B * B)        return  -self.rho_surf * dBdp/(B * B)
182                
183      def getFormationFactor(self, p):      def getFormationFactor(self, p):
184        return self.tab.getValue(p)        return self.tab.getValue(p)
# Line 268  class Well(object): Line 269  class Well(object):
269     """     """
270     WATER="Water"     WATER="Water"
271     GAS="Gas"     GAS="Gas"
272     def __init__(self, name, domain, Q=0., schedule=[0.], phase="Water", BHP_limit=1.*U.atm, *args, **kwargs):     def __init__(self, name, domain, Q=[0.], schedule=[0.], phase="Water", BHP_limit=1.*U.atm, *args, **kwargs):
273         """         """
274         set up well         set up well
275         """         """
276           if not len(schedule) == len(Q):
277               raise ValueError,"length of schedule and Q must match."
278         self.__schedule=schedule         self.__schedule=schedule
        self.__phase=phase  
279         self.__Q = Q         self.__Q = Q
280           self.__phase=phase
281         self.__BHP_limit=BHP_limit         self.__BHP_limit=BHP_limit
282         self.name=name         self.name=name
283         self.domain=domain         self.domain=domain
# Line 297  class Well(object): Line 300  class Well(object):
300         """         """
301         raise NotImplementedError         raise NotImplementedError
302        
303     def getFlowRate(self):     def getFlowRate(self,t):
304        """        """
305        returns the flow rate        returns the flow rate
306        """        """
307        return self.__Q        if t <= self.__schedule[0]:
308             out = self.__Q[0]
309          elif t>= self.__schedule[-1]:
310             out=self.__Q[-1]
311          else:
312            for i in range(1, len(self.__schedule)) :
313                if t < self.__schedule[i]:
314                    out = ( self.__Q[i]-self.__Q[i-1] )/(self.__schedule[i] - self.__schedule[i-1]) * ( t - self.__schedule[i-1] ) + self.__Q[i-1]
315                    break
316          return out
317            
318     def getBHPLimit(self):     def getBHPLimit(self):
319        """        """
# Line 331  class Well(object): Line 343  class Well(object):
343        """        """
344        return self.__phase        return self.__phase
345                
    def isOpen(self, t):  
      """  
      returns True is the well is open at time t  
      """  
      out = False  
      t0=min(t, self.__schedule[0])  
      i=0  
      while t > t0:  
        if out:  
      out=False  
        else:  
      out=True    
        i+=1  
        if i < len(self.__schedule):  
       t0=self.__schedule[i]  
        else:  
       t0=t  
      return out  
         
346     def setWaterProduction(self, v):     def setWaterProduction(self, v):
347        self.water_production=v        self.water_production=v
348     def getWaterProduction(self):     def getWaterProduction(self):
# Line 364  class VerticalPeacemanWell(Well): Line 357  class VerticalPeacemanWell(Well):
357     defines a well using Peaceman's formula     defines a well using Peaceman's formula
358     """     """
359     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],     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],
360               perm=[1.*U.cPoise,1.*U.cPoise, 1.*U.cPoise], phase=Well.WATER, s=0, decay_factor = 5.):               perm=[1.*U.cPoise,1.*U.cPoise, 1.*U.cPoise], phase=Well.WATER, s=0, decay_factor = 5):
361               # reset_factor=0.1):               # reset_factor=0.1):
362         """         """
363         set up well         set up well
# Line 385  class VerticalPeacemanWell(Well): Line 378  class VerticalPeacemanWell(Well):
378    
379         self.__D=D         self.__D=D
380         self.r_el=r_el         self.r_el=r_el
381         self.X0=X0[:self.domain.getDim()]         if USE_NODAL_WELL:
382               self.X0=Locator(ContinuousFunction(self.domain),X0[:self.domain.getDim()]).getX()
383           else:
384               self.X0=X0[:self.domain.getDim()]
385                
386         x=Function(domain).getX()         if USE_NODAL_WELL:
387         self.chi = 1.             x=domain.getX()
388         decay_factor             self.chi = 1.
389         for i in range(domain.getDim()):             for i in range(domain.getDim()):
390          r=decay_factor * D[i]                   self.chi = self.chi * exp(-((x[i]-self.X0[i])/(0.1*D[i]))**2)
391          self.chi = self.chi * clip( (D[i]/2.+r-abs(X0[i]-x[i]))/r, maxval=1., minval=0.)             #saveVTK("t.vtu", chi=self.chi)
392          #if reset_factor > 0:             #1/0
393          #     f=maximum(0., abs(x[i]-X0[i])-D[i]/2*(1.-reset_factor))         else:
394          #     self.chi = self.chi * exp(-(f/(reset_factor*D[i]/2))**2/2)            x=Function(domain).getX()
395          #else:            self.chi = 1.
396              #     self.chi = self.chi * whereNegative(abs(x[i]-X0[i])-D[i]/2)            if decay_factor > 0:
397                   for i in range(domain.getDim()):
398                    r=decay_factor * D[i]
399                    self.chi = self.chi * clip( (D[i]/2.+r-abs(X0[i]-x[i]))/r, maxval=1., minval=0.)
400              else:
401                   for i in range(domain.getDim()):
402                    r=D[i]
403                    self.chi = self.chi * whereNonPositive(abs(X0[i]-x[i])-D[i]/2.)
404    
405         l=integrate(self.chi)         l=integrate(self.chi)
406         self.chi*=1./l         self.chi*=1./l
407    
408                
409         #self.chi=0.00000001         #self.chi=0.00000001
410     def getLocation(self):     def getLocation(self):
# Line 435  class DualPorosity(object): Line 439  class DualPorosity(object):
439                        perm_m_0=None, perm_m_1=None, perm_m_2=None,                        perm_m_0=None, perm_m_1=None, perm_m_2=None,
440                        k_w =None, k_g=None, mu_w =None, mu_g =None,                        k_w =None, k_g=None, mu_w =None, mu_g =None,
441                        rho_w =None, rho_g=None,                        rho_w =None, rho_g=None,
442                        wells=[]):                        wells=[], g=9.81*U.m/U.sec**2):
443        """        """
444        set up        set up
445                
# Line 495  class DualPorosity(object): Line 499  class DualPorosity(object):
499        self.rho_g = rho_g            self.rho_g = rho_g    
500        self.wells=wells        self.wells=wells
501        self.t =0        self.t =0
502          self.g=g
503                
504        self.__iter_max=1        self.__iter_max=1
505        self.__rtol=1.e-4        self.__rtol=1.e-4
# Line 540  class PorosityOneHalfModel(DualPorosity) Line 545  class PorosityOneHalfModel(DualPorosity)
545               perm_f_0=None, perm_f_1=None, perm_f_2=None,               perm_f_0=None, perm_f_1=None, perm_f_2=None,
546               k_w =None, k_g=None, mu_w =None, mu_g =None,               k_w =None, k_g=None, mu_w =None, mu_g =None,
547               rho_w =None, rho_g=None, sigma=0, A_mg=0, f_rg=1.,                 rho_w =None, rho_g=None, sigma=0, A_mg=0, f_rg=1.,  
548                 wells=[]):                 wells=[], g=9.81*U.m/U.sec**2):
549       """       """
550       set up       set up
551            
# Line 586  class PorosityOneHalfModel(DualPorosity) Line 591  class PorosityOneHalfModel(DualPorosity)
591                    perm_m_0=None , perm_m_1=None, perm_m_2=None,                    perm_m_0=None , perm_m_1=None, perm_m_2=None,
592                    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,
593                    rho_w =rho_w, rho_g=rho_g,                    rho_w =rho_w, rho_g=rho_g,
594                    wells=wells)                    wells=wells, g=g)
595       self.L_g=L_g       self.L_g=L_g
596       self.sigma = sigma       self.sigma = sigma
597       self.A_mg = A_mg       self.A_mg = A_mg
# Line 606  class PorosityOneHalfModel(DualPorosity) Line 611  class PorosityOneHalfModel(DualPorosity)
611        def getOldState(self):        def getOldState(self):
612       return self.u_old[0], self.u_old[1],  self.u_old[2]       return self.u_old[0], self.u_old[1],  self.u_old[2]
613    
614        def setInitialState(self, p=1.*U.atm, S_fg=0,  C_mg=None):        def setInitialState(self, p_top=1.*U.atm, p_bottom= 1.*U.atm, S_fg=0,  c_mg=None):
615          """          """
616          sets the initial state          sets the initial state
617                    
618          :param p: pressure          :param p: pressure
619          :param S_fg: gas saturation in fractured rock          :param S_fg: gas saturation in fractured rock
620          :param C_mg: gas concentration in coal matrix. if not given it is calculated          :param c_mg: gas concentration in coal matrix at standart conditions. if not given it is calculated
621              using the  gas adsorption curve.              using the  gas adsorption curve.
622          """              """    
623          self.u=self.__pde.createSolution()          self.u=self.__pde.createSolution()
624          self.u[0]=p          if self.u.getDomain().getDim() == 2:
625               self.u[0]=(p_top + p_bottom) /2
626            else:
627               z=self.u.getDomain().getX()[0]
628               z_top=sup(z)
629               z_bottom=inf(z)
630               self.u[0]=(p_bottom-p_top)/(z_bottom-z_top)*(z-z_top) + p_top
631          self.u[1]=S_fg          self.u[1]=S_fg
632          if C_mg == None:          if c_mg == None:
633            self.u[2]= self.L_g(p)*self.rho_g.getFormationFactor(p)            self.u[2]= self.L_g(self.u[0])
634          else:          else:
635            self.u[2]=C_mg            self.u[2]=c_mg
636                
637        def solvePDE(self, dt):        def solvePDE(self, dt):
638            
639       C_couple=1       C_couple=1
640            
641       p_f, S_fg, C_mg=self.getState()       p_f, S_fg, c_mg=self.getState()
642       p_f_old, S_fg_old, C_mg_old=self.getOldState()       p_f_old, S_fg_old, c_mg_old=self.getOldState()
643    
644           S_fw=1-S_fg           S_fw=1-S_fg
645    
# Line 636  class PorosityOneHalfModel(DualPorosity) Line 647  class PorosityOneHalfModel(DualPorosity)
647            print "p_f range = ",inf(p_f),sup(p_f)            print "p_f range = ",inf(p_f),sup(p_f)
648            print "S_fg range = ",inf(S_fg),sup(S_fg)            print "S_fg range = ",inf(S_fg),sup(S_fg)
649            print "S_fw range = ",inf(S_fw),sup(S_fw)            print "S_fw range = ",inf(S_fw),sup(S_fw)
650            print "C_mg range = ",inf(C_mg),sup(C_mg)            print "c_mg range = ",inf(c_mg),sup(c_mg)
651    
652           k_fw=self.k_w(S_fw)           k_fw=self.k_w(S_fw)
653           if self.verbose: print "k_fw range = ",inf(k_fw),sup(k_fw)           if self.verbose: print "k_fw range = ",inf(k_fw),sup(k_fw)
# Line 660  class PorosityOneHalfModel(DualPorosity) Line 671  class PorosityOneHalfModel(DualPorosity)
671           if self.verbose: print "rho_fw range = ",inf(rho_fw),sup(rho_fw)," (slope %e,%e)"%(inf(drho_fwdp), sup(drho_fwdp))           if self.verbose: print "rho_fw range = ",inf(rho_fw),sup(rho_fw)," (slope %e,%e)"%(inf(drho_fwdp), sup(drho_fwdp))
672    
673           rho_fg = self.rho_g.getValue(p_f)           rho_fg = self.rho_g.getValue(p_f)
674             rho_g_surf = self.rho_g.rho_surf
675       drho_fgdp  = self.rho_g.getValueDifferential(p_f)       drho_fgdp  = self.rho_g.getValueDifferential(p_f)
          if self.verbose: print "rho_fg range = ",inf(rho_fg),sup(rho_fg)," (slope %e,%e)"%(inf(drho_fgdp), sup(drho_fgdp))  
             
   
   
      L_g_0       = self.L_g(p_f)  
      FF_g=self.rho_g.getFormationFactor(p_f)  
      L_g = L_g_0 * FF_g  
   
676           if self.verbose:           if self.verbose:
677            print "L_g_0 range = ",inf(L_g_0),sup(L_g_0)                  print "rho_fg range = ",inf(rho_fg),sup(rho_fg)," (slope %e,%e)"%(inf(drho_fgdp), sup(drho_fgdp))
678            print "FF_g range = ",inf(FF_g),sup(FF_g)                  print "rho_fg surf = ",rho_g_surf
679            print "L_g range = ",inf(L_g),sup(L_g)            
680             L_g = self.L_g(p_f)
681             dL_gdp = self.L_g.getValueDifferential(p_f)
682             if self.verbose: print "L_g range = ",inf(L_g),sup(L_g)," (slope %e,%e)"%(inf(dL_gdp), sup(dL_gdp))
683            
684       A_fw = rho_fw * k_fw/mu_fw       A_fw = rho_fw * k_fw/mu_fw
685       A_fg = rho_fg * k_fg/mu_fg       A_fg = rho_fg * k_fg/mu_fg
686            
687       m = whereNegative(L_g - C_mg)       m = whereNegative(L_g - c_mg)
688       H = 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 )
      print "XXXX H =",H  
      print "XXXX self.sigma  =",self.sigma  
      print "XXXX self.A_mg  =",self.A_mg  
689            
690       E_fpp= S_fw * (  rho_fw * dphi_fdp + phi_f  * drho_fwdp )       E_fpp= S_fw * (  rho_fw * dphi_fdp + phi_f  * drho_fwdp )
691       E_fps=  -  phi_f * rho_fw       E_fps=  -  phi_f * rho_fw
692       E_fsp= S_fg * rho_fg * dphi_fdp + (phi_f * S_fg + C_mg) * drho_fgdp       E_fsp= S_fg * ( rho_fg * dphi_fdp + phi_f * drho_fgdp )
693       E_fss= phi_f * rho_fg       E_fss= phi_f * rho_fg
694            
695       F_fw=0.       F_fw=0.
# Line 700  class PorosityOneHalfModel(DualPorosity) Line 703  class PorosityOneHalfModel(DualPorosity)
703            BHP_limit_I=I.getBHPLimit()            BHP_limit_I=I.getBHPLimit()
704            p_f_I=integrate(p_f*I.chi)            p_f_I=integrate(p_f*I.chi)
705            S_fg_I=integrate(S_fg*I.chi)            S_fg_I=integrate(S_fg*I.chi)
706              if not self.t > 0:
707                   print "AAA", self.t/U.day, p_f_I/U.psi, inf(p_f)/U.psi, S_fg_I, integrate(c_mg*I.chi)/U.Mscf*U.ft**3      
708              A_fw_I= self.k_w(1-S_fg_I)/self.mu_w(p_f_I)
709              A_fg_I= self.k_g(S_fg_I)/self.mu_g(p_f_I)
710                        
711            A_fw_I= self.rho_w.getValue(p_f_I) * self.k_w(1-S_fg_I)/self.mu_w(p_f_I)            if self.verbose: print "well %s is open."%I.name
712            A_fg_I= self.rho_g.getValue(p_f_I) * self.k_g(S_fg_I)/self.mu_g(p_f_I)                if I.getPhase() == Well.WATER:
713                          q_I = I.getFlowRate(self.t+dt)
           if I.isOpen(self.t+dt):  
         if self.verbose: print "well %s is open."%I.name  
         if I.getPhase() == Well.WATER:  
             q_I = self.rho_w.rho_0 * I.getFlowRate()  
714              p_I_ref=p_f_I-q_I/(A_fw_I * Pi_I)              p_I_ref=p_f_I-q_I/(A_fw_I * Pi_I)
715          else:            else:
716              q_I = self.rho_g.rho_0 * I.getFlowRate()              q_I = I.getFlowRate(self.t+dt)
717              p_I_ref=p_f_I-q_I/(A_fg_I * Pi_I)              p_I_ref=p_f_I-q_I/(A_fg_I * Pi_I)
               
         print "ZZZ =",p_f_I, q_I, self.rho_w.rho_0, I.getFlowRate(), A_fw_I, Pi_I, q_I/(A_fw_I * Pi_I)  
718    
719          if BHP_limit_I > p_I_ref:            if BHP_limit_I > p_I_ref:
720            BHP_I=BHP_limit_I            BHP_I=BHP_limit_I
721            D_fw = D_fw + Pi_I * A_fw_I *        chi_I            D_fw = D_fw + Pi_I * A_fw_I *        chi_I
722            F_fw = F_fw - Pi_I * A_fw_I * BHP_I *chi_I            F_fw = F_fw - Pi_I * A_fw_I * BHP_I *chi_I
723            D_fg = D_fg + Pi_I * A_fg_I *        chi_I            D_fg = D_fg + Pi_I * A_fg_I *        chi_I
724            F_fg = F_fg - Pi_I * A_fg_I * BHP_I *chi_I            F_fg = F_fg - Pi_I * A_fg_I * BHP_I *chi_I
725          else:            else:
726            BHP_I=p_I_ref            BHP_I=p_I_ref
727            if I.getPhase() == Well.WATER:            if I.getPhase() == Well.WATER:
728                F_fw = F_fw +  q_I  * chi_I                F_fw = F_fw +  q_I  * chi_I
729                D_fg = D_fg + Pi_I * A_fg_I *         chi_I                F_fg = F_fg + A_fg_I/A_fw_I * q_I * chi_I
               F_fg = F_fg - Pi_I * A_fg_I * BHP_I * chi_I  
730            else:            else:
731                F_fg = F_fg +  q_I  * chi_I                F_fg = F_fg +  q_I  * chi_I
732                D_fw = D_fw + Pi_I * A_fw_I *        chi_I                F_fw = F_fw + A_fw_I/A_fg_I * q_I *chi_I
               F_fw = F_fw - Pi_I * A_fw_I * BHP_I *chi_I  
733                        
           else:  
           if self.verbose: print "well %s is closed."%I.name  
           BHP_I=p_f_I  
734                        
735            if self.verbose: print "well %s:BHP = %e"%(I.name, BHP_I)            if self.verbose: print "well %s:BHP = %e"%(I.name, BHP_I)
736            I.setBHP(BHP_I)            I.setBHP(BHP_I)
737    
738       F_fp_Y = - F_fw       F_fp_Y = - rho_fw * F_fw
739       F_fs_Y = - F_fg       F_fs_Y = - rho_fg * F_fg
740       D_fpp =  D_fw       D_fpp =  rho_fw * D_fw
741       D_fsp =  D_fg       D_fsp =  rho_fg * D_fg
742    
743            
744       if self.domain.getDim() == 3:       if self.domain.getDim() == 3:
# Line 752  class PorosityOneHalfModel(DualPorosity) Line 748  class PorosityOneHalfModel(DualPorosity)
748          F_fp_X = 0 * kronecker(self.domain)[1]          F_fp_X = 0 * kronecker(self.domain)[1]
749          F_fs_X = 0 * kronecker(self.domain)[1]          F_fs_X = 0 * kronecker(self.domain)[1]
750                    
751           F_mg_Y = -H * L_g           F_mg_Y = H * L_g
752    
753    
754       D=self.__pde.createCoefficient("D")       D=self.__pde.createCoefficient("D")
755       A=self.__pde.createCoefficient("A")       A=self.__pde.createCoefficient("A")
756       Y=self.__pde.createCoefficient("Y")       Y=self.__pde.createCoefficient("Y")
757       X=self.__pde.createCoefficient("X")       X=self.__pde.createCoefficient("X")
758            
      dtXI = dt*self.XI  
      dtcXI = dt*(1-self.XI)  
759    
760       D[0,0]=E_fpp + dtXI * D_fpp       D[0,0]=E_fpp + dt * D_fpp
761       D[0,1]=E_fps       D[0,1]=E_fps
762       D[1,0]=E_fsp + dtXI * D_fsp       D[1,0]=E_fsp + dt * D_fsp
763       D[1,1]=E_fss       D[1,1]=E_fss
764       D[1,2]=rho_fg * C_couple       #D[1,2]=rho_g_surf * C_couple
765       #D[2,0]= - dtXI * B * dL_gdp       D[1,2]= rho_g_surf # rho_g_surf
766       D[2,2]= 1 - dtXI * H       D[0,2]= -dt * H * dL_gdp * 0
767         D[2,2]= 1 + dt * H
768            
769            
770       for i in range(self.domain.getDim()):       for i in range(self.domain.getDim()):
771          A[0,i,0,i] = dtXI * A_fw * self.perm_f[i]          A[0,i,0,i] = dt * A_fw * self.perm_f[i]
772          A[1,i,1,i] = dtXI * A_fg * self.perm_f[i]          A[1,i,1,i] = dt * A_fg * self.perm_f[i]
773    
774             X[0,:]=  dt * F_fp_X
775         X[1,:]=  dt * F_fs_X
776    
777       g= grad(p_f_old) * dtcXI * self.perm_f[0:self.domain.getDim()]       Y[0] = E_fpp *  p_f_old + E_fps * S_fg_old +                                  dt * F_fp_Y
778           X[0,:]=  - A_fw * g  + dt * F_fp_X       Y[1] = E_fsp *  p_f_old + E_fss * S_fg_old + C_couple * rho_g_surf * c_mg_old + dt * F_fs_Y
779       X[1,:]=  - A_fg * g  + dt * F_fs_X       Y[2] = c_mg_old                                                             + dt * ( F_mg_Y -  H * dL_gdp * p_f *0 )
780        
      Y[0] = E_fpp *  p_f_old + E_fps * S_fg_old +                                dt * F_fp_Y - dtcXI * D_fpp * p_f_old  
      Y[1] = E_fsp *  p_f_old + E_fss * S_fg_old + C_couple * rho_fg * C_mg_old + dt * F_fs_Y - dtcXI * D_fsp * p_f_old  
      Y[2] = C_mg_old                                                           + dt * F_mg_Y + dtcXI * H * C_mg_old  
       
      print "HHH Y[0] =", Y[0]  
      print "HHH A_fw = ",A_fw  
      print "HHH A_fg= ",A_fg  
      print "HHH F_fg = ",F_fg  
      print "HHH F_fw = ",F_fw  
      print "HHH perm_f = ",self.perm_f  
      print "HHH k = ",(A_fw*self.perm_f[0])/D[0,0]  
      print "HHH k = ",(A_fw*self.perm_f[1])/D[0,0]  
      print "HHH X[0,:]= ",X[0,:]  
781       print "HHH D[0,0] = ",D[0,0]       print "HHH D[0,0] = ",D[0,0]
782       print "HHH D[0,1] = ",D[0,1]       print "HHH D[0,1] = ",D[0,1]
783       print "HHH D[0,2] = ",D[0,2]       print "HHH D[0,2] = ",D[0,2]
784         print "HHH D[1,0] = ",D[1,0]
785         print "HHH D[1,1] = ",D[1,1]
786         print "HHH D[1,2] = ",D[1,2]
787         print "HHH D[2,0] = ",D[2,0]
788         print "HHH D[2,1] = ",D[2,1]
789         print "HHH D[2,2] = ",D[2,2]
790         print "HHH A_fw = ",A_fw
791         print "HHH A_fg = ",A_fg
792         print "HHH A[0,0,0,0] = ",A[0,0,0,0]
793         print "HHH A[0,1,0,1] = ",A[0,1,0,1]
794         print "HHH A[1,0,1,0] = ",A[1,0,1,0]
795         print "HHH A[1,1,1,1] = ",A[1,1,1,1]
796         print "HHH Y[0] ",Y[0]
797         print "HHH Y[1] ",Y[1]
798         print "HHH Y[2] ",Y[2]
799             print "HHH F_fp_Y ",F_fp_Y
800             print "HHH F_fs_Y ",F_fs_Y
801             print "HHH F_mg_Y ",F_mg_Y
802         print "HHH H = ",H
803    
804       self.__pde.setValue(A=A, D=D, X=X, Y=Y)       self.__pde.setValue(A=A, D=D, X=X, Y=Y)
805            
# Line 804  class PorosityOneHalfModel(DualPorosity) Line 808  class PorosityOneHalfModel(DualPorosity)
808       # we should do this in a better way to detect values that are totally wrong       # we should do this in a better way to detect values that are totally wrong
809       u2[1]=clip(u2[1],minval=0, maxval=1)       u2[1]=clip(u2[1],minval=0, maxval=1)
810       u2[2]=clip(u2[2],minval=0)       u2[2]=clip(u2[2],minval=0)
811    
812             #saveVTK("t.vtu", chi = u2[0] )
813    
814             print "XXX dc_mg_,t = ", (u2[2]-c_mg_old)/dt
815             print "XXX d_Sfg_,t = ", (u2[1]- S_fg_old)/dt
816            
817       # update well production       # update well production
818       p_f=u2[0]       p_f=u2[0]
819       S_fg=u2[1]       S_fg=u2[1]
820       for I in self.wells:       for I in self.wells:
821    
822            Pi_I = I.getProductivityIndex()            Pi_I = I.getProductivityIndex()
823            q_I = I.getFlowRate()            q_I = I.getFlowRate(self.t+dt)
824            p_f_I=integrate(p_f*I.chi)            p_f_I=integrate(p_f*I.chi)
825            S_fg_I=integrate(S_fg*I.chi)            S_fg_I=integrate(S_fg*I.chi)
826            rho_fw_I = self.rho_w.getValue(p_f_I)            A_fw_I= self.k_w(1-S_fg_I)/self.mu_w(p_f_I)
827            rho_fg_I = self.rho_g.getValue(p_f_I)                  A_fg_I= self.k_g(S_fg_I)/self.mu_g(p_f_I)
828            A_fw_I= rho_fw_I * self.k_w(1-S_fg_I)/self.mu_w(p_f_I)            print "AAA", (self.t+dt)/U.day, p_f_I/U.psi, inf(p_f)/U.psi, S_fg_I, integrate(u2[2]*I.chi)/U.Mscf*U.ft**3
           A_fg_I= rho_fg_I * self.k_g(S_fg_I)/self.mu_g(p_f_I)  
829    
830                BHP_limit_I=I.getBHPLimit()                BHP_limit_I=I.getBHPLimit()
831            BHP_I=I.getBHP()            BHP_I=I.getBHP()
832                        
833              if self.verbose: print "reservoir pressure for well %s = %e gas)."%(I.name, p_f_I)
834                        
835              if BHP_limit_I < BHP_I:
           if self.verbose: print "reservoir pressure for well %s = %e gas rho= %e)."%(I.name, p_f_I, rho_fg_I)  
             
           if I.isOpen(self.t+dt):  
         if BHP_limit_I < BHP_I:  
836            if I.getPhase() == Well.WATER:            if I.getPhase() == Well.WATER:
837                q_I_w =  q_I                q_I_w =  q_I
838                q_I_g = Pi_I * A_fg_I * (p_f_I - BHP_I)  /self.rho_g.rho_0                q_I_g = Pi_I * A_fg_I * (p_f_I - BHP_I)  
839            else:            else:
840                q_I_g =  q_I                q_I_g =  q_I
841                q_I_w = Pi_I * A_fw_I * (p_f_I - BHP_I)/ self.rho_w.rho_0                q_I_w = Pi_I * A_fw_I * (p_f_I - BHP_I)
         else:  
             q_I_w = Pi_I * A_fw_I * (p_f_I - BHP_I)/ self.rho_w.rho_0  
             q_I_g = Pi_I * A_fg_I * (p_f_I - BHP_I)  /self.rho_g.rho_0  
842            else:            else:
843            q_I_g =0              q_I_w = Pi_I * A_fw_I * (p_f_I - BHP_I)
844            q_I_w =0              q_I_g = Pi_I * A_fg_I * (p_f_I - BHP_I)
845            I.setWaterProduction(q_I_w)            I.setWaterProduction(q_I_w)
846            I.setGasProduction(q_I_g)            I.setGasProduction(q_I_g)
847       return u2       return u2

Legend:
Removed from v.3509  
changed lines
  Added in v.3510

  ViewVC Help
Powered by ViewVC 1.1.26