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

revision 3510 by gross, Fri May 13 06:09:49 2011 UTC revision 3515 by gross, Thu May 19 08:20:57 2011 UTC
# Line 23  __url__="https://launchpad.net/escript-f Line 23  __url__="https://launchpad.net/escript-f
23  from esys.escript.linearPDEs import LinearPDE  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  from esys.escript.pdetools import Locator
26  from esys.escript import sqrt, log, whereNegative, sup, inf, whereNonPositive, integrate, Function, kronecker, grad, Lsup, clip, maximum, exp, ContinuousFunction  from esys.escript import *
27  from esys.weipa import saveVTK  from esys.weipa import saveVTK
28  import math  import math
29
# Line 304  class Well(object): Line 304  class Well(object):
304        """        """
305        returns the flow rate        returns the flow rate
306        """        """
307        if t <= self.__schedule[0]:        out=0.
308           out = self.__Q[0]        for i in range(len(self.__schedule)):
309        elif t>= self.__schedule[-1]:           if t <= self.__schedule[i]:
310           out=self.__Q[-1]               out=self.__Q[i]
311        else:               break
for i in range(1, len(self.__schedule)) :
if t < self.__schedule[i]:
out = ( self.__Q[i]-self.__Q[i-1] )/(self.__schedule[i] - self.__schedule[i-1]) * ( t - self.__schedule[i-1] ) + self.__Q[i-1]
break
312        return out        return out
313
314     def getBHPLimit(self):     def getBHPLimit(self):
# Line 343  class Well(object): Line 339  class Well(object):
339        """        """
340        return self.__phase        return self.__phase
341
def setWaterProduction(self, v):
self.water_production=v
def getWaterProduction(self):
return self.water_production
def setGasProduction(self, v):
self.gas_production=v
def getGasProduction(self):
return self.gas_production

342  class VerticalPeacemanWell(Well):  class VerticalPeacemanWell(Well):
343     """     """
344     defines a well using Peaceman's formula     defines a well using Peaceman's formula
# Line 378  class VerticalPeacemanWell(Well): Line 365  class VerticalPeacemanWell(Well):
365
366         self.__D=D         self.__D=D
367         self.r_el=r_el         self.r_el=r_el
368         if USE_NODAL_WELL:         self.locator=Locator(ContinuousFunction(self.domain),X0[:self.domain.getDim()])
369             self.X0=Locator(ContinuousFunction(self.domain),X0[:self.domain.getDim()]).getX()         self.X0=self.locator.getX()
370         else:         # this needs to be done differently
371             self.X0=X0[:self.domain.getDim()]         self.chi=whereZero(length(self.domain.getX()-self.X0))
372                 self.chi*=1./integrate(self.chi)
373         if USE_NODAL_WELL:         if self.domain.getDim() <3:
374             x=domain.getX()            self.chi*=1/self.__D[2]
self.chi = 1.
for i in range(domain.getDim()):
self.chi = self.chi * exp(-((x[i]-self.X0[i])/(0.1*D[i]))**2)
#saveVTK("t.vtu", chi=self.chi)
#1/0
else:
x=Function(domain).getX()
self.chi = 1.
if decay_factor > 0:
for i in range(domain.getDim()):
r=decay_factor * D[i]
self.chi = self.chi * clip( (D[i]/2.+r-abs(X0[i]-x[i]))/r, maxval=1., minval=0.)
else:
for i in range(domain.getDim()):
r=D[i]
self.chi = self.chi * whereNonPositive(abs(X0[i]-x[i])-D[i]/2.)

l=integrate(self.chi)
self.chi*=1./l
375

#self.chi=0.00000001
376     def getLocation(self):     def getLocation(self):
377         return self.X0         return self.X0
378
def getGeometry(self):
"""
returns a function representing the geometry of the well.
typically a Gaussian profile around the location of the well.

:note: needs to be overwritten
"""
if self.domain.getDim() <3:
return self.chi/self.__D[2]
else:
return self.chi

