/[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 3494 by gross, Tue Apr 5 04:12:08 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  class InterpolationTable(object):      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    class WaterViscosity(MaterialProperty):
129      """
130      defines viscosity of water
131      
132      mu=mu_ref /(1 + X + X*X/2)
133      
134      with X=-C*(p-p_ref)
135      """
136    
137      def __init__(self, mu_ref = 0.3*U.cPoise, p_ref = 1.*U.atm , C = 0./U.bar):
138        """
139        set up
140        """
141        self.mu_ref = mu_ref
142        self.p_ref = p_ref
143        self.C=C
144      def getValue(self, p):
145          """
146          returns the viscosity for given pressure p
147          """
148          X= -self.C * ( p - self.p_ref )
149          return self.mu_ref/(1+ X * (1+ X/2) )  
150          
151    class GasDensity(MaterialPropertyWithDifferential):
152        """
153        set water density as
154          
155              rho = gravity * rho_air_s /B(p)
156              
157        where B is given by an interpolation table
158        """
159        def __init__(self, p, B, gravity=1., rho_air_s=1.2041*U.kg/U.m**3):
160          self.rho_ref =rho_air_s * gravity
161          self.rho_0 =rho_air_s * gravity
162          self.tab=InterpolationTable(x=p, y=B)
163        
164        def getValue(self, p):
165          """
166          returns the density for given pressure p
167          """
168          return self.rho_ref/self.tab.getValue(p)
169          
170        def getValueDifferential(self,  p):
171          """
172          """
173          B    = self.tab.getValue(p)
174          dBdp = self.tab.getValueDifferential(p)
175          return  -self.rho_ref * dBdp/(B * B)
176    
177    class InterpolationTable(MaterialPropertyWithDifferential):
178     """     """
179     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
180     """     """
# Line 72  class InterpolationTable(object): Line 184  class InterpolationTable(object):
184        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
185        the value for x is set to y[0] for        the value for x is set to y[0] for
186        """        """
187          MaterialPropertyWithDifferential.__init__(self)
188        if len(x) != len(y):        if len(x) != len(y):
189       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."
190        if len(x) > 0:        if len(x) < 1 :
191       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."
192                
193        x_ref=x[0]        x_ref=x[0]
# Line 85  class InterpolationTable(object): Line 198  class InterpolationTable(object):
198        self.__x = x        self.__x = x
199        self.__y = y        self.__y = y
200        self.__obeyBounds = obeyBounds        self.__obeyBounds = obeyBounds
201          
202     def getValue(self, x):     def getValue(self, x):
203        """        """
204        returns the interpolated values of x        returns the interpolated values of x
205        """        """
206        x0=self.__x[0]        X=self.__x
207          Y=self.__y
208          
209          x0=X[0]
210        m0=whereNegative(x-x0)        m0=whereNegative(x-x0)
211        if self.__obeyBounds:        if self.__obeyBounds:
212       if sup(m0) > 0:       if sup(m0) > 0:
213          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])
214        out=self.__x[0]        out=self.__x[0]
215        for i in range(1,len(self.__x)):        for i in range(1,len(X)):
216        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]  
217        out = out * m0 + z * (1-m0)        out = out * m0 + z * (1-m0)
218        m0=whereNegative(x-x[i])        m0=whereNonPositive(x-X[i])
219              
220          if self.__obeyBounds:
221            if inf(m0) < 1 :
222               raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
223          else:
224            out = out * m0 + V[-1] * (1-m0)
225          return out
226    
227       def getValueDifferential(self, x):
228          X=self.__x
229          Y=self.__y
230    
231          x0=X[0]
232          m0=whereNegative(x-x0)
233          if self.__obeyBounds:
234         if sup(m0) > 0:
235            raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
236          out=0.
237          for i in range(1,len(X)):
238          z=(Y[i]-Y[i-1])/(X[i]-X[i-1])
239          out = out * m0 + z * (1-m0)
240          m0=whereNonPositive(x-X[i])
241            
242        if self.__obeyBounds:        if self.__obeyBounds:
243          if inf(m0) < 1:          if inf(m0) < 1:
244             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])
245        else:        else:
246          out = out * m0 + y[-1] * (1-m0)          out = out * m0
247        return out        return out  
248          
249    
250  class Well(object):  class Well(object):
251     """     """
# Line 117  class Well(object): Line 256  class Well(object):
256     """     """
257     WATER="Water"     WATER="Water"
258     GAS="Gas"     GAS="Gas"
259     def __init__(self,  *args, **kwargs):     def __init__(self, name, domain, Q=0., schedule=[0.], phase="Water", BHP_limit=1.*U.atm, *args, **kwargs):
260         """         """
261         set up well         set up well
262         """         """
263         pass         self.__schedule=schedule
264           self.__phase=phase
265           self.__Q = Q
266           self.__BHP_limit=BHP_limit
267           self.name=name
268           self.domain=domain
269          
270     def getGeometry(self):     def getGeometry(self):
271         """         """
272         returns a function representing the geometry of the well.         returns a function representing the geometry of the well.
# Line 144  class Well(object): Line 289  class Well(object):
289        """        """
290        returns the flow rate        returns the flow rate
291        """        """
292        return 0.        return self.__Q
293        
294       def getBHPLimit(self):
295          """
296          return bottom-hole pressure limit
297          
298          :note: needs to be overwritten
299          """
300          return self.__BHP_limit
301                
302     def getBHP(self):     def getBHP(self):
303        """        """
# Line 152  class Well(object): Line 305  class Well(object):
305                
306        :note: needs to be overwritten        :note: needs to be overwritten
307        """        """
308        raise NotImplementedError        return self.__BHP
309      
310       def setBHP(self, BHP):
311          """
312          sets bottom-hole pressure
313          """
314          self.__BHP= BHP
315                
316     def getPhase(self):     def getPhase(self):
317        """        """
318        returns the pahse the well works on        returns the pahse the well works on
319        """        """
320        return self.WATER        return self.__phase
321          
322       def isOpen(self, t):
323         """
324         returns True is the well is open at time t
325         """
326         out = False
327         t0=min(t, self.__schedule[0])
328         i=0
329         while t > t0:
330           if out:
331         out=False
332           else:
333         out=True  
334           i+=1
335           if i < len(self.__schedule):
336          t0=self.__schedule[i]
337           else:
338          t0=t
339         return out
340    
341  class PeacemanWell(Well):  class VerticalPeacemanWell(Well):
342     """     """
343     defines a well using Peaceman's formula     defines a well using Peaceman's formula
344     """     """
345     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],
346                 perm=[1.*U.cPoise,1.*U.cPoise, 1.*U.cPoise], phase=Well.WATER, s=0):
347           """
348           set up well
349          
350           :param BHP: ottom-hole pressure
351           :param Q: flow rate
352           :param r: well radius
353           :param X: location
354           :param D: dimension of well block
355           :param perm: permeabilities
356           :param s: skin factor
357           """
358           Well.__init__(self, name, domain, Q=Q, schedule=schedule, phase=phase,BHP_limit=BHP_limit)
359           r_el=0.28 * sqrt( sqrt(perm[1]/perm[0]) * D[0]**2 +  sqrt(perm[0]/perm[1]) * D[1]**2 )\
360                             / ( (perm[1]/perm[0])**0.25 + (perm[1]/perm[0])**0.25 )
361           self.__PI = 2 * math.pi * D[2] * sqrt(perm[1]*perm[0]) / (log(r_el/r) + s)
362           self.__r = r
363    
364           self.__D=D
365           self.r_el=r_el
366           self.X0=X0[:self.domain.getDim()]
367          
368           x=Function(domain).getX()
369           self.chi = 1.
370           for i in range(domain.getDim()):
371            self.chi = self.chi * whereNegative(abs(x[i]-X0[i])-D[i]/2)
372    
373           self.chi*=1./(D[0]*D[1]*D[2])
374          
375          
376           #self.chi=0.00000001
377       def getLocation(self):
378           return self.X0
379          
380       def getGeometry(self):
381           """
382           returns a function representing the geometry of the well.
383           typically a Gaussian profile around the location of the well.
384          
385           :note: needs to be overwritten
386           """
387           return self.chi
388    
389       def getProductivityIndex(self):
390           """
391           returns the productivity index of the well.
392           """
393           return self.__PI
394        
395    
396  class DualPorosity(object):  class DualPorosity(object):
397     """     """
# Line 188  class DualPorosity(object): Line 415  class DualPorosity(object):
415        :type phi_m: `MaterialPropertyWithDifferential` or None        :type phi_m: `MaterialPropertyWithDifferential` or None
416        :param phi: total porosity if None phi=1.        :param phi: total porosity if None phi=1.
417        :type phi: scalar or None        :type phi: scalar or None
418        :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
419        :type S_fg: `MaterialPropertyWithDifferential`        :type S_fg: `MaterialPropertyWithDifferential`
420        :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
421        :type S_mg: `MaterialPropertyWithDifferential`        :type S_mg: `MaterialPropertyWithDifferential`
# Line 235  class DualPorosity(object): Line 462  class DualPorosity(object):
462        self.rho_w = rho_w        self.rho_w = rho_w
463        self.rho_g = rho_g            self.rho_g = rho_g    
464        self.wells=wells        self.wells=wells
465          self.t =0
466          
467          self.__iter_max=1
468          self.__rtol=1.e-4
469          self.verbose=False
470          self.XI=0.5
471       def setIterationControl(self, iter_max=None, rtol=None, verbose=None, xi=None):
472         """
473         sets parameters to control iteration process
474         """
475         if iter_max !=None: self.__iter_max=iter_max
476         if rtol !=None: self.__rtol = rtol
477         if verbose !=None: self.verbose=verbose
478         if xi !=None: self.XI=xi
479        
480       def update(self, dt):
481             self.u, self.u_old = self.u.copy(), self.u
482             n=0
483             rerr=1.
484             while n < self.__iter_max and rerr > self.__rtol:
485                u=self.solvePDE(dt)
486            norm_u=Lsup(u)
487            norm_e=Lsup(u-self.u)
488            
489            if norm_u > 0:
490                rerr=norm_e/norm_u
491            else:
492                rerr=norm_e
493            if self.verbose: print "iteration step %d: relative change = %e"%(n, rerr)
494            n+=1
495            self.u=u
496             print "iteration completed."
497             self.t+=dt
498    
499    
500  class PorosityOneHalfModel(DualPorosity):  class PorosityOneHalfModel(DualPorosity):
# Line 247  class PorosityOneHalfModel(DualPorosity) Line 504  class PorosityOneHalfModel(DualPorosity)
504        This corresponds to the coal bed methan model in the ECLIPSE code.        This corresponds to the coal bed methan model in the ECLIPSE code.
505        """        """
506    
507        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,
508               perm_f_0=None, perm_f_1=None, perm_f_2=None,               perm_f_0=None, perm_f_1=None, perm_f_2=None,
509               k_w =None, k_g=None, mu_w =None, mu_g =None,               k_w =None, k_g=None, mu_w =None, mu_g =None,
510               rho_w =None, rho_g=None,               rho_w =None, rho_g=None, sigma=0, A_mg=0, f_rg=1.,  
511                 wells=[]):                 wells=[]):
512       """       """
513       set up       set up
# Line 263  class PorosityOneHalfModel(DualPorosity) Line 520  class PorosityOneHalfModel(DualPorosity)
520       :type phi: scalar or None       :type phi: scalar or None
521       :param L_g: gas adsorption as function of gas pressure       :param L_g: gas adsorption as function of gas pressure
522       :type L_g: `MaterialPropertyWithDifferential`       :type L_g: `MaterialPropertyWithDifferential`
523       :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
524       :type S_fg: `MaterialPropertyWithDifferential`       :type S_fg: `MaterialPropertyWithDifferential`
525       :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
526       :type perm_f_0: scalar       :type perm_f_0: scalar
# Line 285  class PorosityOneHalfModel(DualPorosity) Line 542  class PorosityOneHalfModel(DualPorosity)
542       :type rho_g: `MaterialPropertyWithDifferential`       :type rho_g: `MaterialPropertyWithDifferential`
543       :param wells : list of wells       :param wells : list of wells
544       :type wells: list of `Well`       :type wells: list of `Well`
545         :param sigma: shape factor for gas matrix diffusion
546         :param A_mg: diffusion constant for gas adsorption
547         :param f_rg: gas re-adsorption factor
548       """       """
549        
550       DualPorosity.__init__(self, domain,       DualPorosity.__init__(self, domain,
551                    phi_f=phi_f, phi_m=None, phi=phi,                    phi_f=phi_f, phi_m=None, phi=phi,
552                    S_fg=S_fg, S_mg=None,                    S_fg=None, S_mg=None,
553                    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,
554                    perm_m_0=None , perm_m_1=None, perm_m_2=None,                    perm_m_0=None , perm_m_1=None, perm_m_2=None,
555                    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,
556                    rho_w =rho_w, rho_g=rho_g,                    rho_w =rho_w, rho_g=rho_g,
557                    wells=wells)                    wells=wells)
558       self.L_g=L_g       self.L_g=L_g
559         self.sigma = sigma
560         self.A_mg = A_mg
561         self.f_rg  = f_rg
562       self.__pde=LinearPDE(self.domain, numEquations=3, numSolutions =3)       self.__pde=LinearPDE(self.domain, numEquations=3, numSolutions =3)
563            
564            
# Line 305  class PorosityOneHalfModel(DualPorosity) Line 568  class PorosityOneHalfModel(DualPorosity)
568       """       """
569       return self.__pde.getSolverOptions()       return self.__pde.getSolverOptions()
570            
571        def solvePDE(self):        def getState(self):
572         return self.u[0], self.u[1],  self.u[2]
573    
574          def getOldState(self):
575         return self.u_old[0], self.u_old[1],  self.u_old[2]
576    
577          def setInitialState(self, p=1.*U.atm, S_fg=0,  C_mg=None):
578            """
579            sets the initial state
580            
581            :param p: pressure
582            :param S_fg: gas saturation in fractured rock
583            :param C_mg: gas concentration in coal matrix. if not given it is calculated
584                using the  gas adsorption curve.
585            """    
586            self.u=self.__pde.createSolution()
587            self.u[0]=p
588            self.u[1]=S_fg
589            if C_mg == None:
590              self.u[2]= self.L_g(p)
591            else:
592              self.u[2]=C_mg
593          
594          def solvePDE(self, dt):
595            
596       p_fw=self.u[0]       C_couple=0
597       p_fg=self.u[1]      
598       p_fc=p_fg-p_fw       p_f, S_fg, C_mg=self.getState()
599       C_mg=self.u[3]       p_f_old, S_fg_old, C_mg_old=self.getOldState()
600          
601                   S_fw=1-S_fg
602       p_fw_old=self.u_old[0]  
603       p_fg_old=self.u_old[1]       if self.verbose:
604       C_mg_old=self.u_old[3]            print "p_f range = ",inf(p_f),sup(p_f)
605              print "S_fg range = ",inf(S_fg),sup(S_fg)
606              print "S_fw range = ",inf(S_fw),sup(S_fw)
607              print "C_mg range = ",inf(C_mg),sup(C_mg)
608    
609             k_fw=self.k_w(S_fw)
610             if self.verbose: print "k_fw range = ",inf(k_fw),sup(k_fw)
611    
612             k_fg=self.k_g(S_fg)
613             if self.verbose: print "k_fg range = ",inf(k_fg),sup(k_fg)
614    
615       phi_f   =self.phi_f.getValue(p_fg)       mu_fw=self.mu_w(p_f)
616       dphi_fdp=self.phi_f.getValueDifferential(p_fg)           if self.verbose: print "mu_fw range = ",inf(mu_fw),sup(mu_fw)
617    
618         mu_fg=self.mu_g(p_f)
619             if self.verbose: print "mu_fg range = ",inf(mu_fg),sup(mu_fg)
620            
621    
622         phi_f   =self.phi_f.getValue(p_f)
623         dphi_fdp=self.phi_f.getValueDifferential(p_f)
624             if self.verbose: print "phi_f range = ",inf(phi_f),sup(phi_f)," (slope %e,%e)"%(inf(dphi_fdp), sup(dphi_fdp))
625            
626       S_fg=  self.S_fg.getValue(p_fc)       rho_fw     = self.rho_w.getValue(p_f)
627       dS_fgdp=   self.S_fg.getValueDifferential(p_fc)       drho_fwdp  = self.rho_w.getValueDifferential(p_f)
628       S_fw=  1-S_fg           if self.verbose: print "rho_fw range = ",inf(rho_fw),sup(rho_fw)," (slope %e,%e)"%(inf(drho_fwdp), sup(drho_fwdp))
629    
630             rho_fg = self.rho_g.getValue(p_f)
631         drho_fgdp  = self.rho_g.getValueDifferential(p_f)
632             if self.verbose: print "rho_fg range = ",inf(rho_fg),sup(rho_fg)," (slope %e,%e)"%(inf(drho_fgdp), sup(drho_fgdp))
633            
634       rho_fw     = self.rho_w.getValue(p_fw)       L_g       = self.L_g.getValue(p_f)
635       drho_fwdp  = self.rho_w.getValueDifferential(p_fw)       dL_gdp =  self.rho_w.getValueDifferential(p_f)
636       rho_fg = self.rho_g.getValue(p_fg)           if self.verbose: print "L_g range = ",inf(L_g),sup(L_g)," (slope %e,%e)"%(inf(dL_gdp), sup(dL_gdp))
      drho_fgdp  = self.rho_g.getValueDifferential(p_fg)  
