/[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 3497 by gross, Wed Apr 6 06:11:56 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 sqrt, log, whereNegative, sup, inf, whereNonPositive, integrate, Function, kronecker, grad, Lsup, clip
27    from esys.weipa import saveVTK
28    import math
29    
30  class MaterialProperty(object):  class MaterialProperty(object):
31     """     """
# Line 61  class MaterialPropertyWithDifferential(M Line 65  class MaterialPropertyWithDifferential(M
65        :remark: needs to be overwritten        :remark: needs to be overwritten
66        """        """
67        raise NotImplementedError        raise NotImplementedError
68    
69    class Porosity(MaterialPropertyWithDifferential):
70        """
71        defines porosity phi as function of pressure
72        
73             phi = phi_p_ref /(1 + X + X**2/2 ) with X= C * (p-p_ref)
74            
75        phi_p_ref is claculted from the initial porosity phi_0 at  pressure p_0
76        """
77        def __init__(self, phi_0, p_0, p_ref=1.*U.atm , C=4.0e-5/U.bar):
78          """
79          """
80          self.phi_p_ref=1 # will be overwritten later
81          self.p_ref=p_ref
82          self.C=C
83          # reset  phi_p_ref to get phi(p_0)=phi_0    
84          self.phi_p_ref=phi_0/self.getValue(p_0)
85          
86        def getValue(self, p):
87          """
88          returns the porosity for given pressure p
89          """
90          X= self.C * ( p - self.p_ref )      
91          return  self.phi_p_ref * (1. + X * (1. + X/2 ) )
92          
93        def getValueDifferential(self, p):
94          """
95          returns the porosity for given pressure p
96          """
97          X= self.C * ( p - self.p_ref )      
98          return  self.phi_p_ref * self.C * (1. + X )
99              
100          
101    class WaterDensity(MaterialPropertyWithDifferential):
102        """
103        set water density as
104          
105              rho = rho_ref (1 + X + X*X/2) with X= C * ( p - p_ref )
106              
107        with rho_ref =rho_s/B_ref * gravity
108        """
109        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):
110          self.rho_0 = rho_s * gravity
111          self.rho_ref = self.rho_0/B_ref
112          self.C=C
113          self.p_ref=p_ref
114        
115        def getValue(self, p):
116          """
117          returns the density for given pressure p
118          """
119          X= self.C * ( p - self.p_ref )
120          return self.rho_ref * (1+ X * (1+ X/2) )  
121          
122        def getValueDifferential(self,  p):
123          """
124          """
125          X= self.C * ( p - self.p_ref )
126          return self.rho_ref * self.C * (1+ X)
127          
128        def getFormationFactor(self, p):
129          return self.rho_0/self.getValue(p)
130    
131    
132    
133    class WaterViscosity(MaterialProperty):
134      """
135      defines viscosity of water
136      
137      mu=mu_ref /(1 + X + X*X/2)
138      
139      with X=-C*(p-p_ref)
140      """
141    
142      def __init__(self, mu_ref = 0.3*U.cPoise, p_ref = 1.*U.atm , C = 0./U.bar):
143        """
144        set up
145        """
146        self.mu_ref = mu_ref
147        self.p_ref = p_ref
148        self.C=C
149      def getValue(self, p):
150          """
151          returns the viscosity for given pressure p
152          """
153          X= -self.C * ( p - self.p_ref )
154          return self.mu_ref/(1+ X * (1+ X/2) )  
155          
156    class GasDensity(MaterialPropertyWithDifferential):
157        """
158        set water density as
159          
160              rho = gravity * rho_air_s /B(p)
161              
162        where B is given by an interpolation table
163        """
164        def __init__(self, p, B, gravity=1., rho_air_s=1.2041*U.kg/U.m**3):
165          self.rho_ref =rho_air_s * gravity
166          self.rho_0 =rho_air_s * gravity
167          self.tab=InterpolationTable(x=p, y=B)
168        
169        def getValue(self, p):
170          """
171          returns the density for given pressure p
172          """
173          return self.rho_ref/self.getFormationFactor(p)
174          
175        def getValueDifferential(self,  p):
176          """
177          """
178          B    = self.getFormationFactor(p)
179          dBdp = self.getFormationFactorDifferential(p)
180          return  -self.rho_ref * dBdp/(B * B)
181          
182        def getFormationFactor(self, p):
183          return self.tab.getValue(p)
184                
185  class InterpolationTable(object):      def getFormationFactorDifferential(self, p):
186          return self.tab.getValueDifferential(p)
187    
188    
189    class InterpolationTable(MaterialPropertyWithDifferential):
190     """     """
191     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
192     """     """
# Line 72  class InterpolationTable(object): Line 196  class InterpolationTable(object):
196        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
197        the value for x is set to y[0] for        the value for x is set to y[0] for
198        """        """
199          MaterialPropertyWithDifferential.__init__(self)
200        if len(x) != len(y):        if len(x) != len(y):
201       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."
202        if len(x) > 0:        if len(x) < 1 :
203       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."
204                
205        x_ref=x[0]        x_ref=x[0]
# Line 85  class InterpolationTable(object): Line 210  class InterpolationTable(object):
210        self.__x = x        self.__x = x
211        self.__y = y        self.__y = y
212        self.__obeyBounds = obeyBounds        self.__obeyBounds = obeyBounds
213          
214     def getValue(self, x):     def getValue(self, x):
215        """        """
216        returns the interpolated values of x        returns the interpolated values of x
217        """        """
218        x0=self.__x[0]        X=self.__x
219          Y=self.__y
220          
221          x0=X[0]
222        m0=whereNegative(x-x0)        m0=whereNegative(x-x0)
223        if self.__obeyBounds:        if self.__obeyBounds:
224       if sup(m0) > 0:       if sup(m0) > 0:
225          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])
226        out=self.__x[0]        out=self.__x[0]
227        for i in range(1,len(self.__x)):        for i in range(1,len(X)):
228        x1=self.__x[i]        z=(Y[i]-Y[i-1])/(X[i]-X[i-1]) * (x-X[i-1]) + Y[i-1]
229        z=(y[i]-y[i-1])/(x[i]-x[i-1]) * (x-x[i-1]) + y[i-1]        out = out * m0 + z * (1-m0)
230          m0=whereNonPositive(x-X[i])
231              
232          if self.__obeyBounds:
233            if inf(m0) < 1 :
234               raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
235          else:
236            out = out * m0 + V[-1] * (1-m0)
237          return out
238    
239       def getValueDifferential(self, x):
240          X=self.__x
241          Y=self.__y
242    
243          x0=X[0]
244          m0=whereNegative(x-x0)
245          if self.__obeyBounds:
246         if sup(m0) > 0:
247            raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
248          out=0.
249          for i in range(1,len(X)):
250          z=(Y[i]-Y[i-1])/(X[i]-X[i-1])
251        out = out * m0 + z * (1-m0)        out = out * m0 + z * (1-m0)
252        m0=whereNegative(x-x[i])        m0=whereNonPositive(x-X[i])
253            
254        if self.__obeyBounds:        if self.__obeyBounds:
255          if inf(m0) < 1:          if inf(m0) < 1:
256             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])
257        else:        else:
258          out = out * m0 + y[-1] * (1-m0)          out = out * m0
259        return out        return out  
260          
261    
262  class Well(object):  class Well(object):
263     """     """
# Line 117  class Well(object): Line 268  class Well(object):
268     """     """
269     WATER="Water"     WATER="Water"
270     GAS="Gas"     GAS="Gas"
271     def __init__(self,  *args, **kwargs):     def __init__(self, name, domain, Q=0., schedule=[0.], phase="Water", BHP_limit=1.*U.atm, *args, **kwargs):
272         """         """
273         set up well         set up well
274         """         """
275         pass         self.__schedule=schedule
276           self.__phase=phase
277           self.__Q = Q
278           self.__BHP_limit=BHP_limit
279           self.name=name
280           self.domain=domain
281          
282     def getGeometry(self):     def getGeometry(self):
283         """         """
284         returns a function representing the geometry of the well.         returns a function representing the geometry of the well.
# Line 144  class Well(object): Line 301  class Well(object):
301        """        """
302        returns the flow rate        returns the flow rate
303        """        """
304        return 0.        return self.__Q
305        
306       def getBHPLimit(self):
307          """
308          return bottom-hole pressure limit
309          
310          :note: needs to be overwritten
311          """
312          return self.__BHP_limit
313                
314     def getBHP(self):     def getBHP(self):
315        """        """
# Line 152  class Well(object): Line 317  class Well(object):
317                
318        :note: needs to be overwritten        :note: needs to be overwritten
319        """        """
320        raise NotImplementedError        return self.__BHP
321      
322       def setBHP(self, BHP):
323          """
324          sets bottom-hole pressure
325          """
326          self.__BHP= BHP
327                
328     def getPhase(self):     def getPhase(self):
329        """        """
330        returns the pahse the well works on        returns the pahse the well works on
331        """        """
332        return self.WATER        return self.__phase
333          
334  class PeacemanWell(Well):     def isOpen(self, t):
335         """
336         returns True is the well is open at time t
337         """
338         out = False
339         t0=min(t, self.__schedule[0])
340         i=0
341         while t > t0:
342           if out:
343         out=False
344           else:
345         out=True  
346           i+=1
347           if i < len(self.__schedule):
348          t0=self.__schedule[i]
349           else:
350          t0=t
351         return out
352          
353       def setWaterProduction(self, v):
354          self.water_production=v
355       def getWaterProduction(self):
356          return self.water_production
357       def setGasProduction(self, v):
358          self.gas_production=v
359       def getGasProduction(self):
360          return self.gas_production
361          
362    class VerticalPeacemanWell(Well):
363     """     """
364     defines a well using Peaceman's formula     defines a well using Peaceman's formula
365     """     """
366     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],
367                 perm=[1.*U.cPoise,1.*U.cPoise, 1.*U.cPoise], phase=Well.WATER, s=0):
368           """
369           set up well
370          
371           :param BHP: ottom-hole pressure
372           :param Q: flow rate
373           :param r: well radius
374           :param X: location
375           :param D: dimension of well block
376           :param perm: permeabilities
377           :param s: skin factor
378           """
379           Well.__init__(self, name, domain, Q=Q, schedule=schedule, phase=phase,BHP_limit=BHP_limit)
380           r_el=0.28 * sqrt( sqrt(perm[1]/perm[0]) * D[0]**2 +  sqrt(perm[0]/perm[1]) * D[1]**2 )\
381                             / ( (perm[1]/perm[0])**0.25 + (perm[1]/perm[0])**0.25 )
382           self.__PI = 2 * math.pi * D[2] * sqrt(perm[1]*perm[0]) / (log(r_el/r) + s)
383           self.__r = r
384    
385           self.__D=D
386           self.r_el=r_el
387           self.X0=X0[:self.domain.getDim()]
388          
389           x=Function(domain).getX()
390           self.chi = 1.
391           for i in range(domain.getDim()):
392            self.chi = self.chi * whereNegative(abs(x[i]-X0[i])-D[i]/2)
393    
394           self.chi*=1./(D[0]*D[1]*D[2])
395          
396          
397           #self.chi=0.00000001
398       def getLocation(self):
399           return self.X0
400          
401       def getGeometry(self):
402           """
403           returns a function representing the geometry of the well.
404           typically a Gaussian profile around the location of the well.
405          
406           :note: needs to be overwritten
407           """
408           return self.chi
409    
410       def getProductivityIndex(self):
411           """
412           returns the productivity index of the well.
413           """
414           return self.__PI
415        
416    
417  class DualPorosity(object):  class DualPorosity(object):
418     """     """
# Line 188  class DualPorosity(object): Line 436  class DualPorosity(object):
436        :type phi_m: `MaterialPropertyWithDifferential` or None        :type phi_m: `MaterialPropertyWithDifferential` or None
437        :param phi: total porosity if None phi=1.        :param phi: total porosity if None phi=1.
438        :type phi: scalar or None        :type phi: scalar or None
439        :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
440        :type S_fg: `MaterialPropertyWithDifferential`        :type S_fg: `MaterialPropertyWithDifferential`
441        :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
442        :type S_mg: `MaterialPropertyWithDifferential`        :type S_mg: `MaterialPropertyWithDifferential`
# Line 235  class DualPorosity(object): Line 483  class DualPorosity(object):
483        self.rho_w = rho_w        self.rho_w = rho_w
484        self.rho_g = rho_g            self.rho_g = rho_g    
485        self.wells=wells        self.wells=wells
486          self.t =0
487          
488          self.__iter_max=1
489          self.__rtol=1.e-4
490          self.verbose=False
491          self.XI=0.5
492       def setIterationControl(self, iter_max=None, rtol=None, verbose=None, xi=None):
493         """
494         sets parameters to control iteration process
495         """
496         if iter_max !=None: self.__iter_max=iter_max
497         if rtol !=None: self.__rtol = rtol
498         if verbose !=None: self.verbose=verbose
499         if xi !=None: self.XI=xi
500        
501       def update(self, dt):
502             self.u, self.u_old = self.u.copy(), self.u
503             n=0
504             rerr=1.
505             while n < self.__iter_max and rerr > self.__rtol:
506                u=self.solvePDE(dt)
507            norm_u=Lsup(u)
508            norm_e=Lsup(u-self.u)
509            
510            if norm_u > 0:
511                rerr=norm_e/norm_u
512            else:
513                rerr=norm_e
514            if self.verbose: print "iteration step %d: relative change = %e"%(n, rerr)
515            n+=1
516            self.u=u
517             print "iteration completed."
518             self.t+=dt
519    
520    
521  class PorosityOneHalfModel(DualPorosity):  class PorosityOneHalfModel(DualPorosity):
# Line 247  class PorosityOneHalfModel(DualPorosity) Line 525  class PorosityOneHalfModel(DualPorosity)
525        This corresponds to the coal bed methan model in the ECLIPSE code.        This corresponds to the coal bed methan model in the ECLIPSE code.
526        """        """
527    
528        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,
529               perm_f_0=None, perm_f_1=None, perm_f_2=None,               perm_f_0=None, perm_f_1=None, perm_f_2=None,
530               k_w =None, k_g=None, mu_w =None, mu_g =None,               k_w =None, k_g=None, mu_w =None, mu_g =None,
531               rho_w =None, rho_g=None,               rho_w =None, rho_g=None, sigma=0, A_mg=0, f_rg=1.,  
532                 wells=[]):                 wells=[]):
533       """       """
534       set up       set up
# Line 263  class PorosityOneHalfModel(DualPorosity) Line 541  class PorosityOneHalfModel(DualPorosity)
541       :type phi: scalar or None       :type phi: scalar or None
542       :param L_g: gas adsorption as function of gas pressure       :param L_g: gas adsorption as function of gas pressure
543       :type L_g: `MaterialPropertyWithDifferential`       :type L_g: `MaterialPropertyWithDifferential`
544       :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
545       :type S_fg: `MaterialPropertyWithDifferential`       :type S_fg: `MaterialPropertyWithDifferential`
546       :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
547       :type perm_f_0: scalar       :type perm_f_0: scalar
# Line 285  class PorosityOneHalfModel(DualPorosity) Line 563  class PorosityOneHalfModel(DualPorosity)
563       :type rho_g: `MaterialPropertyWithDifferential`       :type rho_g: `MaterialPropertyWithDifferential`
564       :param wells : list of wells       :param wells : list of wells
565       :type wells: list of `Well`       :type wells: list of `Well`
566         :param sigma: shape factor for gas matrix diffusion
567         :param A_mg: diffusion constant for gas adsorption
568         :param f_rg: gas re-adsorption factor
569       """       """
570        
571       DualPorosity.__init__(self, domain,       DualPorosity.__init__(self, domain,
572                    phi_f=phi_f, phi_m=None, phi=phi,                    phi_f=phi_f, phi_m=None, phi=phi,
573                    S_fg=S_fg, S_mg=None,                    S_fg=None, S_mg=None,
574                    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,
575                    perm_m_0=None , perm_m_1=None, perm_m_2=None,                    perm_m_0=None , perm_m_1=None, perm_m_2=None,
576                    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,
577                    rho_w =rho_w, rho_g=rho_g,                    rho_w =rho_w, rho_g=rho_g,
578                    wells=wells)                    wells=wells)
579       self.L_g=L_g       self.L_g=L_g
580         self.sigma = sigma
581         self.A_mg = A_mg
582         self.f_rg  = f_rg
583       self.__pde=LinearPDE(self.domain, numEquations=3, numSolutions =3)       self.__pde=LinearPDE(self.domain, numEquations=3, numSolutions =3)
584            
585            
# Line 305  class PorosityOneHalfModel(DualPorosity) Line 589  class PorosityOneHalfModel(DualPorosity)
589       """       """
590       return self.__pde.getSolverOptions()       return self.__pde.getSolverOptions()
591            
592        def solvePDE(self):        def getState(self):
593             return self.u[0], self.u[1],  self.u[2]
      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]  