379     def getProductivityIndex(self):     def getProductivityIndex(self):
380         """         """
381         returns the productivity index of the well.         returns the productivity index of the well.
# Line 515  class DualPorosity(object): Line 469  class DualPorosity(object):
469       if xi !=None: self.XI=xi       if xi !=None: self.XI=xi
470
471     def update(self, dt):     def update(self, dt):
472           self.u, self.u_old = self.u.copy(), self.u           self.u, self.u_old =tuple([ v.copy() for v in self.u ] ), self.u
473           n=0           n=0
474           rerr=1.           converged=False
475           while n < self.__iter_max and rerr > self.__rtol:           while n < self.__iter_max and not converged:
476              u=self.solvePDE(dt)              u=self.solvePDE(dt)
477          norm_u=Lsup(u)              if self.verbose: print "iteration step %d:"
478          norm_e=Lsup(u-self.u)              converged=True
479                for i in range(len(u)):
480                   if isinstance(u[i], Data):
481                 norm_u=Lsup(u[i])
482                 norm_e=Lsup(u[i]-self.u[i])
483                   else:
484                     norm_e=0.
485                     norm_u=1.
486
487          if norm_u > 0:             if norm_u > 0:
488              rerr=norm_e/norm_u                 rerr=norm_e/norm_u
489          else:             else:
490              rerr=norm_e                 rerr=norm_e
491          if self.verbose: print "iteration step %d: relative change = %e"%(n, rerr)                 if rerr>self.__rtol: converged=False
492               if self.verbose: print "   comp %i: relative change = %e"%(i, rerr)
493          n+=1          n+=1
494          self.u=u          self.u=u
495           print "iteration completed."           print "iteration completed."
# Line 606  class PorosityOneHalfModel(DualPorosity) Line 568  class PorosityOneHalfModel(DualPorosity)
568       return self.__pde.getSolverOptions()       return self.__pde.getSolverOptions()
569
570        def getState(self):        def getState(self):
571       return self.u[0], self.u[1],  self.u[2]       return self.u
572
573        def getOldState(self):        def getOldState(self):
574       return self.u_old[0], self.u_old[1],  self.u_old[2]       return self.u_old
575
576        def setInitialState(self, p_top=1.*U.atm, p_bottom= 1.*U.atm, S_fg=0,  c_mg=None):        def setInitialState(self, p_top=1.*U.atm, p_bottom= 1.*U.atm, S_fg=0,  c_mg=None):
577          """          """
# Line 620  class PorosityOneHalfModel(DualPorosity) Line 582  class PorosityOneHalfModel(DualPorosity)
582          :param c_mg: gas concentration in coal matrix at standart conditions. if not given it is calculated          :param c_mg: gas concentration in coal matrix at standart conditions. if not given it is calculated
583              using the  gas adsorption curve.              using the  gas adsorption curve.
584          """              """
585          self.u=self.__pde.createSolution()          if self.domain.getDim() == 2:
586          if self.u.getDomain().getDim() == 2:             p_init=Scalar((p_top + p_bottom) /2, Solution(self.domain))
self.u[0]=(p_top + p_bottom) /2
587          else:          else:
588             z=self.u.getDomain().getX()[0]             z=self.u.getDomain().getX()[0]
589             z_top=sup(z)             z_top=sup(z)
590             z_bottom=inf(z)             z_bottom=inf(z)
591             self.u[0]=(p_bottom-p_top)/(z_bottom-z_top)*(z-z_top) + p_top             p_init=(p_bottom-p_top)/(z_bottom-z_top)*(z-z_top) + p_top
592          self.u[1]=S_fg
593                S_fg_init=Scalar(S_fg, Solution(self.domain))
594
595          if c_mg == None:          if c_mg == None:
596            self.u[2]= self.L_g(self.u[0])                c_mg_init=self.L_g(p_init)
597          else:          else:
598            self.u[2]=c_mg                c_mg_init=Scalar(c_mg, Solution(self.domain))
599
600                wells_init={}
601                for I in self.wells:
602                    BHP_limit_I=I.getBHPLimit()
603                    p_f_I=I.locator(p_init)
604                q_I = I.getFlowRate(0)
605                    if abs(q_I)>0:
606                      raise ValueError,"initial flux must be zero."
607
608                    if p_f_I < BHP_limit_I:  # fix me!
609                      raise ValueError,"initial pressure must be greater than BHP limit."
610                    wells_init[I.name]={ "bhp" : p_f_I, "q_gas": 0., "q_water": 0.}
611
612                self.u=(p_init,S_fg_init, c_mg_init, wells_init)
613
614        def solvePDE(self, dt):        def solvePDE(self, dt):
615
616       C_couple=1       p_f, S_fg, c_mg, wells_state =self.getState()
617             p_f_old, S_fg_old, c_mg_old, wells_sate_old=self.getOldState()
p_f, S_fg, c_mg=self.getState()
p_f_old, S_fg_old, c_mg_old=self.getOldState()
618
619           S_fw=1-S_fg           S_fw=1-S_fg
620
# Line 648  class PorosityOneHalfModel(DualPorosity) Line 623  class PorosityOneHalfModel(DualPorosity)
623            print "S_fg range = ",inf(S_fg),sup(S_fg)            print "S_fg range = ",inf(S_fg),sup(S_fg)
624            print "S_fw range = ",inf(S_fw),sup(S_fw)            print "S_fw range = ",inf(S_fw),sup(S_fw)
625            print "c_mg range = ",inf(c_mg),sup(c_mg)            print "c_mg range = ",inf(c_mg),sup(c_mg)
626                  print "wells state =",wells_state
627
628           k_fw=self.k_w(S_fw)           k_fw=self.k_w(S_fw)
629           if self.verbose: print "k_fw range = ",inf(k_fw),sup(k_fw)           if self.verbose: print "k_fw range = ",inf(k_fw),sup(k_fw)
630
631
632           k_fg=self.k_g(S_fg)           k_fg=self.k_g(S_fg)
633           if self.verbose: print "k_fg range = ",inf(k_fg),sup(k_fg)           if self.verbose: print "k_fg range = ",inf(k_fg),sup(k_fg)
634
# Line 696  class PorosityOneHalfModel(DualPorosity) Line 673  class PorosityOneHalfModel(DualPorosity)
673       F_fg=0.       F_fg=0.
674       D_fw=0.       D_fw=0.
675       D_fg=0.       D_fg=0.
676
677       for I in self.wells:           for I in self.wells:
678            chi_I= I.getGeometry()               Pi_I = I.getProductivityIndex()
679            Pi_I = I.getProductivityIndex()               p_f_I=I.locator(p_f)
680            BHP_limit_I=I.getBHPLimit()               S_fg_I=I.locator(S_fg)
681            p_f_I=integrate(p_f*I.chi)               BHP_I=wells_state[I.name]["bhp"]
682            S_fg_I=integrate(S_fg*I.chi)           if self.verbose: print "well %s: BHP = %e"%(I.name, BHP_I)
683            if not self.t > 0:
684                 print "AAA", self.t/U.day, p_f_I/U.psi, inf(p_f)/U.psi, S_fg_I, integrate(c_mg*I.chi)/U.Mscf*U.ft**3                 D_fw = D_fw + Pi_I *        I.chi
685            A_fw_I= self.k_w(1-S_fg_I)/self.mu_w(p_f_I)           F_fw = F_fw + Pi_I * BHP_I *I.chi
686            A_fg_I= self.k_g(S_fg_I)/self.mu_g(p_f_I)           D_fg = D_fg + Pi_I *        I.chi
687                       F_fg = F_fg + Pi_I * BHP_I *I.chi
688            if self.verbose: print "well %s is open."%I.name
689                if I.getPhase() == Well.WATER:       F_fp_Y = A_fw * F_fw
690              q_I = I.getFlowRate(self.t+dt)       F_fs_Y = A_fg * F_fg
691              p_I_ref=p_f_I-q_I/(A_fw_I * Pi_I)       D_fpp =  A_fw * D_fw
692            else:       D_fsp =  A_fg * D_fg
q_I = I.getFlowRate(self.t+dt)
p_I_ref=p_f_I-q_I/(A_fg_I * Pi_I)

