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

revision 3478 by gross, Wed Mar 23 04:02:39 2011 UTC revision 3497 by gross, Wed Apr 6 06:11:56 2011 UTC
# Line 1  Line 1
1  #######################################################  ######################################################
2  #  #
3  # Copyright (c) 2003-2010 by University of Queensland  # Copyright (c) 2003-2010 by University of Queensland
4  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
22
23  from esys.finley import Rectangle, Brick  from esys.escript.linearPDEs import LinearPDE
24  from esys.escript import unitsSI as U  from esys.escript import unitsSI as U
25    from esys.escript.pdetools import Locator
26    from esys.escript import sqrt, log, whereNegative, sup, inf, whereNonPositive, integrate, Function, kronecker, grad, Lsup, clip
27    from esys.weipa import saveVTK
28    import math
29
30  class MaterialProperty(object):  class MaterialProperty(object):
31     """     """
# Line 61  class MaterialPropertyWithDifferential(M Line 65  class MaterialPropertyWithDifferential(M
65        :remark: needs to be overwritten        :remark: needs to be overwritten
66        """        """
67        raise NotImplementedError        raise NotImplementedError
68
69    class Porosity(MaterialPropertyWithDifferential):
70        """
71        defines porosity phi as function of pressure
72
73             phi = phi_p_ref /(1 + X + X**2/2 ) with X= C * (p-p_ref)
74
75        phi_p_ref is claculted from the initial porosity phi_0 at  pressure p_0
76        """
77        def __init__(self, phi_0, p_0, p_ref=1.*U.atm , C=4.0e-5/U.bar):
78          """
79          """
80          self.phi_p_ref=1 # will be overwritten later
81          self.p_ref=p_ref
82          self.C=C
83          # reset  phi_p_ref to get phi(p_0)=phi_0
84          self.phi_p_ref=phi_0/self.getValue(p_0)
85
86        def getValue(self, p):
87          """
88          returns the porosity for given pressure p
89          """
90          X= self.C * ( p - self.p_ref )
91          return  self.phi_p_ref * (1. + X * (1. + X/2 ) )
92
93        def getValueDifferential(self, p):
94          """
95          returns the porosity for given pressure p
96          """
97          X= self.C * ( p - self.p_ref )
98          return  self.phi_p_ref * self.C * (1. + X )
99
100
101    class WaterDensity(MaterialPropertyWithDifferential):
102        """
103        set water density as
104
105              rho = rho_ref (1 + X + X*X/2) with X= C * ( p - p_ref )
106
107        with rho_ref =rho_s/B_ref * gravity
108        """
109        def __init__(self, B_ref=1., p_ref = 1.*U.atm, gravity=1.,  C = 0./U.bar, rho_s= 998.2 * U.kg/U.m**3):
110          self.rho_0 = rho_s * gravity
111          self.rho_ref = self.rho_0/B_ref
112          self.C=C
113          self.p_ref=p_ref
114
115        def getValue(self, p):
116          """
117          returns the density for given pressure p
118          """
119          X= self.C * ( p - self.p_ref )
120          return self.rho_ref * (1+ X * (1+ X/2) )
121
122        def getValueDifferential(self,  p):
123          """
124          """
125          X= self.C * ( p - self.p_ref )
126          return self.rho_ref * self.C * (1+ X)
127
128        def getFormationFactor(self, p):
129          return self.rho_0/self.getValue(p)
130
131
132
133    class WaterViscosity(MaterialProperty):
134      """
135      defines viscosity of water
136
137      mu=mu_ref /(1 + X + X*X/2)
138
139      with X=-C*(p-p_ref)
140      """
141
142      def __init__(self, mu_ref = 0.3*U.cPoise, p_ref = 1.*U.atm , C = 0./U.bar):
143        """
144        set up
145        """
146        self.mu_ref = mu_ref
147        self.p_ref = p_ref
148        self.C=C
149      def getValue(self, p):
150          """
151          returns the viscosity for given pressure p
152          """
153          X= -self.C * ( p - self.p_ref )
154          return self.mu_ref/(1+ X * (1+ X/2) )
155
156    class GasDensity(MaterialPropertyWithDifferential):
157        """
158        set water density as
159
160              rho = gravity * rho_air_s /B(p)
161
162        where B is given by an interpolation table
163        """
164        def __init__(self, p, B, gravity=1., rho_air_s=1.2041*U.kg/U.m**3):
165          self.rho_ref =rho_air_s * gravity
166          self.rho_0 =rho_air_s * gravity
167          self.tab=InterpolationTable(x=p, y=B)
168
169        def getValue(self, p):
170          """
171          returns the density for given pressure p
172          """
173          return self.rho_ref/self.getFormationFactor(p)
174
175        def getValueDifferential(self,  p):
176          """
177          """
178          B    = self.getFormationFactor(p)
179          dBdp = self.getFormationFactorDifferential(p)
180          return  -self.rho_ref * dBdp/(B * B)
181
182        def getFormationFactor(self, p):
183          return self.tab.getValue(p)
184
185  class InterpolationTable(object):      def getFormationFactorDifferential(self, p):
186          return self.tab.getValueDifferential(p)
187
188
189    class InterpolationTable(MaterialPropertyWithDifferential):
190     """     """
191     a simple 1D interpolation table for escript Data object with non-equidistant nodes     a simple 1D interpolation table for escript Data object with non-equidistant nodes
192     """     """
# Line 72  class InterpolationTable(object): Line 196  class InterpolationTable(object):
196        the interpolation argument is below min(x) or larger than max(x). Otherwise        the interpolation argument is below min(x) or larger than max(x). Otherwise
197        the value for x is set to y[0] for        the value for x is set to y[0] for
198        """        """
199          MaterialPropertyWithDifferential.__init__(self)
200        if len(x) != len(y):        if len(x) != len(y):
201       raise ValueError,"length of interpolation nodes and value lists need to be identical."       raise ValueError,"length of interpolation nodes and value lists need to be identical."
202        if len(x) > 0:        if len(x) < 1 :
203       raise ValueError,"length of interpolation nodes a list needs to at least one."       raise ValueError,"length of interpolation nodes a list needs to at least one."
204
205        x_ref=x[0]        x_ref=x[0]
# Line 85  class InterpolationTable(object): Line 210  class InterpolationTable(object):
210        self.__x = x        self.__x = x
211        self.__y = y        self.__y = y
212        self.__obeyBounds = obeyBounds        self.__obeyBounds = obeyBounds
213
214     def getValue(self, x):     def getValue(self, x):
215        """        """
216        returns the interpolated values of x        returns the interpolated values of x
217        """        """
218        x0=self.__x[0]        X=self.__x
219          Y=self.__y
220
221          x0=X[0]
222        m0=whereNegative(x-x0)        m0=whereNegative(x-x0)
223        if self.__obeyBounds:        if self.__obeyBounds:
224       if sup(m0) > 0:       if sup(m0) > 0:
225          raise ValueError,"interpolation argument out of range [%e, %e]"%(x[0],x[-1])          raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
226        out=self.__x[0]        out=self.__x[0]
227        for i in range(1,len(self.__x)):        for i in range(1,len(X)):
228        x1=self.__x[i]        z=(Y[i]-Y[i-1])/(X[i]-X[i-1]) * (x-X[i-1]) + Y[i-1]
229        z=(y[i]-y[i-1])/(x[i]-x[i-1]) * (x-x[i-1]) + y[i-1]        out = out * m0 + z * (1-m0)
230          m0=whereNonPositive(x-X[i])
231
232          if self.__obeyBounds:
233            if inf(m0) < 1 :
234               raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
235          else:
236            out = out * m0 + V[-1] * (1-m0)
237          return out
238
239       def getValueDifferential(self, x):
240          X=self.__x
241          Y=self.__y
242
243          x0=X[0]
244          m0=whereNegative(x-x0)
245          if self.__obeyBounds:
246         if sup(m0) > 0:
247            raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
248          out=0.
249          for i in range(1,len(X)):
250          z=(Y[i]-Y[i-1])/(X[i]-X[i-1])
251        out = out * m0 + z * (1-m0)        out = out * m0 + z * (1-m0)
252        m0=whereNegative(x-x[i])        m0=whereNonPositive(x-X[i])
253
254        if self.__obeyBounds:        if self.__obeyBounds:
255          if inf(m0) < 1:          if inf(m0) < 1:
256             raise ValueError,"interpolation argument out of range [%e, %e]"%(x[0],x[-1])             raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
257        else:        else:
258          out = out * m0 + y[-1] * (1-m0)          out = out * m0
259        return out        return out
260
261
262  class Well(object):  class Well(object):
263     """     """
# Line 117  class Well(object): Line 268  class Well(object):
268     """     """
269     WATER="Water"     WATER="Water"
270     GAS="Gas"     GAS="Gas"
271     def __init__(self,  *args, **kwargs):     def __init__(self, name, domain, Q=0., schedule=[0.], phase="Water", BHP_limit=1.*U.atm, *args, **kwargs):
272         """         """
273         set up well         set up well
274         """         """
275         pass         self.__schedule=schedule
276           self.__phase=phase
277           self.__Q = Q
278           self.__BHP_limit=BHP_limit
279           self.name=name
280           self.domain=domain
281
282     def getGeometry(self):     def getGeometry(self):
283         """         """
284         returns a function representing the geometry of the well.         returns a function representing the geometry of the well.
# Line 144  class Well(object): Line 301  class Well(object):
301        """        """
302        returns the flow rate        returns the flow rate
303        """        """
304        return 0.        return self.__Q
305
306       def getBHPLimit(self):
307          """
308          return bottom-hole pressure limit
309
310          :note: needs to be overwritten
311          """
312          return self.__BHP_limit
313
314     def getBHP(self):     def getBHP(self):
315        """        """
# Line 152  class Well(object): Line 317  class Well(object):
317
318        :note: needs to be overwritten        :note: needs to be overwritten
319        """        """
320        raise NotImplementedError        return self.__BHP
321
322       def setBHP(self, BHP):
323          """
324          sets bottom-hole pressure
325          """
326          self.__BHP= BHP
327
328     def getPhase(self):     def getPhase(self):
329        """        """
330        returns the pahse the well works on        returns the pahse the well works on
331        """        """
332        return self.WATER        return self.__phase
333
334  class PeacemanWell(Well):     def isOpen(self, t):
335         """
336         returns True is the well is open at time t
337         """
338         out = False
339         t0=min(t, self.__schedule[0])
340         i=0
341         while t > t0:
342           if out:
343         out=False
344           else:
345         out=True
346           i+=1
347           if i < len(self.__schedule):
348          t0=self.__schedule[i]
349           else:
350          t0=t
351         return out
352
353       def setWaterProduction(self, v):
354          self.water_production=v
355       def getWaterProduction(self):
356          return self.water_production
357       def setGasProduction(self, v):
358          self.gas_production=v
359       def getGasProduction(self):
360          return self.gas_production
361
362    class VerticalPeacemanWell(Well):
363     """     """
364     defines a well using Peaceman's formula     defines a well using Peaceman's formula
365     """     """
366     pass     def __init__(self,name, domain, schedule = [0.], BHP_limit=1.*U.atm, Q=0, r=10*U.cm, X0=[0.,0.,0.], D=[1.*U.km,1.*U.km, 1.*U.m],
367                 perm=[1.*U.cPoise,1.*U.cPoise, 1.*U.cPoise], phase=Well.WATER, s=0):
368           """
369           set up well
370
371           :param BHP: ottom-hole pressure
372           :param Q: flow rate
374           :param X: location
375           :param D: dimension of well block
376           :param perm: permeabilities
377           :param s: skin factor
378           """
379           Well.__init__(self, name, domain, Q=Q, schedule=schedule, phase=phase,BHP_limit=BHP_limit)
380           r_el=0.28 * sqrt( sqrt(perm[1]/perm[0]) * D[0]**2 +  sqrt(perm[0]/perm[1]) * D[1]**2 )\
381                             / ( (perm[1]/perm[0])**0.25 + (perm[1]/perm[0])**0.25 )
382           self.__PI = 2 * math.pi * D[2] * sqrt(perm[1]*perm[0]) / (log(r_el/r) + s)
383           self.__r = r
384
385           self.__D=D
386           self.r_el=r_el
387           self.X0=X0[:self.domain.getDim()]
388
389           x=Function(domain).getX()
390           self.chi = 1.
391           for i in range(domain.getDim()):
392            self.chi = self.chi * whereNegative(abs(x[i]-X0[i])-D[i]/2)
393
394           self.chi*=1./(D[0]*D[1]*D[2])
395
396
397           #self.chi=0.00000001
398       def getLocation(self):
399           return self.X0
400
401       def getGeometry(self):
402           """
403           returns a function representing the geometry of the well.
404           typically a Gaussian profile around the location of the well.
405
406           :note: needs to be overwritten
407           """
408           return self.chi
409
410       def getProductivityIndex(self):
411           """
412           returns the productivity index of the well.
413           """
414           return self.__PI
415
416
417  class DualPorosity(object):  class DualPorosity(object):
418     """     """
# Line 188  class DualPorosity(object): Line 436  class DualPorosity(object):
436        :type phi_m: `MaterialPropertyWithDifferential` or None        :type phi_m: `MaterialPropertyWithDifferential` or None
437        :param phi: total porosity if None phi=1.        :param phi: total porosity if None phi=1.
438        :type phi: scalar or None        :type phi: scalar or None
439        :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=p_fg-p_fw        :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=S_fg- p_f
440        :type S_fg: `MaterialPropertyWithDifferential`        :type S_fg: `MaterialPropertyWithDifferential`
441        :param S_mg: gas saturation in the coal matrix as function of the capillary pressure p_mc=p_mg-p_mw        :param S_mg: gas saturation in the coal matrix as function of the capillary pressure p_mc=p_mg-p_mw
442        :type S_mg: `MaterialPropertyWithDifferential`        :type S_mg: `MaterialPropertyWithDifferential`
# Line 235  class DualPorosity(object): Line 483  class DualPorosity(object):
483        self.rho_w = rho_w        self.rho_w = rho_w
484        self.rho_g = rho_g            self.rho_g = rho_g
485        self.wells=wells        self.wells=wells
486          self.t =0
487
488          self.__iter_max=1
489          self.__rtol=1.e-4
490          self.verbose=False
491          self.XI=0.5
492       def setIterationControl(self, iter_max=None, rtol=None, verbose=None, xi=None):
493         """
494         sets parameters to control iteration process
495         """
496         if iter_max !=None: self.__iter_max=iter_max
497         if rtol !=None: self.__rtol = rtol
498         if verbose !=None: self.verbose=verbose
499         if xi !=None: self.XI=xi
500
501       def update(self, dt):
502             self.u, self.u_old = self.u.copy(), self.u
503             n=0
504             rerr=1.
505             while n < self.__iter_max and rerr > self.__rtol:
506                u=self.solvePDE(dt)
507            norm_u=Lsup(u)
508            norm_e=Lsup(u-self.u)
509
510            if norm_u > 0:
511                rerr=norm_e/norm_u
512            else:
513                rerr=norm_e
514            if self.verbose: print "iteration step %d: relative change = %e"%(n, rerr)
515            n+=1
516            self.u=u
517             print "iteration completed."
518             self.t+=dt
519
520
521  class PorosityOneHalfModel(DualPorosity):  class PorosityOneHalfModel(DualPorosity):
# Line 247  class PorosityOneHalfModel(DualPorosity) Line 525  class PorosityOneHalfModel(DualPorosity)
525        This corresponds to the coal bed methan model in the ECLIPSE code.        This corresponds to the coal bed methan model in the ECLIPSE code.
526        """        """
527
528        def __init__(self, domain, phi_f=None, phi=None, L_g=None, S_fg=None,        def __init__(self, domain, phi_f=None, phi=None, L_g=None,
529               perm_f_0=None, perm_f_1=None, perm_f_2=None,               perm_f_0=None, perm_f_1=None, perm_f_2=None,
530               k_w =None, k_g=None, mu_w =None, mu_g =None,               k_w =None, k_g=None, mu_w =None, mu_g =None,
531               rho_w =None, rho_g=None,               rho_w =None, rho_g=None, sigma=0, A_mg=0, f_rg=1.,
532                 wells=[]):                 wells=[]):
533       """       """
534       set up       set up
# Line 263  class PorosityOneHalfModel(DualPorosity) Line 541  class PorosityOneHalfModel(DualPorosity)
541       :type phi: scalar or None       :type phi: scalar or None
542       :param L_g: gas adsorption as function of gas pressure       :param L_g: gas adsorption as function of gas pressure
543       :type L_g: `MaterialPropertyWithDifferential`       :type L_g: `MaterialPropertyWithDifferential`
544       :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=p_fg-p_fw       :param S_fg: gas saturation in the fractured rock as function of the capillary pressure p_fc=S_fg- p_f
545       :type S_fg: `MaterialPropertyWithDifferential`       :type S_fg: `MaterialPropertyWithDifferential`
546       :param perm_f_0: permeability in the x_0 direction in the fractured media       :param perm_f_0: permeability in the x_0 direction in the fractured media
547       :type perm_f_0: scalar       :type perm_f_0: scalar
# Line 285  class PorosityOneHalfModel(DualPorosity) Line 563  class PorosityOneHalfModel(DualPorosity)
563       :type rho_g: `MaterialPropertyWithDifferential`       :type rho_g: `MaterialPropertyWithDifferential`
564       :param wells : list of wells       :param wells : list of wells
565       :type wells: list of `Well`       :type wells: list of `Well`
566         :param sigma: shape factor for gas matrix diffusion
567         :param A_mg: diffusion constant for gas adsorption
568         :param f_rg: gas re-adsorption factor
569       """       """
570
571       DualPorosity.__init__(self, domain,       DualPorosity.__init__(self, domain,
572                    phi_f=phi_f, phi_m=None, phi=phi,                    phi_f=phi_f, phi_m=None, phi=phi,
573                    S_fg=S_fg, S_mg=None,                    S_fg=None, S_mg=None,
574                    perm_f_0=perm_f_0, perm_f_1=perm_f_1, perm_f_2=perm_f_2,                    perm_f_0=perm_f_0, perm_f_1=perm_f_1, perm_f_2=perm_f_2,
575                    perm_m_0=None , perm_m_1=None, perm_m_2=None,                    perm_m_0=None , perm_m_1=None, perm_m_2=None,
576                    k_w =k_w, k_g=k_g, mu_w =mu_w, mu_g =mu_g,                    k_w =k_w, k_g=k_g, mu_w =mu_w, mu_g =mu_g,
577                    rho_w =rho_w, rho_g=rho_g,                    rho_w =rho_w, rho_g=rho_g,
578                    wells=wells)                    wells=wells)
579       self.L_g=L_g       self.L_g=L_g
580         self.sigma = sigma
581         self.A_mg = A_mg
582         self.f_rg  = f_rg
583       self.__pde=LinearPDE(self.domain, numEquations=3, numSolutions =3)       self.__pde=LinearPDE(self.domain, numEquations=3, numSolutions =3)
584
585
# Line 305  class PorosityOneHalfModel(DualPorosity) Line 589  class PorosityOneHalfModel(DualPorosity)
589       """       """
590       return self.__pde.getSolverOptions()       return self.__pde.getSolverOptions()
591
592        def solvePDE(self):        def getState(self):
593             return self.u[0], self.u[1],  self.u[2]
p_fw=self.u[0]
p_fg=self.u[1]
p_fc=p_fg-p_fw
C_mg=self.u[3]