594    
595          def getOldState(self):
596         return self.u_old[0], self.u_old[1],  self.u_old[2]
597    
598       phi_f   =self.phi_f.getValue(p_fg)        def setInitialState(self, p=1.*U.atm, S_fg=0,  C_mg=None):
599       dphi_fdp=self.phi_f.getValueDifferential(p_fg)          """
600            sets the initial state
601            
602            :param p: pressure
603            :param S_fg: gas saturation in fractured rock
604            :param C_mg: gas concentration in coal matrix. if not given it is calculated
605                using the  gas adsorption curve.
606            """    
607            self.u=self.__pde.createSolution()
608            self.u[0]=p
609            self.u[1]=S_fg
610            if C_mg == None:
611              self.u[2]= self.L_g(p)*self.rho_g.getFormationFactor(p)
612            else:
613              self.u[2]=C_mg
614          
615          def solvePDE(self, dt):
616            
617       S_fg=  self.S_fg.getValue(p_fc)       C_couple=1
      dS_fgdp=   self.S_fg.getValueDifferential(p_fc)  
      S_fw=  1-S_fg  
618            
619       rho_fw     = self.rho_w.getValue(p_fw)       p_f, S_fg, C_mg=self.getState()
620       drho_fwdp  = self.rho_w.getValueDifferential(p_fw)       p_f_old, S_fg_old, C_mg_old=self.getOldState()
621       rho_fg = self.rho_g.getValue(p_fg)  
622       drho_fgdp  = self.rho_g.getValueDifferential(p_fg)           S_fw=1-S_fg
623    
624         if self.verbose:
625              print "p_f range = ",inf(p_f),sup(p_f)
626              print "S_fg range = ",inf(S_fg),sup(S_fg)
627              print "S_fw range = ",inf(S_fw),sup(S_fw)
628              print "C_mg range = ",inf(C_mg),sup(C_mg)
629    
630             k_fw=self.k_w(S_fw)
631             if self.verbose: print "k_fw range = ",inf(k_fw),sup(k_fw)
632    
633             k_fg=self.k_g(S_fg)
634             if self.verbose: print "k_fg range = ",inf(k_fg),sup(k_fg)
635    
636         mu_fw=self.mu_w(p_f)
637             if self.verbose: print "mu_fw range = ",inf(mu_fw),sup(mu_fw)
638    
639         mu_fg=self.mu_g(p_f)
640             if self.verbose: print "mu_fg range = ",inf(mu_fg),sup(mu_fg)
641            
642    
643         phi_f   =self.phi_f.getValue(p_f)
644         dphi_fdp=self.phi_f.getValueDifferential(p_f)
645             if self.verbose: print "phi_f range = ",inf(phi_f),sup(phi_f)," (slope %e,%e)"%(inf(dphi_fdp), sup(dphi_fdp))
646            
647       L_g       = self.getValue(p_fg)       rho_fw     = self.rho_w.getValue(p_f)
648       dL_gdp_fg =  self.rho_w.getValueDifferential(p_fg)       drho_fwdp  = self.rho_w.getValueDifferential(p_f)
649             if self.verbose: print "rho_fw range = ",inf(rho_fw),sup(rho_fw)," (slope %e,%e)"%(inf(drho_fwdp), sup(drho_fwdp))
650    
651             rho_fg = self.rho_g.getValue(p_f)
652         drho_fgdp  = self.rho_g.getValueDifferential(p_f)
653             if self.verbose: print "rho_fg range = ",inf(rho_fg),sup(rho_fg)," (slope %e,%e)"%(inf(drho_fgdp), sup(drho_fgdp))
654            
655       A_fw = rho_fw * k_fw/mu_fw       L_g_0       = self.L_g(p_f)
656       A_fg = rho_fg * k_fg/mu_fg       dL_g_0dp    = self.L_g.getValueDifferential(p_f)
657         FF_g=self.rho_g.getFormationFactor(p_f)
658         dFF_gdp=self.rho_g.getFormationFactorDifferential(p_f)
659            
660       m = whereNegative(L_g - C_mg)         L_g = L_g_0 * FF_g
661       B = self.sigma * self.a_mg * (m + (1-m) * self.f_rg * S_fg )       dL_gdp =  0 * (dL_g_0dp * FF_g + L_g_0 * dFF_gdp)
662    
663             if self.verbose:
664              print "L_g_0 range = ",inf(L_g_0),sup(L_g_0)
665              print "L_g range = ",inf(L_g),sup(L_g)," (slope %e,%e)"%(inf(dL_gdp), sup(dL_gdp))
666            
667            
668       E_fww= phi_f * ( rho_fw * dS_fgdp + S_fw * drho_fwdp )       A_fw = rho_fw * k_fw/mu_fw
669       E_fwg= rho_fw * ( S_fw * dphi_fdp - phi_f * dS_fgdp )       A_fg = rho_fg * k_fg/mu_fg
      E_fgw= - phi_f * rho_fg * dS_fgdp  
      E_fgg= S_fg * rho_fg  * dphi_fdp + phi_f * rho_fg * dS_fgdp + phi_f  * S_fg * drho_fgdp  
   
   
