# Diff of /trunk/finley/test/python/coalgas.py

revision 3478 by gross, Wed Mar 23 04:02:39 2011 UTC revision 3502 by gross, Thu Apr 28 05:06:24 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)
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, maximum, exp
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        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        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 + X[-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, decay_factor = 5.):
368                 # reset_factor=0.1):
369           """
370           set up well
371
372           :param BHP: ottom-hole pressure
373           :param Q: flow rate
375           :param X: location
376           :param D: dimension of well block
377           :param perm: permeabilities
378           :param s: skin factor
379           """
380           Well.__init__(self, name, domain, Q=Q, schedule=schedule, phase=phase,BHP_limit=BHP_limit)
381           r_el=0.28 * sqrt( sqrt(perm[1]/perm[0]) * D[0]**2 +  sqrt(perm[0]/perm[1]) * D[1]**2 )\
382                             / ( (perm[1]/perm[0])**0.25 + (perm[1]/perm[0])**0.25 )
383           self.__PI = 2 * math.pi * D[2] * sqrt(perm[1]*perm[0]) / (log(r_el/r) + s)
384           self.__r = r
385
386           self.__D=D
387           self.r_el=r_el
388           self.X0=X0[:self.domain.getDim()]
389
390           x=Function(domain).getX()
391           self.chi = 1.
392           decay_factor
393           for i in range(domain.getDim()):
394            r=decay_factor * D[i]
395            self.chi = self.chi * clip( (D[i]/2.+r-abs(X0[i]-x[i]))/r, maxval=1., minval=0.)
396            #if reset_factor > 0:
397            #     f=maximum(0., abs(x[i]-X0[i])-D[i]/2*(1.-reset_factor))
398            #     self.chi = self.chi * exp(-(f/(reset_factor*D[i]/2))**2/2)
399            #else:
400                #     self.chi = self.chi * whereNegative(abs(x[i]-X0[i])-D[i]/2)
401
402           l=integrate(self.chi)
403           self.chi*=1./l
404
405           #self.chi=0.00000001
406       def getLocation(self):
407           return self.X0
408
409       def getGeometry(self):
410           """
411           returns a function representing the geometry of the well.
412           typically a Gaussian profile around the location of the well.
413
414           :note: needs to be overwritten
415           """
416           if self.domain.getDim() <3:
417              return self.chi/self.__D[2]
418           else:
419          return self.chi
420
421       def getProductivityIndex(self):
422           """
423           returns the productivity index of the well.
424           """
425           return self.__PI
426
427
428  class DualPorosity(object):  class DualPorosity(object):
429     """     """
# Line 188  class DualPorosity(object): Line 447  class DualPorosity(object):
447        :type phi_m: `MaterialPropertyWithDifferential` or None        :type phi_m: `MaterialPropertyWithDifferential` or None
448        :param phi: total porosity if None phi=1.        :param phi: total porosity if None phi=1.
449        :type phi: scalar or None        :type phi: scalar or None
450        :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
451        :type S_fg: `MaterialPropertyWithDifferential`        :type S_fg: `MaterialPropertyWithDifferential`
452        :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
453        :type S_mg: `MaterialPropertyWithDifferential`        :type S_mg: `MaterialPropertyWithDifferential`
# Line 235  class DualPorosity(object): Line 494  class DualPorosity(object):
494        self.rho_w = rho_w        self.rho_w = rho_w
495        self.rho_g = rho_g            self.rho_g = rho_g
496        self.wells=wells        self.wells=wells
497          self.t =0
498
499          self.__iter_max=1
500          self.__rtol=1.e-4
501          self.verbose=False
502          self.XI=0.5
503       def setIterationControl(self, iter_max=None, rtol=None, verbose=None, xi=None):
504         """
505         sets parameters to control iteration process
506         """
507         if iter_max !=None: self.__iter_max=iter_max
508         if rtol !=None: self.__rtol = rtol
509         if verbose !=None: self.verbose=verbose
510         if xi !=None: self.XI=xi
511
512       def update(self, dt):
513             self.u, self.u_old = self.u.copy(), self.u
514             n=0
515             rerr=1.
516             while n < self.__iter_max and rerr > self.__rtol:
517                u=self.solvePDE(dt)
518            norm_u=Lsup(u)
519            norm_e=Lsup(u-self.u)
520
521            if norm_u > 0:
522                rerr=norm_e/norm_u
523            else:
524                rerr=norm_e
525            if self.verbose: print "iteration step %d: relative change = %e"%(n, rerr)
526            n+=1
527            self.u=u
528             print "iteration completed."
529             self.t+=dt
530
531
532  class PorosityOneHalfModel(DualPorosity):  class PorosityOneHalfModel(DualPorosity):
# Line 247  class PorosityOneHalfModel(DualPorosity) Line 536  class PorosityOneHalfModel(DualPorosity)
536        This corresponds to the coal bed methan model in the ECLIPSE code.        This corresponds to the coal bed methan model in the ECLIPSE code.
537        """        """
538
539        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,
540               perm_f_0=None, perm_f_1=None, perm_f_2=None,               perm_f_0=None, perm_f_1=None, perm_f_2=None,
541               k_w =None, k_g=None, mu_w =None, mu_g =None,               k_w =None, k_g=None, mu_w =None, mu_g =None,
542               rho_w =None, rho_g=None,               rho_w =None, rho_g=None, sigma=0, A_mg=0, f_rg=1.,
543                 wells=[]):                 wells=[]):
544       """       """
545       set up       set up
# Line 263  class PorosityOneHalfModel(DualPorosity) Line 552  class PorosityOneHalfModel(DualPorosity)
552       :type phi: scalar or None       :type phi: scalar or None
553       :param L_g: gas adsorption as function of gas pressure       :param L_g: gas adsorption as function of gas pressure
554       :type L_g: `MaterialPropertyWithDifferential`       :type L_g: `MaterialPropertyWithDifferential`
555       :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
556       :type S_fg: `MaterialPropertyWithDifferential`       :type S_fg: `MaterialPropertyWithDifferential`
557       :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
558       :type perm_f_0: scalar       :type perm_f_0: scalar
# Line 285  class PorosityOneHalfModel(DualPorosity) Line 574  class PorosityOneHalfModel(DualPorosity)
574       :type rho_g: `MaterialPropertyWithDifferential`       :type rho_g: `MaterialPropertyWithDifferential`
575       :param wells : list of wells       :param wells : list of wells
576       :type wells: list of `Well`       :type wells: list of `Well`
577         :param sigma: shape factor for gas matrix diffusion
578         :param A_mg: diffusion constant for gas adsorption
579         :param f_rg: gas re-adsorption factor
580       """       """
581
582       DualPorosity.__init__(self, domain,       DualPorosity.__init__(self, domain,
583                    phi_f=phi_f, phi_m=None, phi=phi,                    phi_f=phi_f, phi_m=None, phi=phi,
584                    S_fg=S_fg, S_mg=None,                    S_fg=None, S_mg=None,
585                    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,
586                    perm_m_0=None , perm_m_1=None, perm_m_2=None,                    perm_m_0=None , perm_m_1=None, perm_m_2=None,
587                    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,
588                    rho_w =rho_w, rho_g=rho_g,                    rho_w =rho_w, rho_g=rho_g,
589                    wells=wells)                    wells=wells)
590       self.L_g=L_g       self.L_g=L_g
591         self.sigma = sigma
592         self.A_mg = A_mg
593         self.f_rg  = f_rg
594       self.__pde=LinearPDE(self.domain, numEquations=3, numSolutions =3)       self.__pde=LinearPDE(self.domain, numEquations=3, numSolutions =3)
595
596
# Line 305  class PorosityOneHalfModel(DualPorosity) Line 600  class PorosityOneHalfModel(DualPorosity)
600       """       """
601       return self.__pde.getSolverOptions()       return self.__pde.getSolverOptions()
602
603        def solvePDE(self):        def getState(self):
604             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]
605
606          def getOldState(self):
607         return self.u_old[0], self.u_old[1],  self.u_old[2]
608
609       phi_f   =self.phi_f.getValue(p_fg)        def setInitialState(self, p=1.*U.atm, S_fg=0,  C_mg=None):
610       dphi_fdp=self.phi_f.getValueDifferential(p_fg)          """
611                sets the initial state
612       S_fg=  self.S_fg.getValue(p_fc)
613       dS_fgdp=   self.S_fg.getValueDifferential(p_fc)          :param p: pressure
614       S_fw=  1-S_fg          :param S_fg: gas saturation in fractured rock
615                :param C_mg: gas concentration in coal matrix. if not given it is calculated
616       rho_fw     = self.rho_w.getValue(p_fw)              using the  gas adsorption curve.
617       drho_fwdp  = self.rho_w.getValueDifferential(p_fw)          """
618       rho_fg = self.rho_g.getValue(p_fg)          self.u=self.__pde.createSolution()
619       drho_fgdp  = self.rho_g.getValueDifferential(p_fg)          self.u[0]=p
620            self.u[1]=S_fg
621            if C_mg == None:
622              self.u[2]= self.L_g(p)*self.rho_g.getFormationFactor(p)
623            else:
624              self.u[2]=C_mg
625
626          def solvePDE(self, dt):
627
628         C_couple=1
629
630         p_f, S_fg, C_mg=self.getState()
631         p_f_old, S_fg_old, C_mg_old=self.getOldState()
632
633             S_fw=1-S_fg
634
635         if self.verbose:
636              print "p_f range = ",inf(p_f),sup(p_f)
637              print "S_fg range = ",inf(S_fg),sup(S_fg)
638              print "S_fw range = ",inf(S_fw),sup(S_fw)
639              print "C_mg range = ",inf(C_mg),sup(C_mg)
640
641             k_fw=self.k_w(S_fw)
642             if self.verbose: print "k_fw range = ",inf(k_fw),sup(k_fw)
643
644             k_fg=self.k_g(S_fg)
645             if self.verbose: print "k_fg range = ",inf(k_fg),sup(k_fg)
646
647         mu_fw=self.mu_w(p_f)
648             if self.verbose: print "mu_fw range = ",inf(mu_fw),sup(mu_fw)
649
650         mu_fg=self.mu_g(p_f)
651             if self.verbose: print "mu_fg range = ",inf(mu_fg),sup(mu_fg)
652
653
654         phi_f   =self.phi_f.getValue(p_f)
655         dphi_fdp=self.phi_f.getValueDifferential(p_f)
656             if self.verbose: print "phi_f range = ",inf(phi_f),sup(phi_f)," (slope %e,%e)"%(inf(dphi_fdp), sup(dphi_fdp))
657
658         rho_fw     = self.rho_w.getValue(p_f)
659         drho_fwdp  = self.rho_w.getValueDifferential(p_f)
660             if self.verbose: print "rho_fw range = ",inf(rho_fw),sup(rho_fw)," (slope %e,%e)"%(inf(drho_fwdp), sup(drho_fwdp))
661
662             rho_fg = self.rho_g.getValue(p_f)
663         drho_fgdp  = self.rho_g.getValueDifferential(p_f)
664             if self.verbose: print "rho_fg range = ",inf(rho_fg),sup(rho_fg)," (slope %e,%e)"%(inf(drho_fgdp), sup(drho_fgdp))
665
666
667
668         L_g_0       = self.L_g(p_f)
669         FF_g=self.rho_g.getFormationFactor(p_f)
670         L_g = L_g_0 * FF_g
671
672             if self.verbose:
673              print "L_g_0 range = ",inf(L_g_0),sup(L_g_0)
674              print "FF_g range = ",inf(FF_g),sup(FF_g)
675              print "L_g range = ",inf(L_g),sup(L_g)
676
L_g       = self.getValue(p_fg)
dL_gdp_fg =  self.rho_w.getValueDifferential(p_fg)
677
678       A_fw = rho_fw * k_fw/mu_fw       A_fw = rho_fw * k_fw/mu_fw
679       A_fg = rho_fg * k_fg/mu_fg       A_fg = rho_fg * k_fg/mu_fg
680
681       m = whereNegative(L_g - C_mg)         m = whereNegative(L_g - C_mg)
682       B = self.sigma * self.a_mg * (m + (1-m) * self.f_rg * S_fg )       H = self.sigma * self.A_mg * (m + (1-m) * self.f_rg * S_fg )
683             print "XXXX H =",H
684             print "XXXX self.sigma  =",self.sigma
685       E_fww= phi_f * ( rho_fw * dS_fgdp + S_fw * drho_fwdp )       print "XXXX self.A_mg  =",self.A_mg
686       E_fwg= rho_fw * ( S_fw * dphi_fdp - phi_f * dS_fgdp )
687       E_fgw= - phi_f * rho_fg * dS_fgdp       E_fpp= S_fw * (  rho_fw * dphi_fdp + phi_f  * drho_fwdp )
688       E_fgg= S_fg * rho_fg  * dphi_fdp + phi_f * rho_fg * dS_fgdp + phi_f  * S_fg * drho_fgdp       E_fps=  -  phi_f * rho_fw
689         E_fsp= S_fg * rho_fg * dphi_fdp + (phi_f * S_fg + C_mg) * drho_fgdp
690         E_fss= phi_f * rho_fg
691
692       Q_w=0.       F_fw=0.
693       Q_g=0.       F_fg=0.
694       R_w=0.       D_fw=0.
695       R_g=0.       D_fg=0.
696
697       for I in self.wells:       for I in self.wells:
698          chi_I= I.getGeometry()            chi_I= I.getGeometry()
699          Pi_I = I.getProductivityIndex()            Pi_I = I.getProductivityIndex()
700          P_I = I.getBHP()            BHP_limit_I=I.getBHPLimit()
701          Q_I = I.getFlowRate()            p_f_I=integrate(p_f*I.chi)
702                      S_fg_I=integrate(S_fg*I.chi)
703          if self.getPhase == Well.WATER:
704             Q_w = Q_w + Q_I        * chi_I            A_fw_I= self.rho_w.getValue(p_f_I) * self.k_w(1-S_fg_I)/self.mu_w(p_f_I)
705             Q_g = Q_g + Pi_I * P_I * chi_I            A_fg_I= self.rho_g.getValue(p_f_I) * self.k_g(S_fg_I)/self.mu_g(p_f_I)
706             R_g = R_g + Pi_I       * chi_I
707          else:            if I.isOpen(self.t+dt):
708             Q_g = Q_g + Q_I        * chi_I          if self.verbose: print "well %s is open."%I.name
709             Q_w = Q_w + Pi_I * P_I * chi_I          if I.getPhase() == Well.WATER:
710             R_w = R_w + Pi_I       * chi_I              q_I = self.rho_w.rho_0 * I.getFlowRate()
711                          p_I_ref=p_f_I-q_I/(A_fw_I * Pi_I)
712       F_fw_Y = A_fw * Q_w          else:
713       F_fg_Y = A_fg * Q_g              q_I = self.rho_g.rho_0 * I.getFlowRate()
714       D_fgg = - A_fg * R_g              p_I_ref=p_f_I-q_I/(A_fg_I * Pi_I)
715       D_fww = - A_fw * R_w
716            print "ZZZ =",p_f_I, q_I, self.rho_w.rho_0, I.getFlowRate(), A_fw_I, Pi_I, q_I/(A_fw_I * Pi_I)
717
718            if BHP_limit_I > p_I_ref:
719              BHP_I=BHP_limit_I
720              D_fw = D_fw + Pi_I * A_fw_I *        chi_I
721              F_fw = F_fw - Pi_I * A_fw_I * BHP_I *chi_I
722              D_fg = D_fg + Pi_I * A_fg_I *        chi_I
723              F_fg = F_fg - Pi_I * A_fg_I * BHP_I *chi_I
724            else:
725              BHP_I=p_I_ref
726              if I.getPhase() == Well.WATER:
727                  F_fw = F_fw +  q_I  * chi_I
728                  D_fg = D_fg + Pi_I * A_fg_I *         chi_I
729                  F_fg = F_fg - Pi_I * A_fg_I * BHP_I * chi_I
730              else:
731                  F_fg = F_fg +  q_I  * chi_I
732                  D_fw = D_fw + Pi_I * A_fw_I *        chi_I
733                  F_fw = F_fw - Pi_I * A_fw_I * BHP_I *chi_I
734
735              else:
736              if self.verbose: print "well %s is closed."%I.name
737              BHP_I=p_f_I
738
739              if self.verbose: print "well %s:BHP = %e"%(I.name, BHP_I)
740              I.setBHP(BHP_I)
741
742         F_fp_Y = - F_fw
743         F_fs_Y = - F_fg
744         D_fpp =  D_fw
745         D_fsp =  D_fg
746
747
748       if self.domain.getDim() == 3:       if self.domain.getDim() == 3:
749          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]
750          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]
751       else:       else:
752          F_fw_X = 0 * kronecker(self.domain)[1]          F_fp_X = 0 * kronecker(self.domain)[1]
753          F_fg_X = 0 * kronecker(self.domain)[1]          F_fs_X = 0 * kronecker(self.domain)[1]
754
755       F_mg = B * ( L_g - dL_gdp_fg * p_fg )           F_mg_Y = -H * L_g
756
757
758       D=self.__pde.getNewCoefficient("D")       D=self.__pde.createCoefficient("D")
759       A=self.__pde.getNewCoefficient("A")       A=self.__pde.createCoefficient("A")
760       Y=self.__pde.getNewCoefficient("Y")       Y=self.__pde.createCoefficient("Y")
761       X=self.__pde.getNewCoefficient("X")       X=self.__pde.createCoefficient("X")
762
763       dtXI = dt*self.XI       dtXI = dt*self.XI
764       dtcXI = dt*(1-self.XI)       dtcXI = dt*(1-self.XI)
765
766         D[0,0]=E_fpp + dtXI * D_fpp
767         D[0,1]=E_fps
768         D[1,0]=E_fsp + dtXI * D_fsp
769         D[1,1]=E_fss
770         D[1,2]=rho_fg * C_couple
771         #D[2,0]= - dtXI * B * dL_gdp
772         D[2,2]= 1 - dtXI * H
773
D[0,0]=E_fww + dtXI * D_fww
D[0,1]=E_fwg
D[1,0]=E_fgw
D[1,1]=E_fgg + dtXI * D_fgg
D[1,2]=rho_fg
D[2,1]= dtXI * B * dL_gdp_fg
D[2,2]= 1 + dtXI * B

