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

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)
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
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