p_fw_old=self.u_old[0]
p_fg_old=self.u_old[1]
C_mg_old=self.u_old[3]
594
595          def getOldState(self):
596         return self.u_old[0], self.u_old[1],  self.u_old[2]
597
598       phi_f   =self.phi_f.getValue(p_fg)        def setInitialState(self, p=1.*U.atm, S_fg=0,  C_mg=None):
599       dphi_fdp=self.phi_f.getValueDifferential(p_fg)          """
600            sets the initial state
601
602            :param p: pressure
603            :param S_fg: gas saturation in fractured rock
604            :param C_mg: gas concentration in coal matrix. if not given it is calculated
605                using the  gas adsorption curve.
606            """
607            self.u=self.__pde.createSolution()
608            self.u[0]=p
609            self.u[1]=S_fg
610            if C_mg == None:
611              self.u[2]= self.L_g(p)*self.rho_g.getFormationFactor(p)
612            else:
613              self.u[2]=C_mg
614
615          def solvePDE(self, dt):
616
617       S_fg=  self.S_fg.getValue(p_fc)       C_couple=1
dS_fgdp=   self.S_fg.getValueDifferential(p_fc)
S_fw=  1-S_fg
618
619       rho_fw     = self.rho_w.getValue(p_fw)       p_f, S_fg, C_mg=self.getState()
620       drho_fwdp  = self.rho_w.getValueDifferential(p_fw)       p_f_old, S_fg_old, C_mg_old=self.getOldState()
621       rho_fg = self.rho_g.getValue(p_fg)
622       drho_fgdp  = self.rho_g.getValueDifferential(p_fg)           S_fw=1-S_fg
623
624         if self.verbose:
625              print "p_f range = ",inf(p_f),sup(p_f)
626              print "S_fg range = ",inf(S_fg),sup(S_fg)
627              print "S_fw range = ",inf(S_fw),sup(S_fw)
628              print "C_mg range = ",inf(C_mg),sup(C_mg)
629
630             k_fw=self.k_w(S_fw)
631             if self.verbose: print "k_fw range = ",inf(k_fw),sup(k_fw)
632
633             k_fg=self.k_g(S_fg)
634             if self.verbose: print "k_fg range = ",inf(k_fg),sup(k_fg)
635
636         mu_fw=self.mu_w(p_f)
637             if self.verbose: print "mu_fw range = ",inf(mu_fw),sup(mu_fw)
638
639         mu_fg=self.mu_g(p_f)
640             if self.verbose: print "mu_fg range = ",inf(mu_fg),sup(mu_fg)
641
642
643         phi_f   =self.phi_f.getValue(p_f)
644         dphi_fdp=self.phi_f.getValueDifferential(p_f)
645             if self.verbose: print "phi_f range = ",inf(phi_f),sup(phi_f)," (slope %e,%e)"%(inf(dphi_fdp), sup(dphi_fdp))
646
647       L_g       = self.getValue(p_fg)       rho_fw     = self.rho_w.getValue(p_f)
648       dL_gdp_fg =  self.rho_w.getValueDifferential(p_fg)       drho_fwdp  = self.rho_w.getValueDifferential(p_f)
649             if self.verbose: print "rho_fw range = ",inf(rho_fw),sup(rho_fw)," (slope %e,%e)"%(inf(drho_fwdp), sup(drho_fwdp))
650
651             rho_fg = self.rho_g.getValue(p_f)
652         drho_fgdp  = self.rho_g.getValueDifferential(p_f)
653             if self.verbose: print "rho_fg range = ",inf(rho_fg),sup(rho_fg)," (slope %e,%e)"%(inf(drho_fgdp), sup(drho_fgdp))
654
655       A_fw = rho_fw * k_fw/mu_fw       L_g_0       = self.L_g(p_f)
656       A_fg = rho_fg * k_fg/mu_fg       dL_g_0dp    = self.L_g.getValueDifferential(p_f)
657         FF_g=self.rho_g.getFormationFactor(p_f)
658         dFF_gdp=self.rho_g.getFormationFactorDifferential(p_f)
659
660       m = whereNegative(L_g - C_mg)         L_g = L_g_0 * FF_g
661       B = self.sigma * self.a_mg * (m + (1-m) * self.f_rg * S_fg )       dL_gdp =  0 * (dL_g_0dp * FF_g + L_g_0 * dFF_gdp)
662
663             if self.verbose:
664              print "L_g_0 range = ",inf(L_g_0),sup(L_g_0)
665              print "L_g range = ",inf(L_g),sup(L_g)," (slope %e,%e)"%(inf(dL_gdp), sup(dL_gdp))
666
667
668       E_fww= phi_f * ( rho_fw * dS_fgdp + S_fw * drho_fwdp )       A_fw = rho_fw * k_fw/mu_fw
669       E_fwg= rho_fw * ( S_fw * dphi_fdp - phi_f * dS_fgdp )       A_fg = rho_fg * k_fg/mu_fg
E_fgw= - phi_f * rho_fg * dS_fgdp
E_fgg= S_fg * rho_fg  * dphi_fdp + phi_f * rho_fg * dS_fgdp + phi_f  * S_fg * drho_fgdp