H_fw = dtcXI * A_fw * grad(p_fw_old)
H_fg = dtcXI * A_fg * grad(p_fg_old)
774
775       for i in range(self.domain.getDim()):       for i in range(self.domain.getDim()):
776          A[0,i,0,i] = dtXI * A_fw * self.perm_f[i]          A[0,i,0,i] = dtXI * A_fw * self.perm_f[i]
777          A[1,i,1,i] = dtXI * A_fg * self.perm_f[i]          A[1,i,1,i] = dtXI * A_fg * self.perm_f[i]
778
779          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()]
780          X[1,i]= dt * F_fg_X - H_fg * self.perm_f[i]           X[0,:]=  - A_fw * g  + dt * F_fp_X
781                 X[1,:]=  - A_fg * g  + dt * F_fs_X
782       Y[0] = E_fww * p_fw_old + E_fwg * p_fg_old +                     dt * F_fw_Y - dtcXI * D_fww * p_fw_old
783       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
784       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
785             Y[2] = C_mg_old                                                           + dt * F_mg_Y + dtcXI * H * C_mg_old
786
787         print "HHH Y[0] =", Y[0]
788         print "HHH A_fw = ",A_fw
789         print "HHH A_fg= ",A_fg
790         print "HHH F_fg = ",F_fg
791         print "HHH F_fw = ",F_fw
792         print "HHH perm_f = ",self.perm_f
793         print "HHH k = ",(A_fw*self.perm_f[0])/D[0,0]
794         print "HHH k = ",(A_fw*self.perm_f[1])/D[0,0]
795         print "HHH X[0,:]= ",X[0,:]
796         print "HHH D[0,0] = ",D[0,0]
797         print "HHH D[0,1] = ",D[0,1]
798         print "HHH D[0,2] = ",D[0,2]
799
800       self.__pde.setValue(A=A, D=D, X=X, Y=Y)       self.__pde.setValue(A=A, D=D, X=X, Y=Y)
801
802       u2 = self.__pde.getSolution()       u2 = self.__pde.getSolution()
803         # this is deals with round-off errors to maintain physical meaningful values
804         # we should do this in a better way to detect values that are totally wrong
805         u2[1]=clip(u2[1],minval=0, maxval=1)
806         u2[2]=clip(u2[2],minval=0)
807
808         # update well production
809         p_f=u2[0]
810         S_fg=u2[1]
811         for I in self.wells:
812              Pi_I = I.getProductivityIndex()
813              q_I = I.getFlowRate()
814              p_f_I=integrate(p_f*I.chi)
815              S_fg_I=integrate(S_fg*I.chi)
816              rho_fw_I = self.rho_w.getValue(p_f_I)
817              rho_fg_I = self.rho_g.getValue(p_f_I)
818              A_fw_I= rho_fw_I * self.k_w(1-S_fg_I)/self.mu_w(p_f_I)
819              A_fg_I= rho_fg_I * self.k_g(S_fg_I)/self.mu_g(p_f_I)
820
821                  BHP_limit_I=I.getBHPLimit()
822              BHP_I=I.getBHP()
823
824
825
826              if self.verbose: print "reservoir pressure for well %s = %e gas rho= %e)."%(I.name, p_f_I, rho_fg_I)
827
828              if I.isOpen(self.t+dt):
829            if BHP_limit_I < BHP_I:
830              if I.getPhase() == Well.WATER:
831                  q_I_w =  q_I
832                  q_I_g = Pi_I * A_fg_I * (p_f_I - BHP_I)  /self.rho_g.rho_0
833              else:
834                  q_I_g =  q_I
835                  q_I_w = Pi_I * A_fw_I * (p_f_I - BHP_I)/ self.rho_w.rho_0
836            else:
837                q_I_w = Pi_I * A_fw_I * (p_f_I - BHP_I)/ self.rho_w.rho_0
838                q_I_g = Pi_I * A_fg_I * (p_f_I - BHP_I)  /self.rho_g.rho_0
839              else:
840              q_I_g =0
841              q_I_w =0
842              I.setWaterProduction(q_I_w)
843              I.setGasProduction(q_I_g)
844       return u2       return u2

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