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

revision 3478 by gross, Wed Mar 23 04:02:39 2011 UTC revision 3527 by gross, Tue May 24 06:59:20 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)
# Line 20  __license__="""Licensed under the Open S Line 20  __license__="""Licensed under the Open S
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 *
27    from esys.weipa import saveVTK
28    import math
29
30    USE_NODAL_WELL = False or True
31
32  class MaterialProperty(object):  class MaterialProperty(object):
33     """     """
# Line 61  class MaterialPropertyWithDifferential(M Line 67  class MaterialPropertyWithDifferential(M
67        :remark: needs to be overwritten        :remark: needs to be overwritten
68        """        """
69        raise NotImplementedError        raise NotImplementedError
70
71    class Porosity(MaterialPropertyWithDifferential):
72        """
73        defines porosity phi as function of pressure
74
75             phi = phi_p_ref /(1 + X + X**2/2 ) with X= C * (p-p_ref)
76
77        phi_p_ref is claculted from the initial porosity phi_0 at  pressure p_0
78        """
79        def __init__(self, phi_0, p_0, p_ref=1.*U.atm , C=4.0e-5/U.bar):
80          """
81          """
82          self.phi_p_ref=1 # will be overwritten later
83          self.p_ref=p_ref
84          self.C=C
85          # reset  phi_p_ref to get phi(p_0)=phi_0
86          self.phi_p_ref=phi_0/self.getValue(p_0)
87
88        def getValue(self, p):
89          """
90          returns the porosity for given pressure p
91          """
92          X= self.C * ( p - self.p_ref )
93          return  self.phi_p_ref * (1. + X * (1. + X/2 ) )
94
95        def getValueDifferential(self, p):
96          """
97          returns the porosity for given pressure p
98          """
99          X= self.C * ( p - self.p_ref )
100          return  self.phi_p_ref * self.C * (1. + X )
101
102
103    class WaterDensity(MaterialPropertyWithDifferential):
104        """
105        set water density as
106
107              rho = rho_surf (1 + X + X*X/2) with X= C * ( p - p_ref )
108
109        with rho_surf =rho_s/B_ref * gravity
110        """
111        def __init__(self, B_ref=1., p_ref = 1.*U.atm, gravity=1.,  C = 0./U.bar, rho_s= 999.014 * U.kg/U.m**3):
112          self.rho_surf = rho_s * gravity
113          self.__rho_0 = self.rho_surf/B_ref
114          self.C=C
115          self.p_ref=p_ref
116
117        def getValue(self, p):
118          """
119          returns the density for given pressure p
120          """
121          X= self.C * ( p - self.p_ref )
122          return self.__rho_0 * (1+ X * (1+ X/2) )
123
124        def getValueDifferential(self,  p):
125          """
126          """
127          X= self.C * ( p - self.p_ref )
128          return self.__rho_0 * self.C * (1+ X)
129
130        def getFormationFactor(self, p):
131          return self.rho_surf/self.getValue(p)
132
133
134
135    class WaterViscosity(MaterialProperty):
136      """
137      defines viscosity of water
138
139      mu=mu_ref /(1 + X + X*X/2)
140
141      with X=-C*(p-p_ref)
142      """
143
144      def __init__(self, mu_ref = 0.3*U.cPoise, p_ref = 1.*U.atm , C = 0./U.bar):
145        """
146        set up
147        """
148        self.mu_ref = mu_ref
149        self.p_ref = p_ref
150        self.C=C
151      def getValue(self, p):
152          """
153          returns the viscosity for given pressure p
154          """
155          X= -self.C * ( p - self.p_ref )
156          return self.mu_ref/(1+ X * (1+ X/2) )
157
158    class GasDensity(MaterialPropertyWithDifferential):
159        """
160        set water density as
161
162              rho = gravity * rho_air_s /B(p)
163
164        where B is given by an interpolation table
165        """
166        def __init__(self, p, B, gravity=1., rho_air_s=1.2041*U.kg/U.m**3):
167          self.rho_surf =rho_air_s * gravity
168          self.tab=InterpolationTable(x=p, y=B)
169
170        def getValue(self, p):
171          """
172          returns the density for given pressure p
173          """
174          return self.rho_surf/self.getFormationFactor(p)
175
176        def getValueDifferential(self,  p):
177          """
178          """
179          B    = self.getFormationFactor(p)
180          dBdp = self.getFormationFactorDifferential(p)
181          return  -self.rho_surf * dBdp/(B * B)
182
183  class InterpolationTable(object):      def getFormationFactor(self, p):
184          return self.tab.getValue(p)
185
186        def getFormationFactorDifferential(self, p):
187          return self.tab.getValueDifferential(p)
188
189
190    class InterpolationTable(MaterialPropertyWithDifferential):
191     """     """
192     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
193     """     """
# Line 72  class InterpolationTable(object): Line 197  class InterpolationTable(object):
197        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
198        the value for x is set to y[0] for        the value for x is set to y[0] for
199        """        """
200          MaterialPropertyWithDifferential.__init__(self)
201        if len(x) != len(y):        if len(x) != len(y):
202       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."
203        if len(x) > 0:        if len(x) < 1 :
204       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."
205
206        x_ref=x[0]        x_ref=x[0]
# Line 85  class InterpolationTable(object): Line 211  class InterpolationTable(object):
211        self.__x = x        self.__x = x
212        self.__y = y        self.__y = y
213        self.__obeyBounds = obeyBounds        self.__obeyBounds = obeyBounds
214
215     def getValue(self, x):     def getValue(self, x):
216        """        """
217        returns the interpolated values of x        returns the interpolated values of x
218        """        """
219        x0=self.__x[0]        X=self.__x
220          Y=self.__y
221
222          x0=X[0]
223        m0=whereNegative(x-x0)        m0=whereNegative(x-x0)
224        if self.__obeyBounds:        if self.__obeyBounds:
225       if sup(m0) > 0:       if sup(m0) > 0:
226          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])
227        out=self.__x[0]        out=self.__x[0]
228        for i in range(1,len(self.__x)):        for i in range(1,len(X)):
229        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]
230        out = out * m0 + z * (1-m0)        out = out * m0 + z * (1-m0)
231        m0=whereNegative(x-x[i])        m0=whereNonPositive(x-X[i])
232
233          if self.__obeyBounds:
234            if inf(m0) < 1 :
235               raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
236          else:
237            out = out * m0 + X[-1] * (1-m0)
238          return out
239
240       def getValueDifferential(self, x):
241          X=self.__x
242          Y=self.__y
243
244          x0=X[0]
245          m0=whereNegative(x-x0)
246          if self.__obeyBounds:
247         if sup(m0) > 0:
248            raise ValueError,"interpolation argument out of range [%e, %e]"%(X[0],X[-1])
249          out=0.
250          for i in range(1,len(X)):
251          z=(Y[i]-Y[i-1])/(X[i]-X[i-1])
252          out = out * m0 + z * (1-m0)
253          m0=whereNonPositive(x-X[i])
254
255        if self.__obeyBounds:        if self.__obeyBounds:
256          if inf(m0) < 1:          if inf(m0) < 1:
257             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])
258        else:        else:
259          out = out * m0 + y[-1] * (1-m0)          out = out * m0
260        return out        return out
261
262
263  class Well(object):  class Well(object):
264     """     """
# Line 117  class Well(object): Line 269  class Well(object):
269     """     """
270     WATER="Water"     WATER="Water"
271     GAS="Gas"     GAS="Gas"
272     def __init__(self,  *args, **kwargs):     def __init__(self, name, domain, Q=[0.], schedule=[0.], phase="Water", BHP_limit=1.*U.atm, X0=[0.,0.,0.], *args, **kwargs):
273         """         """
274         set up well         set up well
275         """         """
276         pass         if not len(schedule) == len(Q):
277     def getGeometry(self):             raise ValueError,"length of schedule and Q must match."
278         """         self.__schedule=schedule
279         returns a function representing the geometry of the well.         self.__Q = Q
280         typically a Gaussian profile around the location of the well.         self.__phase=phase
281           self.__BHP_limit=BHP_limit
282           self.name=name
283           self.domain=domain
284           self.locator=Locator(DiracDeltaFunctions(self.domain),X0[:self.domain.getDim()])
285           self.X0=self.locator.getX()
286
287         :note: needs to be overwritten     def getLocation(self):
288         """         return self.X0
raise NotImplementedError
289
290     def getProductivityIndex(self):     def getProductivityIndex(self):
291         """         """
# Line 140  class Well(object): Line 296  class Well(object):
296         """         """
297         raise NotImplementedError         raise NotImplementedError
298
299     def getFlowRate(self):     def getFlowRate(self,t):
300        """        """
301        returns the flow rate        returns the flow rate
302        """        """
303        return 0.        out=0.
304                for i in range(len(self.__schedule)):
305     def getBHP(self):           if t <= self.__schedule[i]:
306                 out=self.__Q[i]
307                 break
308          return out
309
310       def getBHPLimit(self):
311        """        """
312        return bottom-hole pressure        return bottom-hole pressure limit
313
314        :note: needs to be overwritten        :note: needs to be overwritten
315        """        """
316        raise NotImplementedError        return self.__BHP_limit
317
318     def getPhase(self):     def getPhase(self):
319        """        """
320        returns the pahse the well works on        returns the pahse the well works on
321        """        """
322        return self.WATER        return self.__phase
323
324  class PeacemanWell(Well):  class VerticalPeacemanWell(Well):
325     """     """
326     defines a well using Peaceman's formula     defines a well using Peaceman's formula
327     """     """
328     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],
329                 perm=[1.*U.cPoise,1.*U.cPoise, 1.*U.cPoise], phase=Well.WATER, s=0, decay_factor = 5):
330                 # reset_factor=0.1):
331           """
332           set up well
333
334           :param BHP: ottom-hole pressure
335           :param Q: flow rate
336           :param r: well radius
337           :param X: location
338           :param D: dimension of well block
339           :param perm: permeabilities
340           :param s: skin factor
341           """
342           Well.__init__(self, name, domain, Q=Q, schedule=schedule, phase=phase,BHP_limit=BHP_limit, X0=X0)
343           r_el=0.28 * sqrt( sqrt(perm[1]/perm[0]) * D[0]**2 +  sqrt(perm[0]/perm[1]) * D[1]**2 )\
344                             / ( (perm[1]/perm[0])**0.25 + (perm[1]/perm[0])**0.25 )
345
346           r_el=0.11271331774384821*D[0]
347           self.__PI = 2 * math.pi * D[2] * sqrt(perm[1]*perm[0]) / (log(r_el/r) + s)
348           self.__r = r
349
350           self.__D=D
351           self.r_el=r_el
352
353           if self.domain.getDim() == 3:
354               self.geo_fac=1
355           else:
356               self.geo_fac=1./D[2]
357
358
359
360       def getProductivityIndex(self):
361           """
362           returns the productivity index of the well.
363           """
364           return self.__PI
365
366
367  class DualPorosity(object):  class DualPorosity(object):
368     """     """
# Line 176  class DualPorosity(object): Line 374  class DualPorosity(object):
374                        perm_m_0=None, perm_m_1=None, perm_m_2=None,                        perm_m_0=None, perm_m_1=None, perm_m_2=None,
375                        k_w =None, k_g=None, mu_w =None, mu_g =None,                        k_w =None, k_g=None, mu_w =None, mu_g =None,
376                        rho_w =None, rho_g=None,                        rho_w =None, rho_g=None,
377                        wells=[]):                        wells=[], g=9.81*U.m/U.sec**2):
378        """        """
379        set up        set up
380
# Line 188  class DualPorosity(object): Line 386  class DualPorosity(object):
386        :type phi_m: `MaterialPropertyWithDifferential` or None        :type phi_m: `MaterialPropertyWithDifferential` or None
387        :param phi: total porosity if None phi=1.        :param phi: total porosity if None phi=1.
388        :type phi: scalar or None        :type phi: scalar or None
389        :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
390        :type S_fg: `MaterialPropertyWithDifferential`        :type S_fg: `MaterialPropertyWithDifferential`
391        :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
392        :type S_mg: `MaterialPropertyWithDifferential`        :type S_mg: `MaterialPropertyWithDifferential`
# Line 235  class DualPorosity(object): Line 433  class DualPorosity(object):
433        self.rho_w = rho_w        self.rho_w = rho_w
434        self.rho_g = rho_g            self.rho_g = rho_g
435        self.wells=wells        self.wells=wells
436          self.t =0
437          self.g=g
438
439          self.__iter_max=1
440          self.__rtol=1.e-4
441          self.verbose=False
442          self.XI=0.5
443       def setIterationControl(self, iter_max=None, rtol=None, verbose=None, xi=None):
444         """
445         sets parameters to control iteration process
446         """
447         if iter_max !=None: self.__iter_max=iter_max
448         if rtol !=None: self.__rtol = rtol
449         if verbose !=None: self.verbose=verbose
450         if xi !=None: self.XI=xi
451
452       def update(self, dt):
453             self.u, self.u_old =tuple([ v.copy() for v in self.u ] ), self.u
454             n=0
455             converged=False
456             while n < self.__iter_max and not converged:
457                u=self.solvePDE(dt)
458                if self.verbose: print "iteration step %d:"%n
459                converged=True
460                for i in range(len(u)):
461                   if isinstance(u[i], Data):
462                 norm_u=Lsup(u[i])
463                 norm_e=Lsup(u[i]-self.u[i])
464                   else:
465                     norm_e=0.
466                     norm_u=1.
467
468               if norm_u > 0:
469                   rerr=norm_e/norm_u
470               else:
471                   rerr=norm_e
472                   if norm_e>self.__rtol * norm_u + 1.e-10: converged=False
473               if self.verbose: print "   comp %i: change = %e (value = %e)"%(i, norm_e,norm_u)
474            n+=1
475            self.u=u
476             print "iteration completed."
477             self.t+=dt
478
479
480  class PorosityOneHalfModel(DualPorosity):  class PorosityOneHalfModel(DualPorosity):
# Line 247  class PorosityOneHalfModel(DualPorosity) Line 484  class PorosityOneHalfModel(DualPorosity)
484        This corresponds to the coal bed methan model in the ECLIPSE code.        This corresponds to the coal bed methan model in the ECLIPSE code.
485        """        """
486
487        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,
488               perm_f_0=None, perm_f_1=None, perm_f_2=None,               perm_f_0=None, perm_f_1=None, perm_f_2=None,
489               k_w =None, k_g=None, mu_w =None, mu_g =None,               k_w =None, k_g=None, mu_w =None, mu_g =None,
490               rho_w =None, rho_g=None,               rho_w =None, rho_g=None, sigma=0, A_mg=0, f_rg=1.,
491                 wells=[]):                 wells=[], g=9.81*U.m/U.sec**2):
492       """       """
493       set up       set up
494
# Line 263  class PorosityOneHalfModel(DualPorosity) Line 500  class PorosityOneHalfModel(DualPorosity)
500       :type phi: scalar or None       :type phi: scalar or None
501       :param L_g: gas adsorption as function of gas pressure       :param L_g: gas adsorption as function of gas pressure
502       :type L_g: `MaterialPropertyWithDifferential`       :type L_g: `MaterialPropertyWithDifferential`
503       :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
504       :type S_fg: `MaterialPropertyWithDifferential`       :type S_fg: `MaterialPropertyWithDifferential`
505       :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
506       :type perm_f_0: scalar       :type perm_f_0: scalar
# Line 285  class PorosityOneHalfModel(DualPorosity) Line 522  class PorosityOneHalfModel(DualPorosity)
522       :type rho_g: `MaterialPropertyWithDifferential`       :type rho_g: `MaterialPropertyWithDifferential`
523       :param wells : list of wells       :param wells : list of wells
524       :type wells: list of `Well`       :type wells: list of `Well`
525         :param sigma: shape factor for gas matrix diffusion
526         :param A_mg: diffusion constant for gas adsorption
527         :param f_rg: gas re-adsorption factor
528       """       """
529
530       DualPorosity.__init__(self, domain,       DualPorosity.__init__(self, domain,
531                    phi_f=phi_f, phi_m=None, phi=phi,                    phi_f=phi_f, phi_m=None, phi=phi,
532                    S_fg=S_fg, S_mg=None,                    S_fg=None, S_mg=None,
533                    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,
534                    perm_m_0=None , perm_m_1=None, perm_m_2=None,                    perm_m_0=None , perm_m_1=None, perm_m_2=None,
535                    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,
536                    rho_w =rho_w, rho_g=rho_g,                    rho_w =rho_w, rho_g=rho_g,
537                    wells=wells)                    wells=wells, g=g)
538       self.L_g=L_g       self.L_g=L_g
539         self.sigma = sigma
540         self.A_mg = A_mg
541         self.f_rg  = f_rg
542       self.__pde=LinearPDE(self.domain, numEquations=3, numSolutions =3)       self.__pde=LinearPDE(self.domain, numEquations=3, numSolutions =3)
543
544
# Line 305  class PorosityOneHalfModel(DualPorosity) Line 548  class PorosityOneHalfModel(DualPorosity)
548       """       """
549       return self.__pde.getSolverOptions()       return self.__pde.getSolverOptions()
550
551        def solvePDE(self):        def getState(self):
552             return self.u
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]
553
554          def getOldState(self):
555         return self.u_old
556
557       phi_f   =self.phi_f.getValue(p_fg)        def setInitialState(self, p_top=1.*U.atm, p_bottom= 1.*U.atm, S_fg=0,  c_mg=None):
558       dphi_fdp=self.phi_f.getValueDifferential(p_fg)          """
559                sets the initial state
560       S_fg=  self.S_fg.getValue(p_fc)
561       dS_fgdp=   self.S_fg.getValueDifferential(p_fc)          :param p: pressure
562       S_fw=  1-S_fg          :param S_fg: gas saturation in fractured rock
563                :param c_mg: gas concentration in coal matrix at standart conditions. if not given it is calculated
564       rho_fw     = self.rho_w.getValue(p_fw)              using the  gas adsorption curve.
565       drho_fwdp  = self.rho_w.getValueDifferential(p_fw)          """
566       rho_fg = self.rho_g.getValue(p_fg)          if self.domain.getDim() == 2:
567       drho_fgdp  = self.rho_g.getValueDifferential(p_fg)             p_init=Scalar((p_top + p_bottom) /2, Solution(self.domain))
568                else:
569       L_g       = self.getValue(p_fg)             z=self.u.getDomain().getX()[0]
570       dL_gdp_fg =  self.rho_w.getValueDifferential(p_fg)             z_top=sup(z)
571                   z_bottom=inf(z)
572               p_init=(p_bottom-p_top)/(z_bottom-z_top)*(z-z_top) + p_top
573
574                S_fg_init=Scalar(S_fg, Solution(self.domain))
575
576            if c_mg == None:
577                  c_mg_init=self.L_g(p_init)
578            else:
579                  c_mg_init=Scalar(c_mg, Solution(self.domain))
580
581                q_gas=Scalar(0., DiracDeltaFunctions(self.domain))
582                q_water=Scalar(0., DiracDeltaFunctions(self.domain))
583                BHP=interpolate(p_init, DiracDeltaFunctions(self.domain))
584
585                self.u=(p_init,S_fg_init, c_mg_init, BHP, q_gas, q_water)
586
587          def solvePDE(self, dt):
588
589         p_f, S_fg, c_mg, BHP_old, q_gas_old, q_water_old =self.getState()
590         p_f_old, S_fg_old, c_mg_old, BHP, q_gas, q_water =self.getOldState()
591
592             S_fw=1-S_fg
593
594         if self.verbose:
595              print "p_f range = ",inf(p_f),sup(p_f)
596              print "S_fg range = ",inf(S_fg),sup(S_fg)
597              print "S_fw range = ",inf(S_fw),sup(S_fw)
598              print "c_mg range = ",inf(c_mg),sup(c_mg)
599                  print "BHP =",BHP
600                  print "q_gas =",q_gas
601                  print "q_water =",q_water
602
603             k_fw=self.k_w(S_fw)
604             if self.verbose: print "k_fw range = ",inf(k_fw),sup(k_fw)
605
606
607             k_fg=self.k_g(S_fg)
608             if self.verbose: print "k_fg range = ",inf(k_fg),sup(k_fg)
609
610         mu_fw=self.mu_w(p_f)
611             if self.verbose: print "mu_fw range = ",inf(mu_fw),sup(mu_fw)
612
613         mu_fg=self.mu_g(p_f)
614             if self.verbose: print "mu_fg range = ",inf(mu_fg),sup(mu_fg)
615
616
617         phi_f   =self.phi_f.getValue(p_f)
618         dphi_fdp=self.phi_f.getValueDifferential(p_f)
619             if self.verbose: print "phi_f range = ",inf(phi_f),sup(phi_f)," (slope %e,%e)"%(inf(dphi_fdp), sup(dphi_fdp))
620
621         rho_fw     = self.rho_w.getValue(p_f)
622         drho_fwdp  = self.rho_w.getValueDifferential(p_f)
623             if self.verbose: print "rho_fw range = ",inf(rho_fw),sup(rho_fw)," (slope %e,%e)"%(inf(drho_fwdp), sup(drho_fwdp))
624
625             rho_fg = self.rho_g.getValue(p_f)
626             rho_g_surf = self.rho_g.rho_surf
627         drho_fgdp  = self.rho_g.getValueDifferential(p_f)
628             if self.verbose:
629                    print "rho_fg range = ",inf(rho_fg),sup(rho_fg)," (slope %e,%e)"%(inf(drho_fgdp), sup(drho_fgdp))
630                    print "rho_fg surf = ",rho_g_surf
631
632         L_g = self.L_g(p_f)
633         dL_gdp = self.L_g.getValueDifferential(p_f)
634             if self.verbose: print "L_g range = ",inf(L_g),sup(L_g)," (slope %e,%e)"%(inf(dL_gdp), sup(dL_gdp))
635
636       A_fw = rho_fw * k_fw/mu_fw       A_fw = rho_fw * k_fw/mu_fw
637       A_fg = rho_fg * k_fg/mu_fg       A_fg = rho_fg * k_fg/mu_fg
638
639       m = whereNegative(L_g - C_mg)         m = whereNegative(L_g - c_mg)
640       B = self.sigma * self.a_mg * (m + (1-m) * self.f_rg * S_fg )       H = self.sigma * self.A_mg * (m + (1-m) * self.f_rg * S_fg )
641
642             E_fpp= S_fw * (  rho_fw * dphi_fdp + phi_f  * drho_fwdp )
643       E_fww= phi_f * ( rho_fw * dS_fgdp + S_fw * drho_fwdp )       E_fps=  -  phi_f * rho_fw
644       E_fwg= rho_fw * ( S_fw * dphi_fdp - phi_f * dS_fgdp )       E_fsp= S_fg * ( rho_fg * dphi_fdp + phi_f * drho_fgdp )
645       E_fgw= - phi_f * rho_fg * dS_fgdp       E_fss= phi_f * rho_fg
646       E_fgg= S_fg * rho_fg  * dphi_fdp + phi_f * rho_fg * dS_fgdp + phi_f  * S_fg * drho_fgdp
647
648
649         F_fw=0.
650         F_fg=0.
651         D_fw=0.
652         D_fg=0.
653
654             prod_index=Scalar(0., DiracDeltaFunctions(self.domain))
655             geo_fac=Scalar(0., DiracDeltaFunctions(self.domain))
656             for I in self.wells:
657             prod_index.setTaggedValue(I.name, I.getProductivityIndex() )
658             geo_fac.setTaggedValue(I.name, I.geo_fac)
659
660         F_fp_Y = A_fw * prod_index * BHP  * geo_fac
661         F_fs_Y = A_fg * prod_index * BHP * geo_fac
662         D_fpp =  A_fw * prod_index * geo_fac
663         D_fsp =  A_fg * prod_index * geo_fac
664
665
Q_w=0.
Q_g=0.
R_w=0.
R_g=0.