637            
      L_g       = self.getValue(p_fg)  
      dL_gdp_fg =  self.rho_w.getValueDifferential(p_fg)  
638            
639       A_fw = rho_fw * k_fw/mu_fw       A_fw = rho_fw * k_fw/mu_fw
640       A_fg = rho_fg * k_fg/mu_fg       A_fg = rho_fg * k_fg/mu_fg
641            
642       m = whereNegative(L_g - C_mg)         m = whereNegative(L_g - C_mg)  
643       B = self.sigma * self.a_mg * (m + (1-m) * self.f_rg * S_fg )       B = self.sigma * self.A_mg * (m + (1-m) * self.f_rg * S_fg )
644            
645            
646       E_fww= phi_f * ( rho_fw * dS_fgdp + S_fw * drho_fwdp )       E_fpp= S_fw * (  rho_fw * dphi_fdp + phi_f  * drho_fwdp )
647       E_fwg= rho_fw * ( S_fw * dphi_fdp - phi_f * dS_fgdp )       E_fps=  -  phi_f * rho_fw
648       E_fgw= - phi_f * rho_fg * dS_fgdp       E_fsp= S_fg *( rho_fg * dphi_fdp - phi_f *  drho_fgdp )
649       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
650        
651         F_fw=0.
652             F_fg=0.
653       Q_w=0.       D_fw=0.
654       Q_g=0.       D_fg=0.
      R_w=0.  
      R_g=0.  
