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

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

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