670            
671       Q_w=0.       m = whereNegative(L_g - C_mg)
672       Q_g=0.       B = self.sigma * self.A_mg * (m + (1-m) * self.f_rg * S_fg )
673       R_w=0.       print "XXXX B =",B
674       R_g=0.       print "XXXX self.sigma  =",self.sigma
675         print "XXXX self.A_mg  =",self.A_mg
676        
677         E_fpp= S_fw * (  rho_fw * dphi_fdp + phi_f  * drho_fwdp )
678         E_fps=  -  phi_f * rho_fw
679         E_fsp= S_fg *( rho_fg * dphi_fdp - phi_f *  drho_fgdp )
680         E_fss= phi_f * rho_fg
681        
682         F_fw=0.
683         F_fg=0.
684         D_fw=0.
685         D_fg=0.
686            
687       for I in self.wells:       for I in self.wells:
688          chi_I= I.getGeometry()            chi_I= I.getGeometry()
689          Pi_I = I.getProductivityIndex()            loc=Locator(Function(self.domain),I.getLocation())
690          P_I = I.getBHP()            Pi_I = I.getProductivityIndex()
691          Q_I = I.getFlowRate()            A_fw_I= loc(A_fw)
692                      A_fg_I= loc(A_fg)
693          if self.getPhase == Well.WATER:            BHP_limit_I=I.getBHPLimit()
694             Q_w = Q_w + Q_I        * chi_I            
695             Q_g = Q_g + Pi_I * P_I * chi_I            if I.isOpen(self.t+dt):
696             R_g = R_g + Pi_I       * chi_I          if self.verbose: print "well %s is open."%I.name
697          else:          if I.getPhase() == Well.WATER:
698             Q_g = Q_g + Q_I        * chi_I              q_I = self.rho_w.rho_0 * I.getFlowRate()
699             Q_w = Q_w + Pi_I * P_I * chi_I              p_I_ref=loc(p_f)-q_I/(A_fw_I * Pi_I)
700             R_w = R_w + Pi_I       * chi_I          else:
701                          q_I = self.rho_g.rho_0 * I.getFlowRate()
702       F_fw_Y = A_fw * Q_w              p_I_ref=loc(p_f)-q_I/(A_fg_I * Pi_I)
703       F_fg_Y = A_fg * Q_g              
704       D_fgg = - A_fg * R_g          print "ZZZ =",loc(p_f), q_I, self.rho_w.rho_0, I.getFlowRate(), A_fw_I, Pi_I, q_I/(A_fw_I * Pi_I)
705       D_fww = - A_fw * R_w  
706            if BHP_limit_I > p_I_ref:
707              BHP_I=BHP_limit_I
708              D_fw = D_fw + Pi_I * A_fw_I *              chi_I
709              F_fw = F_fw - Pi_I * A_fw_I * BHP_limit_I *chi_I
710              D_fg = D_fg + Pi_I * A_fg_I *              chi_I
711              F_fg = F_fg - Pi_I * A_fg_I * BHP_limit_I *chi_I
712            else:
713              if I.getPhase() == Well.WATER:
714                  F_fw = F_fw +  q_I  * chi_I
715                  F_fg = F_fg +  A_fg_I/A_fw_I *  q_I *chi_I
716              else:
717                  F_fg = F_fg +  q_I  * chi_I
718                  F_fw = F_fw +  A_fw_I/A_fg_I *  q_I *chi_I
719              BHP_I=p_I_ref
720              else:
721              if self.verbose: print "well %s is closed."%I.name
722              BHP_I=loc(p_f)
723              
724              if self.verbose: print "well %s:BHP = %e"%(I.name, BHP_I)
725              I.setBHP(BHP_I)
726    
727         F_fp_Y = - F_fw
728         F_fs_Y = - F_fg
729         D_fpp =  D_fw
730         D_fsp =  D_fg
731    
732            
733       if self.domain.getDim() == 3:       if self.domain.getDim() == 3:
734          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]
735          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]
736       else:       else:
737          F_fw_X = 0 * kronecker(self.domain)[1]          F_fp_X = 0 * kronecker(self.domain)[1]
738          F_fg_X = 0 * kronecker(self.domain)[1]          F_fs_X = 0 * kronecker(self.domain)[1]
739                    
740       F_mg = B * ( L_g - dL_gdp_fg * p_fg )       #F_mg_Y = B * ( L_g - dL_gdp * p_f )
741             F_mg_Y = B * L_g
742    
743        
744       D=self.__pde.getNewCoefficient("D")       D=self.__pde.createCoefficient("D")
745       A=self.__pde.getNewCoefficient("A")       A=self.__pde.createCoefficient("A")
746       Y=self.__pde.getNewCoefficient("Y")       Y=self.__pde.createCoefficient("Y")
747       X=self.__pde.getNewCoefficient("X")       X=self.__pde.createCoefficient("X")
748            
749       dtXI = dt*self.XI       dtXI = dt*self.XI
750       dtcXI = dt*(1-self.XI)       dtcXI = dt*(1-self.XI)
751        
752       D[0,0]=E_fww + dtXI * D_fww       D[0,0]=E_fpp + dtXI * D_fpp
753       D[0,1]=E_fwg       D[0,1]=E_fps
754       D[1,0]=E_fgw       D[1,0]=E_fsp + dtXI * D_fsp
755       D[1,1]=E_fgg + dtXI * D_fgg       D[1,1]=E_fss
756       D[1,2]=rho_fg       D[1,2]=rho_fg * C_couple
757       D[2,1]= dtXI * B * dL_gdp_fg       #D[2,0]= - dtXI * B * dL_gdp
758       D[2,2]= 1 + dtXI * B       D[2,2]= 1 + dtXI * B
759            
      H_fw = dtcXI * A_fw * grad(p_fw_old)  
      H_fg = dtcXI * A_fg * grad(p_fg_old)  