655            
656       for I in self.wells:       for I in self.wells:
657          chi_I= I.getGeometry()  
658          Pi_I = I.getProductivityIndex()            chi_I= I.getGeometry()
659          P_I = I.getBHP()            loc=Locator(Function(self.domain),I.getLocation())
660          Q_I = I.getFlowRate()            Pi_I = I.getProductivityIndex()
661                      A_fw_I= loc(A_fw)
662          if self.getPhase == Well.WATER:            A_fg_I= loc(A_fg)
663             Q_w = Q_w + Q_I        * chi_I            BHP_limit_I=I.getBHPLimit()
664             Q_g = Q_g + Pi_I * P_I * chi_I            
665             R_g = R_g + Pi_I       * chi_I            if I.isOpen(self.t+dt):
666          else:          if self.verbose: print "well %s is open."%I.name
667             Q_g = Q_g + Q_I        * chi_I          if I.getPhase() == Well.WATER:
668             Q_w = Q_w + Pi_I * P_I * chi_I              q_I = self.rho_w.rho_0 * I.getFlowRate()
669             R_w = R_w + Pi_I       * chi_I              p_I_ref=loc(p_f)-q_I/(A_fw_I * Pi_I)
670                      else:
671       F_fw_Y = A_fw * Q_w              q_I = self.rho_g.rho_0 * I.getFlowRate()
672       F_fg_Y = A_fg * Q_g              p_I_ref=loc(p_f)-q_I/(A_fg_I * Pi_I)
673       D_fgg = - A_fg * R_g              
674       D_fww = - A_fw * R_w          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)
675            1/0
676            if BHP_limit_I > p_I_ref:
677              BHP_I=BHP_limit_I
678              D_fw = D_fw + Pi_I * A_fw_I *              chi_I
679              F_fw = F_fw - Pi_I * A_fw_I * BHP_limit_I *chi_I
680              D_fg = D_fg + Pi_I * A_fg_I *              chi_I
681              F_fg = F_fg - Pi_I * A_fg_I * BHP_limit_I *chi_I
682            else:
683              if I.getPhase() == Well.WATER:
684                  F_fw = F_fw +  q_I  * chi_I
685                  F_fg = F_fg +  A_fg_I/A_fw_I *  q_I *chi_I
686              else:
687                  F_fg = F_fg +  q_I  * chi_I
688                  F_fw = F_fw +  A_fw_I/A_fg_I *  q_I *chi_I
689              BHP_I=p_I_ref
690              else:
691              if self.verbose: print "well %s is closed."%I.name
692              BHP_I=loc(p_f)
693              
694              if self.verbose: print "well %s:BHP = %e"%(I.name, BHP_I)
695              I.setBHP(BHP_I)
696    
697         F_fp_Y = - F_fw
698         F_fs_Y = - F_fg
699         D_fpp =  D_fw
700         D_fsp =  D_fg
701    
702            
703       if self.domain.getDim() == 3:       if self.domain.getDim() == 3:
704          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]
705          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]
706       else:       else:
707          F_fw_X = 0 * kronecker(self.domain)[1]          F_fp_X = 0 * kronecker(self.domain)[1]
708          F_fg_X = 0 * kronecker(self.domain)[1]          F_fs_X = 0 * kronecker(self.domain)[1]
709                    
710       F_mg = B * ( L_g - dL_gdp_fg * p_fg )       F_mg_Y = B * ( L_g - dL_gdp * p_f )
711    
712        
713       D=self.__pde.getNewCoefficient("D")  
714       A=self.__pde.getNewCoefficient("A")       D=self.__pde.createCoefficient("D")
715       Y=self.__pde.getNewCoefficient("Y")       A=self.__pde.createCoefficient("A")
716       X=self.__pde.getNewCoefficient("X")       Y=self.__pde.createCoefficient("Y")
717         X=self.__pde.createCoefficient("X")
718            
719       dtXI = dt*self.XI       dtXI = dt*self.XI
720       dtcXI = dt*(1-self.XI)       dtcXI = dt*(1-self.XI)
721        
722       D[0,0]=E_fww + dtXI * D_fww       D[0,0]=E_fpp + dtXI * D_fpp
723       D[0,1]=E_fwg       D[0,1]=E_fps
724       D[1,0]=E_fgw       D[1,0]=E_fsp + dtXI * D_fsp
725       D[1,1]=E_fgg + dtXI * D_fgg       D[1,1]=E_fss
726       D[1,2]=rho_fg       D[1,2]=rho_fg * C_couple
727       D[2,1]= dtXI * B * dL_gdp_fg       D[2,0]= - dtXI * B * dL_gdp
728       D[2,2]= 1 + dtXI * B       D[2,2]= 1 + dtXI * B
729            
      H_fw = dtcXI * A_fw * grad(p_fw_old)  
      H_fg = dtcXI * A_fg * grad(p_fg_old)  