for I in self.wells:
chi_I= I.getGeometry()
Pi_I = I.getProductivityIndex()
P_I = I.getBHP()
Q_I = I.getFlowRate()

if self.getPhase == Well.WATER:
Q_w = Q_w + Q_I        * chi_I
Q_g = Q_g + Pi_I * P_I * chi_I
R_g = R_g + Pi_I       * chi_I
else:
Q_g = Q_g + Q_I        * chi_I
Q_w = Q_w + Pi_I * P_I * chi_I
R_w = R_w + Pi_I       * chi_I

F_fw_Y = A_fw * Q_w
F_fg_Y = A_fg * Q_g
D_fgg = - A_fg * R_g
D_fww = - A_fw * R_w

666       if self.domain.getDim() == 3:       if self.domain.getDim() == 3:
667          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]
668          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]
669       else:       else:
670          F_fw_X = 0 * kronecker(self.domain)[1]          F_fp_X = 0 * kronecker(self.domain)[1]
671          F_fg_X = 0 * kronecker(self.domain)[1]          F_fs_X = 0 * kronecker(self.domain)[1]
672
673       F_mg = B * ( L_g - dL_gdp_fg * p_fg )           F_mg_Y = H * L_g
674
675
676         D=self.__pde.createCoefficient("D")
677         A=self.__pde.createCoefficient("A")
678         Y=self.__pde.createCoefficient("Y")
679         X=self.__pde.createCoefficient("X")
680
681       D=self.__pde.getNewCoefficient("D")       d_dirac=self.__pde.createCoefficient("d_dirac")
682       A=self.__pde.getNewCoefficient("A")       y_dirac=self.__pde.createCoefficient("y_dirac")
683       Y=self.__pde.getNewCoefficient("Y")
684       X=self.__pde.getNewCoefficient("X")
685
686       dtXI = dt*self.XI       D[0,0]=E_fpp
687       dtcXI = dt*(1-self.XI)       D[0,1]=E_fps
688             d_dirac[0,0]=dt * D_fpp
D[0,0]=E_fww + dtXI * D_fww
D[0,1]=E_fwg
D[1,0]=E_fgw
D[1,1]=E_fgg + dtXI * D_fgg
D[1,2]=rho_fg
D[2,1]= dtXI * B * dL_gdp_fg
D[2,2]= 1 + dtXI * B
689
690       H_fw = dtcXI * A_fw * grad(p_fw_old)       D[1,0]=E_fsp
691       H_fg = dtcXI * A_fg * grad(p_fg_old)       D[1,1]=E_fss
692         D[1,2]= rho_g_surf
693         d_dirac[1,0]=dt * D_fsp
694
695       for i in range(self.domain.getDim()):       D[0,2]= -dt * H * dL_gdp * 0
696          A[0,i,0,i] = dtXI * A_fw * self.perm_f[i]       D[2,2]= 1 + dt * H
A[1,i,1,i] = dtXI * A_fg * self.perm_f[i]