670
671       Q_w=0.       m = whereNegative(L_g - C_mg)
672       Q_g=0.       B = self.sigma * self.A_mg * (m + (1-m) * self.f_rg * S_fg )
673       R_w=0.       print "XXXX B =",B
674       R_g=0.       print "XXXX self.sigma  =",self.sigma
675         print "XXXX self.A_mg  =",self.A_mg
676
677         E_fpp= S_fw * (  rho_fw * dphi_fdp + phi_f  * drho_fwdp )
678         E_fps=  -  phi_f * rho_fw
679         E_fsp= S_fg *( rho_fg * dphi_fdp - phi_f *  drho_fgdp )
680         E_fss= phi_f * rho_fg
681
682         F_fw=0.
683         F_fg=0.
684         D_fw=0.
685         D_fg=0.
686
687       for I in self.wells:       for I in self.wells:
688          chi_I= I.getGeometry()            chi_I= I.getGeometry()
689          Pi_I = I.getProductivityIndex()            loc=Locator(Function(self.domain),I.getLocation())
690          P_I = I.getBHP()            Pi_I = I.getProductivityIndex()
691          Q_I = I.getFlowRate()            A_fw_I= loc(A_fw)
692                      A_fg_I= loc(A_fg)
693          if self.getPhase == Well.WATER:            BHP_limit_I=I.getBHPLimit()
694             Q_w = Q_w + Q_I        * chi_I
695             Q_g = Q_g + Pi_I * P_I * chi_I            if I.isOpen(self.t+dt):
696             R_g = R_g + Pi_I       * chi_I          if self.verbose: print "well %s is open."%I.name
697          else:          if I.getPhase() == Well.WATER:
698             Q_g = Q_g + Q_I        * chi_I              q_I = self.rho_w.rho_0 * I.getFlowRate()
699             Q_w = Q_w + Pi_I * P_I * chi_I              p_I_ref=loc(p_f)-q_I/(A_fw_I * Pi_I)
700             R_w = R_w + Pi_I       * chi_I          else:
701                          q_I = self.rho_g.rho_0 * I.getFlowRate()
702       F_fw_Y = A_fw * Q_w              p_I_ref=loc(p_f)-q_I/(A_fg_I * Pi_I)
703       F_fg_Y = A_fg * Q_g
704       D_fgg = - A_fg * R_g          print "ZZZ =",loc(p_f), q_I, self.rho_w.rho_0, I.getFlowRate(), A_fw_I, Pi_I, q_I/(A_fw_I * Pi_I)
705       D_fww = - A_fw * R_w
706            if BHP_limit_I > p_I_ref:
707              BHP_I=BHP_limit_I
708              D_fw = D_fw + Pi_I * A_fw_I *              chi_I
709              F_fw = F_fw - Pi_I * A_fw_I * BHP_limit_I *chi_I
710              D_fg = D_fg + Pi_I * A_fg_I *              chi_I
711              F_fg = F_fg - Pi_I * A_fg_I * BHP_limit_I *chi_I
712            else:
713              if I.getPhase() == Well.WATER:
714                  F_fw = F_fw +  q_I  * chi_I
715                  F_fg = F_fg +  A_fg_I/A_fw_I *  q_I *chi_I
716              else:
717                  F_fg = F_fg +  q_I  * chi_I
718                  F_fw = F_fw +  A_fw_I/A_fg_I *  q_I *chi_I
719              BHP_I=p_I_ref
720              else:
721              if self.verbose: print "well %s is closed."%I.name
722              BHP_I=loc(p_f)
723
724              if self.verbose: print "well %s:BHP = %e"%(I.name, BHP_I)
725              I.setBHP(BHP_I)
726
727         F_fp_Y = - F_fw
728         F_fs_Y = - F_fg
729         D_fpp =  D_fw
730         D_fsp =  D_fg
731
732
733       if self.domain.getDim() == 3:       if self.domain.getDim() == 3:
734          F_fw_X = A_fw * self.g * rho_fw * self.perm_f[2] * kronecker(self.domain)[2]          F_fp_X = ( A_fw * self.g * rho_fw * self.perm_f[2] ) * kronecker(self.domain)[2]
735          F_fg_X = A_fg * self.g * rho_fg * self.perm_f[2] * kronecker(self.domain)[2]          F_fs_X = ( A_fg * self.g * rho_fg * self.perm_f[2] ) * kronecker(self.domain)[2]
736       else:       else:
737          F_fw_X = 0 * kronecker(self.domain)[1]          F_fp_X = 0 * kronecker(self.domain)[1]
738          F_fg_X = 0 * kronecker(self.domain)[1]          F_fs_X = 0 * kronecker(self.domain)[1]
739
740       F_mg = B * ( L_g - dL_gdp_fg * p_fg )       #F_mg_Y = B * ( L_g - dL_gdp * p_f )
741             F_mg_Y = B * L_g
742
743
744       D=self.__pde.getNewCoefficient("D")       D=self.__pde.createCoefficient("D")
745       A=self.__pde.getNewCoefficient("A")       A=self.__pde.createCoefficient("A")
746       Y=self.__pde.getNewCoefficient("Y")       Y=self.__pde.createCoefficient("Y")
747       X=self.__pde.getNewCoefficient("X")       X=self.__pde.createCoefficient("X")
748
749       dtXI = dt*self.XI       dtXI = dt*self.XI
750       dtcXI = dt*(1-self.XI)       dtcXI = dt*(1-self.XI)
751
752       D[0,0]=E_fww + dtXI * D_fww       D[0,0]=E_fpp + dtXI * D_fpp
753       D[0,1]=E_fwg       D[0,1]=E_fps
754       D[1,0]=E_fgw       D[1,0]=E_fsp + dtXI * D_fsp
755       D[1,1]=E_fgg + dtXI * D_fgg       D[1,1]=E_fss
756       D[1,2]=rho_fg       D[1,2]=rho_fg * C_couple
757       D[2,1]= dtXI * B * dL_gdp_fg       #D[2,0]= - dtXI * B * dL_gdp
758       D[2,2]= 1 + dtXI * B       D[2,2]= 1 + dtXI * B
759
H_fw = dtcXI * A_fw * grad(p_fw_old)
H_fg = dtcXI * A_fg * grad(p_fg_old)
760
761       for i in range(self.domain.getDim()):       for i in range(self.domain.getDim()):
762          A[0,i,0,i] = dtXI * A_fw * self.perm_f[i]          A[0,i,0,i] = dtXI * A_fw * self.perm_f[i]
763          A[1,i,1,i] = dtXI * A_fg * self.perm_f[i]          A[1,i,1,i] = dtXI * A_fg * self.perm_f[i]
764
765          X[0,i]= dt * F_fw_X - H_fw * self.perm_f[i]       g= grad(p_f_old) * dtcXI * self.perm_f[0:self.domain.getDim()]
766          X[1,i]= dt * F_fg_X - H_fg * self.perm_f[i]           X[0,:]=  - A_fw * g  + dt * F_fp_X
767                 X[1,:]=  - A_fg * g  + dt * F_fs_X
768       Y[0] = E_fww * p_fw_old + E_fwg * p_fg_old +                     dt * F_fw_Y - dtcXI * D_fww * p_fw_old
769       Y[1] = E_fgw * p_fw_old + E_fgg * p_fg_old + rho_fg * C_mg_old + dt * F_fg_Y - dtcXI * D_fgg * p_fg_old       Y[0] = E_fpp *  p_f_old + E_fps * S_fg_old +                                dt * F_fp_Y - dtcXI * D_fpp * p_f_old
770       Y[2] = C_mg_old                                                + dt * F_mg_Y - dtcXI * ( dL_gdp_fg * p_fg_old + C_mg_old)       Y[1] = E_fsp *  p_f_old + E_fss * S_fg_old + C_couple * rho_fg * C_mg_old + dt * F_fs_Y - dtcXI * D_fsp * p_f_old
771         Y[2] = C_mg_old                                                           + dt * F_mg_Y - dtcXI * B * C_mg_old  # - dL_gdp * p_f_old)
772
773         print "HHH Y[0] =", Y[0]
774         print "HHH A_fw = ",A_fw
775         print "HHH A_fg= ",A_fg
776         print "HHH F_fg = ",F_fg
777         print "HHH F_fw = ",F_fw
778         print "HHH perm_f = ",self.perm_f
779         print "HHH k = ",(A_fw*self.perm_f[0])/D[0,0]
780         print "HHH k = ",(A_fw*self.perm_f[1])/D[0,0]
781         print "HHH X[0,:]= ",X[0,:]
782         print "HHH D[0,0] = ",D[0,0]
783         print "HHH D[0,1] = ",D[0,1]
784         print "HHH D[0,2] = ",D[0,2]
785
786         print "hhh B= ",B
787
788
789
790
791         #print "HHH Y[1] = ",Y[1]
792         #print "HHH X[1,:] = ",X[1,:]
793         #print "HHH D[1,0] = ",D[1,0]
794         #print "HHH D[1,1] = ",D[1,1]
795         #print "HHH D[1,2] = ",D[1,2]
796         #print "HHH A_fg =",  A_fg
797
798
799
800
801         #1/0
802       self.__pde.setValue(A=A, D=D, X=X, Y=Y)       self.__pde.setValue(A=A, D=D, X=X, Y=Y)
803
804       u2 = self.__pde.getSolution()       u2 = self.__pde.getSolution()
805         # this is deals with round-off errors to maintain physical meaningful values
806         # we should do this in a better way to detect values that are totally wrong
807         u2[1]=clip(u2[1],minval=0, maxval=1)
808         u2[2]=clip(u2[2],minval=0)
809
810         # update well production
811         p_f=u2[0]
812         for I in self.wells:
813              loc=Locator(Function(self.domain),I.getLocation())
814              Pi_I = I.getProductivityIndex()
815              q_I = I.getFlowRate()
816              A_fw_I= loc(A_fw)
817              A_fg_I= loc(A_fg)
818              p_f_I=loc(A_fg)
819              BHP_limit_I=I.getBHPLimit()
820              BHP=I.getBHP()
821
822
823              rho_fw_I = loc(rho_fw)
824              rho_fg_I = loc(rho_fg)
825              A_fg = rho_fg * k_fg/mu_fg
826
827              if I.isOpen(self.t+dt):
828            if BHP_limit_I < BHP:
829              if I.getPhase() == Well.WATER:
830                  q_I_w =  q_I
831                  q_I_g =  A_fg_I/A_fw_I *  q_I
832              else:
833                  q_I_g =  q_I
834                  q_I_w =  A_fw_I/A_fg_I *  q_I
835            else:
836                q_I_w = Pi_I * A_fw_I/rho_fw_I * (p_f_I - BHP_I)
837                q_I_g = Pi_I * A_fg_I/rho_fg_I * (p_f_I - BHP_I)
838              else:
839              q_I_g =0
840              q_I_w =0
841              I.setWaterProduction(q_I_w * rho_fw_I/self.rho_w.rho_0)
842              I.setGasProduction(q_I_g * rho_fg_I/self.rho_g.rho_0 )
843       return u2       return u2

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