730            
731       for i in range(self.domain.getDim()):       for i in range(self.domain.getDim()):
732          A[0,i,0,i] = dtXI * A_fw * self.perm_f[i]          A[0,i,0,i] = dtXI * A_fw * self.perm_f[i]
733          A[1,i,1,i] = dtXI * A_fg * self.perm_f[i]          A[1,i,1,i] = dtXI * A_fg * self.perm_f[i]
734            
735          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()]
736          X[1,i]= dt * F_fg_X - H_fg * self.perm_f[i]           X[0,:]=  - A_fw * g  + dt * F_fp_X
737                 X[1,:]=  - A_fg * g  + dt * F_fs_X
738       Y[0] = E_fww * p_fw_old + E_fwg * p_fg_old +                     dt * F_fw_Y - dtcXI * D_fww * p_fw_old  
739       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
740       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
741         Y[2] = C_mg_old                                                           + dt * F_mg_Y - dtcXI * B * ( C_mg_old - dL_gdp * p_f_old)
742        
743         print "HHH Y[0] =", Y[0]
744         print "HHH A_fw = ",A_fw
745         print "HHH A_fg= ",A_fg
746         print "HHH F_fg = ",F_fg
747         print "HHH F_fw = ",F_fw
748         print "HHH perm_f = ",self.perm_f
749         print "HHH k = ",(A_fw*self.perm_f[0])/D[0,0]
750         print "HHH k = ",(A_fw*self.perm_f[1])/D[0,0]
751         print "HHH X[0,:]= ",X[0,:]
752         print "HHH D[0,0] = ",D[0,0]
753         print "HHH D[0,1] = ",D[0,1]
754         print "HHH D[0,2] = ",D[0,2]
755        
756    
757          
758            
759         #print "HHH Y[1] = ",Y[1]
760         #print "HHH X[1,:] = ",X[1,:]
761         #print "HHH D[1,0] = ",D[1,0]
762         #print "HHH D[1,1] = ",D[1,1]
763         #print "HHH D[1,2] = ",D[1,2]
764         #print "HHH A_fg =",  A_fg
765            
766    
767    
768    
769         #1/0
770       self.__pde.setValue(A=A, D=D, X=X, Y=Y)       self.__pde.setValue(A=A, D=D, X=X, Y=Y)
771            
772       u2 = self.__pde.getSolution()       u2 = self.__pde.getSolution()
773         # this is deals with round-off errors to maintain physical meaningful values
774         # we should do this in a better way to detect values that are totally wrong
775         u2[1]=clip(u2[1],minval=0, maxval=1)
776         u2[2]=clip(u2[2],minval=0)
777         print "p _new =", u2[0]
778         #1/0
779       return u2       return u2

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

  ViewVC Help
Powered by ViewVC 1.1.26