if BHP_limit_I > p_I_ref:
BHP_I=BHP_limit_I
D_fw = D_fw + Pi_I * A_fw_I *        chi_I
F_fw = F_fw - Pi_I * A_fw_I * BHP_I *chi_I
D_fg = D_fg + Pi_I * A_fg_I *        chi_I
F_fg = F_fg - Pi_I * A_fg_I * BHP_I *chi_I
else:
BHP_I=p_I_ref
if I.getPhase() == Well.WATER:
F_fw = F_fw +  q_I  * chi_I
F_fg = F_fg + A_fg_I/A_fw_I * q_I * chi_I
else:
F_fg = F_fg +  q_I  * chi_I
F_fw = F_fw + A_fw_I/A_fg_I * q_I *chi_I

if self.verbose: print "well %s:BHP = %e"%(I.name, BHP_I)
I.setBHP(BHP_I)

F_fp_Y = - rho_fw * F_fw
F_fs_Y = - rho_fg * F_fg
D_fpp =  rho_fw * D_fw
D_fsp =  rho_fg * D_fg
693
694
695       if self.domain.getDim() == 3:       if self.domain.getDim() == 3:
# Line 761  class PorosityOneHalfModel(DualPorosity) Line 712  class PorosityOneHalfModel(DualPorosity)
712       D[0,1]=E_fps       D[0,1]=E_fps
713       D[1,0]=E_fsp + dt * D_fsp       D[1,0]=E_fsp + dt * D_fsp
714       D[1,1]=E_fss       D[1,1]=E_fss
#D[1,2]=rho_g_surf * C_couple
715       D[1,2]= rho_g_surf # rho_g_surf       D[1,2]= rho_g_surf # rho_g_surf
716       D[0,2]= -dt * H * dL_gdp * 0       D[0,2]= -dt * H * dL_gdp * 0
717       D[2,2]= 1 + dt * H       D[2,2]= 1 + dt * H
# Line 774  class PorosityOneHalfModel(DualPorosity) Line 724  class PorosityOneHalfModel(DualPorosity)
724           X[0,:]=  dt * F_fp_X           X[0,:]=  dt * F_fp_X
725       X[1,:]=  dt * F_fs_X       X[1,:]=  dt * F_fs_X
726
727       Y[0] = E_fpp *  p_f_old + E_fps * S_fg_old +                                  dt * F_fp_Y       Y[0] = E_fpp *  p_f_old + E_fps * S_fg_old +                         dt * F_fp_Y
728       Y[1] = E_fsp *  p_f_old + E_fss * S_fg_old + C_couple * rho_g_surf * c_mg_old + dt * F_fs_Y       Y[1] = E_fsp *  p_f_old + E_fss * S_fg_old + rho_g_surf * c_mg_old + dt * F_fs_Y
729       Y[2] = c_mg_old                                                             + dt * ( F_mg_Y -  H * dL_gdp * p_f *0 )       Y[2] = c_mg_old                                                    + dt * ( F_mg_Y -  H * dL_gdp * p_f *0 )
730
731       print "HHH D[0,0] = ",D[0,0]       print "HHH D[0,0] = ",D[0,0]
732       print "HHH D[0,1] = ",D[0,1]       print "HHH D[0,1] = ",D[0,1]
# Line 806  class PorosityOneHalfModel(DualPorosity) Line 756  class PorosityOneHalfModel(DualPorosity)
756       u2 = self.__pde.getSolution()       u2 = self.__pde.getSolution()
757       # this is deals with round-off errors to maintain physical meaningful values       # this is deals with round-off errors to maintain physical meaningful values
758       # we should do this in a better way to detect values that are totally wrong       # we should do this in a better way to detect values that are totally wrong
u2[1]=clip(u2[1],minval=0, maxval=1)
u2[2]=clip(u2[2],minval=0)

