/[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 3487 by gross, Mon Mar 28 03:44:22 2011 UTC
# 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 import sqrt, log, whereNegative, sup, inf
26    import math
27    
28  class MaterialProperty(object):  class MaterialProperty(object):
29     """     """
# Line 61  class MaterialPropertyWithDifferential(M Line 63  class MaterialPropertyWithDifferential(M
63        :remark: needs to be overwritten        :remark: needs to be overwritten
64        """        """
65        raise NotImplementedError        raise NotImplementedError
66    
67    class Porosity(MaterialPropertyWithDifferential):
68        """
69        defines porosity phi as function of pressure
70        
71             phi = phi_p_ref /(1 + X + X**2/2 ) with X= C * (p-p_ref)
72            
73        phi_p_ref is claculted from the initial porosity phi_0 at  pressure p_0
74        """
75        def __init__(self, phi_0, p_0, p_ref=1.*U.atm , C=4.0e-5/U.bar):
76          """
77          """
78          self.phi_p_ref=1 # will be overwritten later
79          self.p_ref=p_ref
80          self.C=C
81          # reset  phi_p_ref to get phi(p_0)=phi_0    
82          self.phi_p_ref=phi_0/self.getValue(p_0)
83                
84  class InterpolationTable(object):      def getValue(self, p):
85          """
86          returns the porosity for given pressure p
87          """
88          X= self.C * ( p - self.p_ref )      
89          return  self.phi_p_ref * (1. + X * (1. + X/2 ) )
90          
91        def getValueDifferential(self, p):
92          """
93          returns the porosity for given pressure p
94          """
95          X= self.C * ( p - self.p_ref )      
96          return  self.phi_p_ref * self.C * (1. + X )
97              
98          
99    class WaterDensity(MaterialPropertyWithDifferential):
100        """
101        set water density as
102          
103              rho = rho_ref (1 + X + X*X/2) with X= C * ( p - p_ref )
104              
105        with rho_ref =rho_s/B_ref * gravity
106        """
107        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):
108          self.rho_ref = rho_s/B_ref * gravity
109          self.C=C
110          self.p_ref=p_ref
111        
112        def getValue(self, p):
113          """
114          returns the density for given pressure p
115          """
116          X= self.C * ( p - self.p_ref )
117          return self.rho_ref * (1+ X * (1+ X/2) )  
118          
119        def getValueDifferential(self,  p):
120          """
121          """
122          X= self.C * ( p - self.p_ref )
123          return self.rho_ref * self.C * (1+ X)
124    
125    class WaterViscosity(MaterialProperty):
126      """
127      defines viscosity of water
128      
129      mu=mu_ref /(1 + X + X*X/2)
130      
131      with X=-C*(p-p_ref)
132      """
133    
134      def __init__(self, mu_ref = 0.3*U.cPoise, p_ref = 1.*U.atm , C = 0./U.bar):
135        """
136        set up
137        """
138        self.mu_ref = mu_ref
139        self.p_ref = p_ref
140        self.C=C
141      def getValue(self, p):
142          """
143          returns the viscosity for given pressure p
144          """
145          X= -self.C * ( p - self.p_ref )
146          return self.mu_ref/(1+ X * (1+ X/2) )  
147          
148    class GasDensity(MaterialPropertyWithDifferential):
149        """
150        set water density as
151          
152              rho = gravity * rho_air_s /B(p)
153              
154        where B is given by an interpolation table
155        """
156        def __init__(self, p, B, gravity=1., rho_air_s=1.2041*U.kg/U.m**3):
157          self.rho_ref =rho_air_s * gravity
158          self.tab=InterpolationTable(x=p, y=B)
159        
160        def getValue(self, p):
161          """
162          returns the density for given pressure p
163          """
164          return self.rho_ref/self.tab.getValue(p)
165          
166        def getValueDifferential(self,  p):
167          """
168          """
169          B    = self.tab.getValue(p)
170          dBdp = self.tab.getValueDifferential(p)
171          return  -self.rho_ref * dBdp/(B * B)
172    
173    class InterpolationTable(MaterialPropertyWithDifferential):
174     """     """
175     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
176     """     """
# Line 72  class InterpolationTable(object): Line 180  class InterpolationTable(object):
180        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
181        the value for x is set to y[0] for        the value for x is set to y[0] for
182        """        """
183          MaterialPropertyWithDifferential.__init__(self)
184        if len(x) != len(y):        if len(x) != len(y):
185       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."
186        if len(x) > 0:        if len(x) < 1 :
187       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."
188                
189        x_ref=x[0]        x_ref=x[0]
# Line 85  class InterpolationTable(object): Line 194  class InterpolationTable(object):
194        self.__x = x        self.__x = x
195        self.__y = y        self.__y = y
196        self.__obeyBounds = obeyBounds        self.__obeyBounds = obeyBounds
197          
198     def getValue(self, x):     def getValue(self, x):
199        """        """
200        returns the interpolated values of x        returns the interpolated values of x
201        """        """
202        x0=self.__x[0]        X=self.__x
203          Y=self.__y
204          
205          x0=X[0]
206        m0=whereNegative(x-x0)        m0=whereNegative(x-x0)
207        if self.__obeyBounds:        if self.__obeyBounds:
208       if sup(m0) > 0:       if sup(m0) > 0:
209          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])
210        out=self.__x[0]        out=self.__x[0]
211        for i in range(1,len(self.__x)):        for i in range(1,len(X)):
212        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]  
213        out = out * m0 + z * (1-m0)        out = out * m0 + z * (1-m0)
214        m0=whereNegative(x-x[i])        m0=whereNegative(x-X[i])
215            
216        if self.__obeyBounds:        if self.__obeyBounds:
217          if inf(m0) < 1:          if inf(m0) < 1:
218             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])
219        else:        else:
220          out = out * m0 + y[-1] * (1-m0)          out = out * m0 + V[-1] * (1-m0)
221        return out        return out
222    
223       def getValueDifferential(self, x):
224          X=self.__x
225          Y=self.__y
226    
227          x0=X[0]
228          m0=whereNegative(x-x0)
229          if self.__obeyBounds:
230         if sup(m0) > 0:
231            raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
232          out=0.
233          for i in range(1,len(X)):
234          z=(Y[i]-Y[i-1])/(X[i]-X[i-1])
235          out = out * m0 + z * (1-m0)
236          m0=whereNegative(x-X[i])
237        
238          if self.__obeyBounds:
239            if inf(m0) < 1:
240               raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
241          else:
242            out = out * m0
243          return out  
244          
245    
246  class Well(object):  class Well(object):
247     """     """
248     generic well     generic well
# Line 117  class Well(object): Line 252  class Well(object):
252     """     """
253     WATER="Water"     WATER="Water"
254     GAS="Gas"     GAS="Gas"
255     def __init__(self,  *args, **kwargs):     def __init__(self,  Q=0., schedule=[0.], phase="Water", *args, **kwargs):
256         """         """
257         set up well         set up well
258         """         """
259         pass         self.__schedule=schedule
260           self.__phase=phase
261           self.__Q = Q
262          
263     def getGeometry(self):     def getGeometry(self):
264         """         """
265         returns a function representing the geometry of the well.         returns a function representing the geometry of the well.
# Line 140  class Well(object): Line 278  class Well(object):
278         """         """
279         raise NotImplementedError         raise NotImplementedError
280        
281     def getFlowRate(self):     def getFlowRate(self, t):
282        """        """
283        returns the flow rate        returns the flow rate
284        """        """
285        return 0.        if self.isOpen(t):
286              return self.__Q
287          else:
288        return 0.
289        
290     def getBHP(self):     def getBHP(self):
291        """        """
292        return bottom-hole pressure        return bottom-hole pressure
# Line 158  class Well(object): Line 299  class Well(object):
299        """        """
300        returns the pahse the well works on        returns the pahse the well works on
301        """        """
302        return self.WATER        return self.__phase
303          
304       def isOpen(self, t):
305         """
306         returns True is the well is open at time t
307         """
308         out = False
309         t0=min(t, self.__schedule[0])
310         i=0
311         while t < t0:
312           if out:
313         out=False
314           else:
315         out=True  
316           i+=1
317           if i < len(self.__schedule[i]):
318          t0=self.__schedule[i]
319           else:
320          t0=t
321         return out
322    
323  class PeacemanWell(Well):  class VerticalPeacemanWell(Well):
324     """     """
325     defines a well using Peaceman's formula     defines a well using Peaceman's formula
326     """     """
327     pass     def __init__(self,name, schedule = [0.], BHP=1.*U.atm, Q=0, r=10*U.cm, X0=[0.,0.,0.], D=[1.*U.km,1.*U.km, 1.*U.m],
328                 perm=[1.*U.cPoise,1.*U.cPoise, 1.*U.cPoise], phase=Well.WATER, s=0):
329           """
330           set up well
331          
332           :param BHP: ottom-hole pressure
333           :param Q: flow rate
334           :param r: well radius
335           :param X: location
336           :param D: dimension of well block
337           :param perm: permeabilities
338           :param s: skin factor
339           """
340           Well.__init__(self, Q=Q, schedule=schedule, phase=phase,)
341           r_el=0.28 * ( (perm[1]/perm[0])**0.5 * D[0]**2 +  (perm[0]/perm[1])**0.5 * D[1]**2 ) \
342                             / ( (perm[1]/perm[0])**0.25 + (perm[1]/perm[0])**0.25 )
343           self.__PI = 2 * math.pi * D[2] * sqrt(perm[1]*perm[0]) / (log(r_el/r) + s)
344           self.__r = r
345    
346           self.__BHP = BHP
347           self.__D=D
348           self.r_el=r_el
349           self.X0=X0
350          
351       def getGeometry(self):
352           """
353           returns a function representing the geometry of the well.
354           typically a Gaussian profile around the location of the well.
355          
356           :note: needs to be overwritten
357           """
358           raise NotImplementedError
359    
360       def getProductivityIndex(self,t):
361           """
362           returns the productivity index of the well.
363           typically a Gaussian profile around the location of the well.
364           """
365           return self.__PI
366          
367       def getBHP(self):
368          """
369          return bottom-hole pressure
370          """
371          return self.__BHP
372    
373  class DualPorosity(object):  class DualPorosity(object):
374     """     """
# Line 188  class DualPorosity(object): Line 392  class DualPorosity(object):
392        :type phi_m: `MaterialPropertyWithDifferential` or None        :type phi_m: `MaterialPropertyWithDifferential` or None
393        :param phi: total porosity if None phi=1.        :param phi: total porosity if None phi=1.
394        :type phi: scalar or None        :type phi: scalar or None
395        :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
396        :type S_fg: `MaterialPropertyWithDifferential`        :type S_fg: `MaterialPropertyWithDifferential`
397        :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
398        :type S_mg: `MaterialPropertyWithDifferential`        :type S_mg: `MaterialPropertyWithDifferential`
# Line 235  class DualPorosity(object): Line 439  class DualPorosity(object):
439        self.rho_w = rho_w        self.rho_w = rho_w
440        self.rho_g = rho_g            self.rho_g = rho_g    
441        self.wells=wells        self.wells=wells
442          self.t =0
443          
444          self.__iter_max=1
445          self.__rtol=1.e-4
446          self.__verbose=False
447       def setIterationControl(self, iter_max=None, rtol=None, verbose=None):
448         """
449         sets parameters to control iteration process
450         """
451         if iter_max !=None: self.__iter_max=iter_max
452         if rtol !=None: self.__rtol = rtol
453         if verbose !=None: self.__verbose=verbose
454        
455       def update(self, dt):
456             self.u, self.u_old = self.u.copy(), self.u
457             n=0
458             rerr=1.
459             while n < self.__iter_max and rerr > self.__rtol:
460                u=self.solvePDE(dt)
461            norm_u=Lsup(u)
462            norm_e=Lsup(u-self.u)
463            if norm_u > 0:
464                rerr=norm_e/norm_u
465            else:
466                rerr=norm_e
467            if self.__verbose: print "iteration step %d: relative change = %e"%(n, rerr)
468            n+=1
469            self.u=u
470             print "iteration completed."
471             self.t+=dt
472    
473    
474  class PorosityOneHalfModel(DualPorosity):  class PorosityOneHalfModel(DualPorosity):
# Line 247  class PorosityOneHalfModel(DualPorosity) Line 478  class PorosityOneHalfModel(DualPorosity)
478        This corresponds to the coal bed methan model in the ECLIPSE code.        This corresponds to the coal bed methan model in the ECLIPSE code.
479        """        """
480    
481        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,
482               perm_f_0=None, perm_f_1=None, perm_f_2=None,               perm_f_0=None, perm_f_1=None, perm_f_2=None,
483               k_w =None, k_g=None, mu_w =None, mu_g =None,               k_w =None, k_g=None, mu_w =None, mu_g =None,
484               rho_w =None, rho_g=None,               rho_w =None, rho_g=None,
# Line 263  class PorosityOneHalfModel(DualPorosity) Line 494  class PorosityOneHalfModel(DualPorosity)
494       :type phi: scalar or None       :type phi: scalar or None
495       :param L_g: gas adsorption as function of gas pressure       :param L_g: gas adsorption as function of gas pressure
496       :type L_g: `MaterialPropertyWithDifferential`       :type L_g: `MaterialPropertyWithDifferential`
497       :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
498       :type S_fg: `MaterialPropertyWithDifferential`       :type S_fg: `MaterialPropertyWithDifferential`
499       :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
500       :type perm_f_0: scalar       :type perm_f_0: scalar
# Line 289  class PorosityOneHalfModel(DualPorosity) Line 520  class PorosityOneHalfModel(DualPorosity)
520        
521       DualPorosity.__init__(self, domain,       DualPorosity.__init__(self, domain,
522                    phi_f=phi_f, phi_m=None, phi=phi,                    phi_f=phi_f, phi_m=None, phi=phi,
523                    S_fg=S_fg, S_mg=None,                    S_fg=None, S_mg=None,
524                    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,
525                    perm_m_0=None , perm_m_1=None, perm_m_2=None,                    perm_m_0=None , perm_m_1=None, perm_m_2=None,
526                    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,
# Line 305  class PorosityOneHalfModel(DualPorosity) Line 536  class PorosityOneHalfModel(DualPorosity)
536       """       """
537       return self.__pde.getSolverOptions()       return self.__pde.getSolverOptions()
538            
539        def solvePDE(self):        def getState(self):
540             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]  
541    
542          def getOldState(self):
543         return self.u_old[0], self.u_old[1],  self.u_old[2]
544    
545       phi_f   =self.phi_f.getValue(p_fg)        def setInitialState(self, p=1.*U.atm, S_fg=0,  C_mg=None):
546       dphi_fdp=self.phi_f.getValueDifferential(p_fg)          """
547            sets the initial state
548            
549            :param p: pressure
550            :param S_fg: gas saturation in fractured rock
551            :param C_mg: gas concentration in coal matrix. if not given it is calculated
552                using the  gas adsorption curve.
553            """    
554            self.u=self.__pde.createSolution()
555            self.u[0]=p
556            self.u[1]=S_fg
557            if C_mg == None:
558              self.u[2]= self.L_g(p)
559            else:
560              self.u[2]=C_mg
561          
562          def solvePDE(self, dt):
563        
564         p_f, S_fg, C_mg=self.getState()
565         p_f_old, S_fg_old, C_mg_old=self.getOldState()
566    
567             S_fw=1-S_fg
568    
569             k_fw=self.k_w(S_fw)
570             k_fg=self.k_g(S_fg)
571             mu_fw=self.mu_w(p_f)
572             mu_fg=self.mu_g(p_f)
573            
574    
575         phi_f   =self.phi_f.getValue(S_fg)
576         dphi_fdp=self.phi_f.getValueDifferential(S_fg)
577            
      S_fg=  self.S_fg.getValue(p_fc)  
      dS_fgdp=   self.S_fg.getValueDifferential(p_fc)  
578       S_fw=  1-S_fg       S_fw=  1-S_fg
579            
580       rho_fw     = self.rho_w.getValue(p_fw)       rho_fw     = self.rho_w.getValue(p_f)
581       drho_fwdp  = self.rho_w.getValueDifferential(p_fw)       drho_fwdp  = self.rho_w.getValueDifferential( p_f)
582       rho_fg = self.rho_g.getValue(p_fg)       rho_fg = self.rho_g.getValue(p_f)
583       drho_fgdp  = self.rho_g.getValueDifferential(p_fg)       drho_fgdp  = self.rho_g.getValueDifferential(p_f)
584        
585         L_g       = self.L_g.getValue(p_f)
586         dL_gdp =  self.rho_w.getValueDifferential(p_f)
587            
      L_g       = self.getValue(p_fg)  
      dL_gdp_fg =  self.rho_w.getValueDifferential(p_fg)  
588            
589       A_fw = rho_fw * k_fw/mu_fw       A_fw = rho_fw * k_fw/mu_fw
590       A_fg = rho_fg * k_fg/mu_fg       A_fg = rho_fg * k_fg/mu_fg
# Line 340  class PorosityOneHalfModel(DualPorosity) Line 593  class PorosityOneHalfModel(DualPorosity)
593       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 )
594            
595            
596       E_fww= phi_f * ( rho_fw * dS_fgdp + S_fw * drho_fwdp )       E_fpp= S_fw * (  rho_fw * dphi_fdp + phi_f  * drho_fwdp )
597       E_fwg= rho_fw * ( S_fw * dphi_fdp - phi_f * dS_fgdp )       E_fps=  -  phi_f * rho_fw
598       E_fgw= - phi_f * rho_fg * dS_fgdp       E_fsp= S_fg *( rho_fg * dphi_fdp - phi_f *  drho_fgdp )
599       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
600    
601    
602            
# Line 367  class PorosityOneHalfModel(DualPorosity) Line 620  class PorosityOneHalfModel(DualPorosity)
620             Q_w = Q_w + Pi_I * P_I * chi_I             Q_w = Q_w + Pi_I * P_I * chi_I
621             R_w = R_w + Pi_I       * chi_I             R_w = R_w + Pi_I       * chi_I
622                        
623       F_fw_Y = A_fw * Q_w       F_fp_Y = A_fw * Q_w
624       F_fg_Y = A_fg * Q_g       F_fs_Y = A_fg * Q_g
625       D_fgg = - A_fg * R_g       D_fss = - A_fg * R_g
626       D_fww = - A_fw * R_w       D_fpp = - A_fw * R_w
627            
628       if self.domain.getDim() == 3:       if self.domain.getDim() == 3:
629          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]
630          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]
631       else:       else:
632          F_fw_X = 0 * kronecker(self.domain)[1]          F_fp_X = 0 * kronecker(self.domain)[1]
633          F_fg_X = 0 * kronecker(self.domain)[1]          F_fs_X = 0 * kronecker(self.domain)[1]
634                    
635       F_mg = B * ( L_g - dL_gdp_fg * p_fg )       F_mg = B * ( L_g - dL_gdp * S_fg )
636    
637            
638       D=self.__pde.getNewCoefficient("D")       D=self.__pde.getNewCoefficient("D")
# Line 390  class PorosityOneHalfModel(DualPorosity) Line 643  class PorosityOneHalfModel(DualPorosity)
643       dtXI = dt*self.XI       dtXI = dt*self.XI
644       dtcXI = dt*(1-self.XI)       dtcXI = dt*(1-self.XI)
645            
646       D[0,0]=E_fww + dtXI * D_fww       D[0,0]=E_fpp + dtXI * D_fpp
647       D[0,1]=E_fwg       D[0,1]=E_fps
648       D[1,0]=E_fgw       D[1,0]=E_fsp
649       D[1,1]=E_fgg + dtXI * D_fgg       D[1,1]=E_fss + dtXI * D_fss
650       D[1,2]=rho_fg       D[1,2]=rho_fg
651       D[2,1]= dtXI * B * dL_gdp_fg       D[2,1]= dtXI * B * dL_gdp
652       D[2,2]= 1 + dtXI * B       D[2,2]= 1 + dtXI * B
653            
654       H_fw = dtcXI * A_fw * grad(p_fw_old)       H_fw = dtcXI * A_fw * grad( p_f_old)
655       H_fg = dtcXI * A_fg * grad(p_fg_old)       H_fg = dtcXI * A_fg * grad(S_fg_old)
656            
657       for i in range(self.domain.getDim()):       for i in range(self.domain.getDim()):
658          A[0,i,0,i] = dtXI * A_fw * self.perm_f[i]          A[0,i,0,i] = dtXI * A_fw * self.perm_f[i]
659          A[1,i,1,i] = dtXI * A_fg * self.perm_f[i]          A[1,i,1,i] = dtXI * A_fg * self.perm_f[i]
660                    
661          X[0,i]= dt * F_fw_X - H_fw * self.perm_f[i]          X[0,i]= dt * F_fp_X - H_fw * self.perm_f[i]
662          X[1,i]= dt * F_fg_X - H_fg * self.perm_f[i]          X[1,i]= dt * F_fs_X - H_fg * self.perm_f[i]
663                    
664       Y[0] = E_fww * p_fw_old + E_fwg * p_fg_old +                     dt * F_fw_Y - dtcXI * D_fww * p_fw_old       Y[0] = E_fpp *  p_f_old + E_fps * S_fg_old +                     dt * F_fp_Y - dtcXI * D_fpp *  p_f_old
665       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[1] = E_fsp *  p_f_old + E_fss * S_fg_old + rho_fg * C_mg_old + dt * F_fs_Y - dtcXI * D_fss * S_fg_old
666       Y[2] = C_mg_old                                                + dt * F_mg_Y - dtcXI * ( dL_gdp_fg * p_fg_old + C_mg_old)       Y[2] = C_mg_old                                                + dt * F_mg_Y - dtcXI * ( dL_gdp * S_fg_old + C_mg_old)
667            
668       self.__pde.setValue(A=A, D=D, X=X, Y=Y)       self.__pde.setValue(A=A, D=D, X=X, Y=Y)
669            

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

  ViewVC Help
Powered by ViewVC 1.1.26