760            
761       for i in range(self.domain.getDim()):       for i in range(self.domain.getDim()):
762          A[0,i,0,i] = dtXI * A_fw * self.perm_f[i]          A[0,i,0,i] = dtXI * A_fw * self.perm_f[i]
763          A[1,i,1,i] = dtXI * A_fg * self.perm_f[i]          A[1,i,1,i] = dtXI * A_fg * self.perm_f[i]
764            
765          X[0,i]= dt * F_fw_X - H_fw * self.perm_f[i]       g= grad(p_f_old) * dtcXI * self.perm_f[0:self.domain.getDim()]
766          X[1,i]= dt * F_fg_X - H_fg * self.perm_f[i]           X[0,:]=  - A_fw * g  + dt * F_fp_X
767                 X[1,:]=  - A_fg * g  + dt * F_fs_X
768       Y[0] = E_fww * p_fw_old + E_fwg * p_fg_old +                     dt * F_fw_Y - dtcXI * D_fww * p_fw_old  
769       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[0] = E_fpp *  p_f_old + E_fps * S_fg_old +                                dt * F_fp_Y - dtcXI * D_fpp * p_f_old
770       Y[2] = C_mg_old                                                + dt * F_mg_Y - dtcXI * ( dL_gdp_fg * p_fg_old + C_mg_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
771         Y[2] = C_mg_old                                                           + dt * F_mg_Y - dtcXI * B * C_mg_old  # - dL_gdp * p_f_old)
772        
773         print "HHH Y[0] =", Y[0]
774         print "HHH A_fw = ",A_fw
775         print "HHH A_fg= ",A_fg
776         print "HHH F_fg = ",F_fg
777         print "HHH F_fw = ",F_fw
778         print "HHH perm_f = ",self.perm_f
779         print "HHH k = ",(A_fw*self.perm_f[0])/D[0,0]
780         print "HHH k = ",(A_fw*self.perm_f[1])/D[0,0]
781         print "HHH X[0,:]= ",X[0,:]
782         print "HHH D[0,0] = ",D[0,0]
783         print "HHH D[0,1] = ",D[0,1]
784         print "HHH D[0,2] = ",D[0,2]
785        
786         print "hhh B= ",B
787            
788    
789          
790            
791         #print "HHH Y[1] = ",Y[1]
792         #print "HHH X[1,:] = ",X[1,:]
793         #print "HHH D[1,0] = ",D[1,0]
794         #print "HHH D[1,1] = ",D[1,1]
795         #print "HHH D[1,2] = ",D[1,2]
796         #print "HHH A_fg =",  A_fg
797        
798    
799    
800    
801         #1/0
802       self.__pde.setValue(A=A, D=D, X=X, Y=Y)       self.__pde.setValue(A=A, D=D, X=X, Y=Y)
803            
804       u2 = self.__pde.getSolution()       u2 = self.__pde.getSolution()
805         # this is deals with round-off errors to maintain physical meaningful values
806         # we should do this in a better way to detect values that are totally wrong
807         u2[1]=clip(u2[1],minval=0, maxval=1)
808         u2[2]=clip(u2[2],minval=0)
809        
810         # update well production
811         p_f=u2[0]
812         for I in self.wells:
813              loc=Locator(Function(self.domain),I.getLocation())
814              Pi_I = I.getProductivityIndex()
815              q_I = I.getFlowRate()
816              A_fw_I= loc(A_fw)
817              A_fg_I= loc(A_fg)
818              p_f_I=loc(A_fg)
819              BHP_limit_I=I.getBHPLimit()
820              BHP=I.getBHP()
821              
822              
823              rho_fw_I = loc(rho_fw)
824              rho_fg_I = loc(rho_fg)
825              A_fg = rho_fg * k_fg/mu_fg
826        
827              if I.isOpen(self.t+dt):
828            if BHP_limit_I < BHP:
829              if I.getPhase() == Well.WATER:
830                  q_I_w =  q_I
831                  q_I_g =  A_fg_I/A_fw_I *  q_I
832              else:
833                  q_I_g =  q_I
834                  q_I_w =  A_fw_I/A_fg_I *  q_I
835            else:
836                q_I_w = Pi_I * A_fw_I/rho_fw_I * (p_f_I - BHP_I)
837                q_I_g = Pi_I * A_fg_I/rho_fg_I * (p_f_I - BHP_I)
838              else:
839              q_I_g =0
840              q_I_w =0
841              I.setWaterProduction(q_I_w * rho_fw_I/self.rho_w.rho_0)
842              I.setGasProduction(q_I_g * rho_fg_I/self.rho_g.rho_0 )
843       return u2       return u2

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

  ViewVC Help
Powered by ViewVC 1.1.26