#saveVTK("t.vtu", chi = u2[0] )

print "XXX dc_mg_,t = ", (u2[2]-c_mg_old)/dt
print "XXX d_Sfg_,t = ", (u2[1]- S_fg_old)/dt

# update well production
759       p_f=u2[0]       p_f=u2[0]
760       S_fg=u2[1]       S_fg=clip(u2[1],minval=0, maxval=1)
761       for I in self.wells:       c_mg=clip(u2[2],minval=0)
762
Pi_I = I.getProductivityIndex()
q_I = I.getFlowRate(self.t+dt)
p_f_I=integrate(p_f*I.chi)
S_fg_I=integrate(S_fg*I.chi)
A_fw_I= self.k_w(1-S_fg_I)/self.mu_w(p_f_I)
A_fg_I= self.k_g(S_fg_I)/self.mu_g(p_f_I)
print "AAA", (self.t+dt)/U.day, p_f_I/U.psi, inf(p_f)/U.psi, S_fg_I, integrate(u2[2]*I.chi)/U.Mscf*U.ft**3
763
764                BHP_limit_I=I.getBHPLimit()       # update well production and BHP
765            BHP_I=I.getBHP()           wells_state_new={}
766                   for I in self.wells:
767            if self.verbose: print "reservoir pressure for well %s = %e gas)."%(I.name, p_f_I)           q_I = I.getFlowRate(self.t+dt)
768                           Pi_I = I.getProductivityIndex()
769            if BHP_limit_I < BHP_I:               p_f_I=I.locator(p_f)
770            if I.getPhase() == Well.WATER:               S_fg_I=I.locator(S_fg)
771                q_I_w =  q_I               BHP_limit_I=I.getBHPLimit()
772                q_I_g = Pi_I * A_fg_I * (p_f_I - BHP_I)
773            else:               A_fw_I= self.k_w(1-S_fg_I)/self.mu_w(p_f_I)*self.rho_w(p_f_I)
774                q_I_g =  q_I               A_fg_I= self.k_g(S_fg_I)/self.mu_g(p_f_I)*self.rho_g(p_f_I)
775                q_I_w = Pi_I * A_fw_I * (p_f_I - BHP_I)           if I.getPhase() == Well.WATER:
776            else:                   BHP= p_f_I - self.rho_w.rho_surf*q_I/(A_fw_I * Pi_I)
777              q_I_w = Pi_I * A_fw_I * (p_f_I - BHP_I)               else:
778              q_I_g = Pi_I * A_fg_I * (p_f_I - BHP_I)                   BHP= p_f_I - self.rho_g.rho_surf*q_I/(A_fg_I * Pi_I)
779            I.setWaterProduction(q_I_w)               BHP=max(BHP, BHP_limit_I)
780            I.setGasProduction(q_I_g)           if self.verbose: print "well %s: BHP = %e"%(I.name, BHP)
781       return u2
782                 wells_state_new[I.name]={ "bhp" : BHP, "q_gas": A_fg_I * Pi_I*(p_f_I-BHP)/self.rho_g.rho_surf, "q_water": A_fw_I * Pi_I*(p_f_I-BHP)/self.rho_w.rho_surf }
783
784             print wells_state_new
785             return (p_f,S_fg, c_mg, wells_state_new)

Legend:
 Removed from v.3510 changed lines Added in v.3515