X[0,i]= dt * F_fw_X - H_fw * self.perm_f[i]
X[1,i]= dt * F_fg_X - H_fg * self.perm_f[i]

Y[0] = E_fww * p_fw_old + E_fwg * p_fg_old +                     dt * F_fw_Y - dtcXI * D_fww * p_fw_old
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[2] = C_mg_old                                                + dt * F_mg_Y - dtcXI * ( dL_gdp_fg * p_fg_old + C_mg_old)
697
698       self.__pde.setValue(A=A, D=D, X=X, Y=Y)
699         for i in range(self.domain.getDim()):
700            A[0,i,0,i] = dt * A_fw * self.perm_f[i]
701            A[1,i,1,i] = dt * A_fg * self.perm_f[i]
702
703             X[0,:]=  dt * F_fp_X
704         X[1,:]=  dt * F_fs_X
705
706         Y[0] = E_fpp *  p_f_old + E_fps * S_fg_old
707         Y[1] = E_fsp *  p_f_old + E_fss * S_fg_old + rho_g_surf * c_mg_old
708         Y[2] = c_mg_old                                                    + dt * ( F_mg_Y -  H * dL_gdp * p_f *0 )
709
710         y_dirac[0] =dt * F_fp_Y
711         y_dirac[1] =dt * F_fs_Y
712
713         print "HHH D[0,0] = ",D[0,0]
714         print "HHH D[0,1] = ",D[0,1]
715         print "HHH D[0,2] = ",D[0,2]
716         print "HHH D[1,0] = ",D[1,0]
717         print "HHH D[1,1] = ",D[1,1]
718         print "HHH D[1,2] = ",D[1,2]
719         print "HHH D[2,0] = ",D[2,0]
720         print "HHH D[2,1] = ",D[2,1]
721         print "HHH D[2,2] = ",D[2,2]
722         print "HHH A_fw = ",A_fw
723         print "HHH A_fg = ",A_fg
724         print "HHH A[0,0,0,0] = ",A[0,0,0,0]
725         print "HHH A[0,1,0,1] = ",A[0,1,0,1]
726         print "HHH A[1,0,1,0] = ",A[1,0,1,0]
727         print "HHH A[1,1,1,1] = ",A[1,1,1,1]
728         print "HHH Y[0] ",Y[0]
729         print "HHH Y[1] ",Y[1]
730         print "HHH Y[2] ",Y[2]
731             print "HHH F_fp_Y ",F_fp_Y
732             print "HHH F_fs_Y ",F_fs_Y
733             print "HHH F_mg_Y ",F_mg_Y
734         print "HHH H = ",H
735
736         self.__pde.setValue(A=A, D=D, X=X, Y=Y, d_dirac=d_dirac , y_dirac=y_dirac)
737
738       u2 = self.__pde.getSolution()       u2 = self.__pde.getSolution()
739       return u2       # this is deals with round-off errors to maintain physical meaningful values
740         # we should do this in a better way to detect values that are totally wrong
741         p_f=u2[0]
742         S_fg=clip(u2[1],minval=0, maxval=1)
743         c_mg=clip(u2[2],minval=0)
744
745
746
747             q=Scalar(0., DiracDeltaFunctions(self.domain))
748             BHP_limit=Scalar(0., DiracDeltaFunctions(self.domain))
750             for I in self.wells:
751             q.setTaggedValue(I.name, I.getFlowRate(self.t+dt))
752             BHP_limit.setTaggedValue(I.name, I.getBHPLimit())
753             if I.getPhase() == Well.WATER:
755
756         p_f_wells=interpolate(p_f, DiracDeltaFunctions(self.domain))
757             S_fg_wells=interpolate(S_fg, DiracDeltaFunctions(self.domain))
758             A_fw_wells= self.k_w(1-S_fg_wells)/self.mu_w(p_f_wells)*self.rho_w(p_f_wells)
759             A_fg_wells= self.k_g(S_fg_wells)/self.mu_g(p_f_wells)*self.rho_g(p_f_wells)
760
761             BHP=clip( p_f_wells - q/prod_index * ( m * self.rho_w.rho_surf/A_fw_wells + (1-m)* self.rho_g.rho_surf/A_fg_wells ), minval = BHP_limit)
762             BHP=clip( p_f_wells - q/prod_index * self.rho_w.rho_surf/A_fw_wells, minval = BHP_limit)
763             q_gas    = prod_index * A_fg_wells * (p_f_wells - BHP)/self.rho_g.rho_surf
764             q_water = prod_index * A_fw_wells * (p_f_wells - BHP)/self.rho_w.rho_surf
765
766             return (p_f,S_fg, c_mg, BHP, q_gas, q_water )

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