/[escript]/trunk/finley/test/python/coalgas.py
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

  ViewVC Help
Powered by ViewVC 1.1.26