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

revision 3478 by gross, Wed Mar 23 04:02:39 2011 UTC revision 3484 by gross, Thu Mar 24 05:42:07 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
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 InterpolationTable(object):  class WaterDensity(MaterialPropertyWithDifferential):
100        """
101        set water density as
102
103              rho = rho_ref (1 + X + X*X/2) with X= C * ( p - p_ref )
104
105        with rho_ref =rho_s/B_ref * gravity
106        """
107        def __init__(self, B_ref=1., p_ref = 1.*U.atm, gravity=1.,  C = 0./U.bar, rho_s= 998.2 * U.kg/U.m**3):
108          self.rho_ref = rho_s/B_ref * gravity
109          self.C=C
110          self.p_ref=p_ref
111
112        def getValue(self, p):
113          """
114          returns the density for given pressure p
115          """
116          X= self.C * ( p - self.p_ref )
117          return self.rho_ref * (1+ X * (1+ X/2) )
118
119        def getValueDifferential(self,  p):
120          """
121          """
122          X= self.C * ( p - self.p_ref )
123          return self.rho_ref * self.C * (1+ X)
124
125    class WaterViscosity(MaterialProperty):
126      """
127      defines viscosity of water
128
129      mu=mu_ref /(1 + X + X*X/2)
130
131      with X=-C*(p-p_ref)
132      """
133
134      def __init__(self, mu_ref = 0.3*U.cPoise, p_ref = 1.*U.atm , C = 0./U.bar):
135        """
136        set up
137        """
138        self.mu_ref = mu_ref
139        self.p_ref = p_ref
140        self.C=C
141      def getValue(self, p):
142          """
143          returns the viscosity for given pressure p
144          """
145          X= -self.C * ( p - self.p_ref )
146          return self.mu_ref/(1+ X * (1+ X/2) )
147
148    class GasDensity(MaterialPropertyWithDifferential):
149        """
150        set water density as
151
152              rho = gravity * rho_air_s /B(p)
153
154        where B is given by an interpolation table
155        """
156        def __init__(self, p, B, gravity=1., rho_air_s=1.2041*U.kg/U.m**3):
157          self.rho_ref =rho_air_s * gravity
158          self.tab=InterpolationTable(x=p, y=B)
159
160        def getValue(self, p):
161          """
162          returns the density for given pressure p
163          """
164          return self.rho_ref/self.tab.getValue(p)
165
166        def getValueDifferential(self,  p):
167          """
168          """
169          B    = self.tab.getValue(p)
170          dBdp = self.tab.getValueDifferential(p)
171          return  -self.rho_ref * dBdp/(B * B)
172
173    class InterpolationTable(MaterialPropertyWithDifferential):
174     """     """
175     a simple 1D interpolation table for escript Data object with non-equidistant nodes     a simple 1D interpolation table for escript Data object with non-equidistant nodes
176     """     """
# Line 72  class InterpolationTable(object): Line 180  class InterpolationTable(object):
180        the interpolation argument is below min(x) or larger than max(x). Otherwise        the interpolation argument is below min(x) or larger than max(x). Otherwise
181        the value for x is set to y[0] for        the value for x is set to y[0] for
182        """        """
183          MaterialPropertyWithDifferential.__init__(self)
184        if len(x) != len(y):        if len(x) != len(y):
185       raise ValueError,"length of interpolation nodes and value lists need to be identical."       raise ValueError,"length of interpolation nodes and value lists need to be identical."
186        if len(x) > 0:        if len(x) < 1 :
187       raise ValueError,"length of interpolation nodes a list needs to at least one."       raise ValueError,"length of interpolation nodes a list needs to at least one."
188
189        x_ref=x[0]        x_ref=x[0]
# Line 85  class InterpolationTable(object): Line 194  class InterpolationTable(object):
194        self.__x = x        self.__x = x
195        self.__y = y        self.__y = y
196        self.__obeyBounds = obeyBounds        self.__obeyBounds = obeyBounds
197
198     def getValue(self, x):     def getValue(self, x):
199        """        """
200        returns the interpolated values of x        returns the interpolated values of x
# Line 108  class InterpolationTable(object): Line 218  class InterpolationTable(object):
218          out = out * m0 + y[-1] * (1-m0)          out = out * m0 + y[-1] * (1-m0)
219        return out        return out
220
221       def getValueDifferential(self, x):
222          x0=self.__x[0]
223          m0=whereNegative(x-x0)
224          if self.__obeyBounds:
225         if sup(m0) > 0:
226            raise ValueError,"interpolation argument out of range [%e, %e]"%(x[0],x[-1])
227          out=0.
228          for i in range(1,len(self.__x)):
229          x1=self.__x[i]
230          z=(y[i]-y[i-1])/(x[i]-x[i-1])
231          out = out * m0 + z * (1-m0)
232          m0=whereNegative(x-x[i])
233
234          if self.__obeyBounds:
235            if inf(m0) < 1:
236               raise ValueError,"interpolation argument out of range [%e, %e]"%(x[0],x[-1])
237          else:
238            out = out * m0
239          return out
240
241
242  class Well(object):  class Well(object):
243     """     """
244     generic well     generic well
# Line 117  class Well(object): Line 248  class Well(object):
248     """     """
249     WATER="Water"     WATER="Water"
250     GAS="Gas"     GAS="Gas"
251     def __init__(self,  *args, **kwargs):     def __init__(self,  Q=0., schedule=[0.], phase="Water", *args, **kwargs):
252         """         """
253         set up well         set up well
254         """         """
255         pass         self.__schedule=schedule
256           self.__phase=phase
257           self.__Q = Q
258
259     def getGeometry(self):     def getGeometry(self):
260         """         """
261         returns a function representing the geometry of the well.         returns a function representing the geometry of the well.
# Line 140  class Well(object): Line 274  class Well(object):
274         """         """
275         raise NotImplementedError         raise NotImplementedError
276
277     def getFlowRate(self):     def getFlowRate(self, t):
278        """        """
279        returns the flow rate        returns the flow rate
280        """        """
281        return 0.        if self.isOpen(t):
282              return self.__Q
283          else:
284        return 0.
285
286     def getBHP(self):     def getBHP(self):
287        """        """
288        return bottom-hole pressure        return bottom-hole pressure
# Line 158  class Well(object): Line 295  class Well(object):
295        """        """
296        returns the pahse the well works on        returns the pahse the well works on
297        """        """
298        return self.WATER        return self.__phase
299
300       def isOpen(self, t):
301         """
302         returns True is the well is open at time t
303         """
304         out = False
305         t0=min(t, self.__schedule[0])
306         i=0
307         while t < t0:
308           if out:
309         out=False
310           else:
311         out=True
312           i+=1
313           if i < len(self.__schedule[i]):
314          t0=self.__schedule[i]
315           else:
316          t0=t
317         return out
318
319  class PeacemanWell(Well):  class VerticalPeacemanWell(Well):
320     """     """
321     defines a well using Peaceman's formula     defines a well using Peaceman's formula
322     """     """
323     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],
324                 perm=[1.*U.cPoise,1.*U.cPoise, 1.*U.cPoise], phase=Well.WATER, s=0):
325           """
326           set up well
327
328           :param BHP: ottom-hole pressure
329           :param Q: flow rate
331           :param X: location
332           :param D: dimension of well block
333           :param perm: permeabilities
334           :param s: skin factor
335           """
336           Well.__init__(self, Q=Q, schedule=schedule, phase=phase,)
337           r_el=0.28 * ( (perm[1]/perm[0])**0.5 * D[0]**2 +  (perm[0]/perm[1])**0.5 * D[1]**2 ) \
338                             / ( (perm[1]/perm[0])**0.25 + (perm[1]/perm[0])**0.25 )
339           self.__PI = 2 * math.pi * D[2] * sqrt(perm[1]*perm[0]) / (log(r_el/r) + s)
340           self.__r = r
341
342           self.__BHP = BHP
343           self.__D=D
344           self.r_el=r_el
345           self.X0=X0
346
347       def getGeometry(self):
348           """
349           returns a function representing the geometry of the well.
350           typically a Gaussian profile around the location of the well.
351
352           :note: needs to be overwritten
353           """
354           raise NotImplementedError
355
356       def getProductivityIndex(self,t):
357           """
358           returns the productivity index of the well.
359           typically a Gaussian profile around the location of the well.
360           """
361           return self.__PI
362
363       def getBHP(self):
364          """
365          return bottom-hole pressure
366          """
367          return self.__BHP
368
369  class DualPorosity(object):  class DualPorosity(object):
370     """     """
# Line 188  class DualPorosity(object): Line 388  class DualPorosity(object):
388        :type phi_m: `MaterialPropertyWithDifferential` or None        :type phi_m: `MaterialPropertyWithDifferential` or None
389        :param phi: total porosity if None phi=1.        :param phi: total porosity if None phi=1.
390        :type phi: scalar or None        :type phi: scalar or None
391        :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
392        :type S_fg: `MaterialPropertyWithDifferential`        :type S_fg: `MaterialPropertyWithDifferential`
393        :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
394        :type S_mg: `MaterialPropertyWithDifferential`        :type S_mg: `MaterialPropertyWithDifferential`
# Line 247  class PorosityOneHalfModel(DualPorosity) Line 447  class PorosityOneHalfModel(DualPorosity)
447        This corresponds to the coal bed methan model in the ECLIPSE code.        This corresponds to the coal bed methan model in the ECLIPSE code.
448        """        """
449
450        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,
451               perm_f_0=None, perm_f_1=None, perm_f_2=None,               perm_f_0=None, perm_f_1=None, perm_f_2=None,
452               k_w =None, k_g=None, mu_w =None, mu_g =None,               k_w =None, k_g=None, mu_w =None, mu_g =None,
453               rho_w =None, rho_g=None,               rho_w =None, rho_g=None,
# Line 263  class PorosityOneHalfModel(DualPorosity) Line 463  class PorosityOneHalfModel(DualPorosity)
463       :type phi: scalar or None       :type phi: scalar or None
464       :param L_g: gas adsorption as function of gas pressure       :param L_g: gas adsorption as function of gas pressure
465       :type L_g: `MaterialPropertyWithDifferential`       :type L_g: `MaterialPropertyWithDifferential`
466       :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
467       :type S_fg: `MaterialPropertyWithDifferential`       :type S_fg: `MaterialPropertyWithDifferential`
468       :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
469       :type perm_f_0: scalar       :type perm_f_0: scalar
# Line 289  class PorosityOneHalfModel(DualPorosity) Line 489  class PorosityOneHalfModel(DualPorosity)
489
490       DualPorosity.__init__(self, domain,       DualPorosity.__init__(self, domain,
491                    phi_f=phi_f, phi_m=None, phi=phi,                    phi_f=phi_f, phi_m=None, phi=phi,
492                    S_fg=S_fg, S_mg=None,                    S_fg=None, S_mg=None,
493                    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,
494                    perm_m_0=None , perm_m_1=None, perm_m_2=None,                    perm_m_0=None , perm_m_1=None, perm_m_2=None,
495                    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 505  class PorosityOneHalfModel(DualPorosity)
505       """       """
506       return self.__pde.getSolverOptions()       return self.__pde.getSolverOptions()
507
508          def setInitialState(p=1.*U.atm, S_fg=0,  C_mg=None):
509            """
510            sets the initial state
511
512            :param p: pressure
513            :param S_fg: gas saturation in fractured rock
514            :param C_mg: gas concentration in coal matrix. if not given it is calculated
515                using the  gas adsorption curve.
516            """
517
518            self.u=self.__pde.getNewCoefficient("u")
519            u[0]=p
520            u[1]=S_fg
521            if C_mg == None:
522              u[2]= self.L_g(p)
523            else:
524              u[2]=C_mg
525
526        def solvePDE(self):        def solvePDE(self):
527
528       p_fw=self.u[0]       p_f=self.u[0]
529       p_fg=self.u[1]       S_fg=self.u[1]
p_fc=p_fg-p_fw
530       C_mg=self.u[3]       C_mg=self.u[3]
531
532
533       p_fw_old=self.u_old[0]       p_f_old=self.u_old[0]
534       p_fg_old=self.u_old[1]       S_fg_old=self.u_old[1]
535       C_mg_old=self.u_old[3]       C_mg_old=self.u_old[3]
536
537
538       phi_f   =self.phi_f.getValue(p_fg)       phi_f   =self.phi_f.getValue(S_fg)
539       dphi_fdp=self.phi_f.getValueDifferential(p_fg)       dphi_fdp=self.phi_f.getValueDifferential(S_fg)
540
S_fg=  self.S_fg.getValue(p_fc)
dS_fgdp=   self.S_fg.getValueDifferential(p_fc)
541       S_fw=  1-S_fg       S_fw=  1-S_fg
542
543       rho_fw     = self.rho_w.getValue(p_fw)       rho_fw     = self.rho_w.getValue(p_f)
544       drho_fwdp  = self.rho_w.getValueDifferential(p_fw)       drho_fwdp  = self.rho_w.getValueDifferential( p_f)
545       rho_fg = self.rho_g.getValue(p_fg)       rho_fg = self.rho_g.getValue(p_f)
546       drho_fgdp  = self.rho_g.getValueDifferential(p_fg)       drho_fgdp  = self.rho_g.getValueDifferential(p_f)
547
548       L_g       = self.getValue(p_fg)       L_g       = self.getValue(p_f)
549       dL_gdp_fg =  self.rho_w.getValueDifferential(p_fg)       dL_gdp =  self.rho_w.getValueDifferential(p_f)
550
551       A_fw = rho_fw * k_fw/mu_fw       A_fw = rho_fw * k_fw/mu_fw
552       A_fg = rho_fg * k_fg/mu_fg       A_fg = rho_fg * k_fg/mu_fg
# Line 340  class PorosityOneHalfModel(DualPorosity) Line 555  class PorosityOneHalfModel(DualPorosity)
555       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 )
556
557
558       E_fww= phi_f * ( rho_fw * dS_fgdp + S_fw * drho_fwdp )       E_fpp= S_fw * (  rho_fw * dphi_fdp + phi_f  * drho_fwdp )
559       E_fwg= rho_fw * ( S_fw * dphi_fdp - phi_f * dS_fgdp )       E_fps=  -  phi_f * rho_fw
560       E_fgw= - phi_f * rho_fg * dS_fgdp       E_fsp= S_fg *( rho_fg * dphi_fdp - phi_f *  drho_fgdp )
561       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
562
563
564
# Line 367  class PorosityOneHalfModel(DualPorosity) Line 582  class PorosityOneHalfModel(DualPorosity)
582             Q_w = Q_w + Pi_I * P_I * chi_I             Q_w = Q_w + Pi_I * P_I * chi_I
583             R_w = R_w + Pi_I       * chi_I             R_w = R_w + Pi_I       * chi_I
584
585       F_fw_Y = A_fw * Q_w       F_fp_Y = A_fw * Q_w
586       F_fg_Y = A_fg * Q_g       F_fs_Y = A_fg * Q_g
587       D_fgg = - A_fg * R_g       D_fss = - A_fg * R_g
588       D_fww = - A_fw * R_w       D_fpp = - A_fw * R_w
589
590       if self.domain.getDim() == 3:       if self.domain.getDim() == 3:
591          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]
592          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]
593       else:       else:
594          F_fw_X = 0 * kronecker(self.domain)[1]          F_fp_X = 0 * kronecker(self.domain)[1]
595          F_fg_X = 0 * kronecker(self.domain)[1]          F_fs_X = 0 * kronecker(self.domain)[1]
596
597       F_mg = B * ( L_g - dL_gdp_fg * p_fg )       F_mg = B * ( L_g - dL_gdp * S_fg )
598
599
600       D=self.__pde.getNewCoefficient("D")       D=self.__pde.getNewCoefficient("D")
# Line 390  class PorosityOneHalfModel(DualPorosity) Line 605  class PorosityOneHalfModel(DualPorosity)
605       dtXI = dt*self.XI       dtXI = dt*self.XI
606       dtcXI = dt*(1-self.XI)       dtcXI = dt*(1-self.XI)
607
608       D[0,0]=E_fww + dtXI * D_fww       D[0,0]=E_fpp + dtXI * D_fpp
609       D[0,1]=E_fwg       D[0,1]=E_fps
610       D[1,0]=E_fgw       D[1,0]=E_fsp
611       D[1,1]=E_fgg + dtXI * D_fgg       D[1,1]=E_fss + dtXI * D_fss
612       D[1,2]=rho_fg       D[1,2]=rho_fg
613       D[2,1]= dtXI * B * dL_gdp_fg       D[2,1]= dtXI * B * dL_gdp
614       D[2,2]= 1 + dtXI * B       D[2,2]= 1 + dtXI * B
615
616       H_fw = dtcXI * A_fw * grad(p_fw_old)       H_fw = dtcXI * A_fw * grad( p_f_old)
617       H_fg = dtcXI * A_fg * grad(p_fg_old)       H_fg = dtcXI * A_fg * grad(S_fg_old)
618
619       for i in range(self.domain.getDim()):       for i in range(self.domain.getDim()):
620          A[0,i,0,i] = dtXI * A_fw * self.perm_f[i]          A[0,i,0,i] = dtXI * A_fw * self.perm_f[i]
621          A[1,i,1,i] = dtXI * A_fg * self.perm_f[i]          A[1,i,1,i] = dtXI * A_fg * self.perm_f[i]
622
623          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]
624          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]
625
626       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
627       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
628       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)
629
630       self.__pde.setValue(A=A, D=D, X=X, Y=Y)       self.__pde.setValue(A=A, D=D, X=X, Y=Y)
631

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