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

revision 3501 by gross, Wed Apr 13 04:07:53 2011 UTC revision 3502 by gross, Thu Apr 28 05:06:24 2011 UTC
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  from esys.escript import sqrt, log, whereNegative, sup, inf, whereNonPositive, integrate, Function, kronecker, grad, Lsup, clip, maximum, exp
27  from esys.weipa import saveVTK  from esys.weipa import saveVTK
28  import math  import math
29
# Line 233  class InterpolationTable(MaterialPropert Line 233  class InterpolationTable(MaterialPropert
233          if inf(m0) < 1 :          if inf(m0) < 1 :
234             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])
235        else:        else:
236          out = out * m0 + V[-1] * (1-m0)          out = out * m0 + X[-1] * (1-m0)
237        return out        return out
238
239     def getValueDifferential(self, x):     def getValueDifferential(self, x):
# Line 364  class VerticalPeacemanWell(Well): Line 364  class VerticalPeacemanWell(Well):
364     defines a well using Peaceman's formula     defines a well using Peaceman's formula
365     """     """
366     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],     def __init__(self,name, domain, schedule = [0.], BHP_limit=1.*U.atm, Q=0, r=10*U.cm, X0=[0.,0.,0.], D=[1.*U.km,1.*U.km, 1.*U.m],
367               perm=[1.*U.cPoise,1.*U.cPoise, 1.*U.cPoise], phase=Well.WATER, s=0):               perm=[1.*U.cPoise,1.*U.cPoise, 1.*U.cPoise], phase=Well.WATER, s=0, decay_factor = 5.):
368                 # reset_factor=0.1):
369         """         """
370         set up well         set up well
371
# Line 388  class VerticalPeacemanWell(Well): Line 389  class VerticalPeacemanWell(Well):
389
390         x=Function(domain).getX()         x=Function(domain).getX()
391         self.chi = 1.         self.chi = 1.
392           decay_factor
393         for i in range(domain.getDim()):         for i in range(domain.getDim()):
394          self.chi = self.chi * whereNegative(abs(x[i]-X0[i])-D[i]/2)          r=decay_factor * D[i]
395            self.chi = self.chi * clip( (D[i]/2.+r-abs(X0[i]-x[i]))/r, maxval=1., minval=0.)
396            #if reset_factor > 0:
397            #     f=maximum(0., abs(x[i]-X0[i])-D[i]/2*(1.-reset_factor))
398            #     self.chi = self.chi * exp(-(f/(reset_factor*D[i]/2))**2/2)
399            #else:
400                #     self.chi = self.chi * whereNegative(abs(x[i]-X0[i])-D[i]/2)
401
402         self.chi*=1./(D[0]*D[1]*D[2])         l=integrate(self.chi)
403                 self.chi*=1./l
404
405         #self.chi=0.00000001         #self.chi=0.00000001
406     def getLocation(self):     def getLocation(self):
# Line 405  class VerticalPeacemanWell(Well): Line 413  class VerticalPeacemanWell(Well):
413
414         :note: needs to be overwritten         :note: needs to be overwritten
415         """         """
416         return self.chi         if self.domain.getDim() <3:
417              return self.chi/self.__D[2]
418           else:
419          return self.chi
420
421     def getProductivityIndex(self):     def getProductivityIndex(self):
422         """         """
# Line 651  class PorosityOneHalfModel(DualPorosity) Line 662  class PorosityOneHalfModel(DualPorosity)
662           rho_fg = self.rho_g.getValue(p_f)           rho_fg = self.rho_g.getValue(p_f)
663       drho_fgdp  = self.rho_g.getValueDifferential(p_f)       drho_fgdp  = self.rho_g.getValueDifferential(p_f)
664           if self.verbose: print "rho_fg range = ",inf(rho_fg),sup(rho_fg)," (slope %e,%e)"%(inf(drho_fgdp), sup(drho_fgdp))           if self.verbose: print "rho_fg range = ",inf(rho_fg),sup(rho_fg)," (slope %e,%e)"%(inf(drho_fgdp), sup(drho_fgdp))
665
666
667
668       L_g_0       = self.L_g(p_f)       L_g_0       = self.L_g(p_f)
669       FF_g=self.rho_g.getFormationFactor(p_f)       FF_g=self.rho_g.getFormationFactor(p_f)
670       L_g = L_g_0 * FF_g       L_g = L_g_0 * FF_g
# Line 683  class PorosityOneHalfModel(DualPorosity) Line 696  class PorosityOneHalfModel(DualPorosity)
696
697       for I in self.wells:       for I in self.wells:
698            chi_I= I.getGeometry()            chi_I= I.getGeometry()
loc=Locator(Function(self.domain),I.getLocation())
699            Pi_I = I.getProductivityIndex()            Pi_I = I.getProductivityIndex()
A_fw_I= loc(A_fw)
A_fg_I= loc(A_fg)
700            BHP_limit_I=I.getBHPLimit()            BHP_limit_I=I.getBHPLimit()
701              p_f_I=integrate(p_f*I.chi)
702              S_fg_I=integrate(S_fg*I.chi)
703
704              A_fw_I= self.rho_w.getValue(p_f_I) * self.k_w(1-S_fg_I)/self.mu_w(p_f_I)
705              A_fg_I= self.rho_g.getValue(p_f_I) * self.k_g(S_fg_I)/self.mu_g(p_f_I)
706
707            if I.isOpen(self.t+dt):            if I.isOpen(self.t+dt):
708          if self.verbose: print "well %s is open."%I.name          if self.verbose: print "well %s is open."%I.name
709          if I.getPhase() == Well.WATER:          if I.getPhase() == Well.WATER:
710              q_I = self.rho_w.rho_0 * I.getFlowRate()              q_I = self.rho_w.rho_0 * I.getFlowRate()
711              p_I_ref=loc(p_f)-q_I/(A_fw_I * Pi_I)              p_I_ref=p_f_I-q_I/(A_fw_I * Pi_I)
712          else:          else:
713              q_I = self.rho_g.rho_0 * I.getFlowRate()              q_I = self.rho_g.rho_0 * I.getFlowRate()
714              p_I_ref=loc(p_f)-q_I/(A_fg_I * Pi_I)              p_I_ref=p_f_I-q_I/(A_fg_I * Pi_I)
715
716          print "ZZZ =",loc(p_f), q_I, self.rho_w.rho_0, I.getFlowRate(), A_fw_I, Pi_I, q_I/(A_fw_I * Pi_I)          print "ZZZ =",p_f_I, q_I, self.rho_w.rho_0, I.getFlowRate(), A_fw_I, Pi_I, q_I/(A_fw_I * Pi_I)
717
718          if BHP_limit_I > p_I_ref:          if BHP_limit_I > p_I_ref:
719            BHP_I=BHP_limit_I            BHP_I=BHP_limit_I
720            D_fw = D_fw + Pi_I * A_fw_I *              chi_I            D_fw = D_fw + Pi_I * A_fw_I *        chi_I
721            F_fw = F_fw - Pi_I * A_fw_I * BHP_limit_I *chi_I            F_fw = F_fw - Pi_I * A_fw_I * BHP_I *chi_I
722            D_fg = D_fg + Pi_I * A_fg_I *              chi_I            D_fg = D_fg + Pi_I * A_fg_I *        chi_I
723            F_fg = F_fg - Pi_I * A_fg_I * BHP_limit_I *chi_I            F_fg = F_fg - Pi_I * A_fg_I * BHP_I *chi_I
724          else:          else:
725              BHP_I=p_I_ref
726            if I.getPhase() == Well.WATER:            if I.getPhase() == Well.WATER:
727                F_fw = F_fw +  q_I  * chi_I                F_fw = F_fw +  q_I  * chi_I
728                F_fg = F_fg +  A_fg_I/A_fw_I *  q_I *chi_I                D_fg = D_fg + Pi_I * A_fg_I *         chi_I
729                  F_fg = F_fg - Pi_I * A_fg_I * BHP_I * chi_I
730            else:            else:
731                F_fg = F_fg +  q_I  * chi_I                F_fg = F_fg +  q_I  * chi_I
732                F_fw = F_fw +  A_fw_I/A_fg_I *  q_I *chi_I                D_fw = D_fw + Pi_I * A_fw_I *        chi_I
733            BHP_I=p_I_ref                F_fw = F_fw - Pi_I * A_fw_I * BHP_I *chi_I
734
735            else:            else:
736            if self.verbose: print "well %s is closed."%I.name            if self.verbose: print "well %s is closed."%I.name
737            BHP_I=loc(p_f)            BHP_I=p_f_I
738
739            if self.verbose: print "well %s:BHP = %e"%(I.name, BHP_I)            if self.verbose: print "well %s:BHP = %e"%(I.name, BHP_I)
740            I.setBHP(BHP_I)            I.setBHP(BHP_I)
# Line 789  class PorosityOneHalfModel(DualPorosity) Line 807  class PorosityOneHalfModel(DualPorosity)
807
808       # update well production       # update well production
809       p_f=u2[0]       p_f=u2[0]
810         S_fg=u2[1]
811       for I in self.wells:       for I in self.wells:
loc=Locator(Function(self.domain),I.getLocation())
812            Pi_I = I.getProductivityIndex()            Pi_I = I.getProductivityIndex()
813            q_I = I.getFlowRate()            q_I = I.getFlowRate()
814            A_fw_I= loc(A_fw)            p_f_I=integrate(p_f*I.chi)
815            A_fg_I= loc(A_fg)            S_fg_I=integrate(S_fg*I.chi)
816            p_f_I=loc(p_f)            rho_fw_I = self.rho_w.getValue(p_f_I)
817            BHP_limit_I=I.getBHPLimit()            rho_fg_I = self.rho_g.getValue(p_f_I)
818            BHP=I.getBHP()            A_fw_I= rho_fw_I * self.k_w(1-S_fg_I)/self.mu_w(p_f_I)
819                        A_fg_I= rho_fg_I * self.k_g(S_fg_I)/self.mu_g(p_f_I)
820
821                  BHP_limit_I=I.getBHPLimit()
822              BHP_I=I.getBHP()
823
824
825            rho_fw_I = self.rho_w.getValue(p_f_I)
rho_fg_I = self.rho_g.getValue(p_f_I)
826            if self.verbose: print "reservoir pressure for well %s = %e gas rho= %e)."%(I.name, p_f_I, rho_fg_I)            if self.verbose: print "reservoir pressure for well %s = %e gas rho= %e)."%(I.name, p_f_I, rho_fg_I)
827
828            if I.isOpen(self.t+dt):            if I.isOpen(self.t+dt):
829          if BHP_limit_I < BHP:          if BHP_limit_I < BHP_I:
830            if I.getPhase() == Well.WATER:            if I.getPhase() == Well.WATER:
831                q_I_w =  q_I                q_I_w =  q_I
832                q_I_g =  A_fg_I/A_fw_I *  q_I                q_I_g = Pi_I * A_fg_I * (p_f_I - BHP_I)  /self.rho_g.rho_0
print "ASDFASDFDF ", q_I_g/rho_fg_I, A_fg_I/A_fw_I
833            else:            else:
834                q_I_g =  q_I                q_I_g =  q_I
835                q_I_w =  A_fw_I/A_fg_I *  q_I                q_I_w = Pi_I * A_fw_I * (p_f_I - BHP_I)/ self.rho_w.rho_0
836          else:          else:
837              q_I_w = Pi_I * A_fw_I * (p_f_I - BHP_I)/ self.rho_w.rho_0              q_I_w = Pi_I * A_fw_I * (p_f_I - BHP_I)/ self.rho_w.rho_0
838              q_I_g = Pi_I * A_fg_I * (p_f_I - BHP_I)  /self.rho_g.rho_0              q_I_g = Pi_I * A_fg_I * (p_f_I - BHP_I)  /self.rho_g.rho_0

Legend:
 Removed from v.3501 changed lines Added in v.3502