/[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 3523 by gross, Thu May 19 08:20:57 2011 UTC revision 3524 by gross, Tue May 24 02:23:00 2011 UTC
# Line 269  class Well(object): Line 269  class Well(object):
269     """     """
270     WATER="Water"     WATER="Water"
271     GAS="Gas"     GAS="Gas"
272     def __init__(self, name, domain, Q=[0.], schedule=[0.], phase="Water", BHP_limit=1.*U.atm, *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         """         """
# Line 281  class Well(object): Line 281  class Well(object):
281         self.__BHP_limit=BHP_limit         self.__BHP_limit=BHP_limit
282         self.name=name         self.name=name
283         self.domain=domain         self.domain=domain
284           self.locator=Locator(ContinuousFunction(self.domain),X0[:self.domain.getDim()])
285           self.X0=self.locator.getX()
286                
287     def getGeometry(self):     def getLocation(self):
288         """         return self.X0      
        returns a function representing the geometry of the well.  
        typically a Gaussian profile around the location of the well.  
         
        :note: needs to be overwritten  
        """  
        raise NotImplementedError  
289    
290     def getProductivityIndex(self):     def getProductivityIndex(self):
291         """         """
# Line 319  class Well(object): Line 315  class Well(object):
315        """        """
316        return self.__BHP_limit        return self.__BHP_limit
317                
    def getBHP(self):  
       """  
       return bottom-hole pressure  
         
       :note: needs to be overwritten  
       """  
       return self.__BHP  
     
    def setBHP(self, BHP):  
       """  
       sets bottom-hole pressure  
       """  
       self.__BHP= BHP  
         
318     def getPhase(self):     def getPhase(self):
319        """        """
320        returns the pahse the well works on        returns the pahse the well works on
# Line 357  class VerticalPeacemanWell(Well): Line 339  class VerticalPeacemanWell(Well):
339         :param perm: permeabilities         :param perm: permeabilities
340         :param s: skin factor         :param s: skin factor
341         """         """
342         Well.__init__(self, name, domain, Q=Q, schedule=schedule, phase=phase,BHP_limit=BHP_limit)         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 )\         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 )                           / ( (perm[1]/perm[0])**0.25 + (perm[1]/perm[0])**0.25 )
345         self.__PI = 2 * math.pi * D[2] * sqrt(perm[1]*perm[0]) / (log(r_el/r) + s)         self.__PI = 2 * math.pi * D[2] * sqrt(perm[1]*perm[0]) / (log(r_el/r) + s)
# Line 365  class VerticalPeacemanWell(Well): Line 347  class VerticalPeacemanWell(Well):
347    
348         self.__D=D         self.__D=D
349         self.r_el=r_el         self.r_el=r_el
        self.locator=Locator(ContinuousFunction(self.domain),X0[:self.domain.getDim()])  
        self.X0=self.locator.getX()  
        # this needs to be done differently  
        self.chi=whereZero(length(self.domain.getX()-self.X0))  
        self.chi*=1./integrate(self.chi)  
        if self.domain.getDim() <3:  
           self.chi*=1/self.__D[2]  
350    
    def getLocation(self):  
        return self.X0  
351                
352     def getProductivityIndex(self):     def getProductivityIndex(self):
353         """         """
# Line 668  class PorosityOneHalfModel(DualPorosity) Line 641  class PorosityOneHalfModel(DualPorosity)
641       E_fps=  -  phi_f * rho_fw       E_fps=  -  phi_f * rho_fw
642       E_fsp= S_fg * ( rho_fg * dphi_fdp + phi_f * drho_fgdp )       E_fsp= S_fg * ( rho_fg * dphi_fdp + phi_f * drho_fgdp )
643       E_fss= phi_f * rho_fg       E_fss= phi_f * rho_fg
644        
645        
646            
647       F_fw=0.       F_fw=0.
648       F_fg=0.       F_fg=0.
649       D_fw=0.       D_fw=0.
650       D_fg=0.       D_fg=0.
651    
652             prod_index=Scalar(0., DiracDeltaFunctions(self.domain))
653           for I in self.wells:           for I in self.wells:
654               Pi_I = I.getProductivityIndex()           prod_index.setTaggedValue(I.name, I.getProductivityIndex() )
655               p_f_I=I.locator(p_f)  
656               S_fg_I=I.locator(S_fg)       F_fp_Y = A_fw * prod_index * wells_state["bhp"]
657               BHP_I=wells_state[I.name]["bhp"]       F_fs_Y = A_fg * prod_index * wells_state["bhp"]
658           if self.verbose: print "well %s: BHP = %e"%(I.name, BHP_I)       D_fpp =  A_fw * prod_index
659         D_fsp =  A_fg * prod_index
          D_fw = D_fw + Pi_I *        I.chi  
          F_fw = F_fw + Pi_I * BHP_I *I.chi  
          D_fg = D_fg + Pi_I *        I.chi  
          F_fg = F_fg + Pi_I * BHP_I *I.chi  
   
      F_fp_Y = A_fw * F_fw  
      F_fs_Y = A_fg * F_fg  
      D_fpp =  A_fw * D_fw  
      D_fsp =  A_fg * D_fg  
660    
661            
662       if self.domain.getDim() == 3:       if self.domain.getDim() == 3:
# Line 706  class PorosityOneHalfModel(DualPorosity) Line 673  class PorosityOneHalfModel(DualPorosity)
673       A=self.__pde.createCoefficient("A")       A=self.__pde.createCoefficient("A")
674       Y=self.__pde.createCoefficient("Y")       Y=self.__pde.createCoefficient("Y")
675       X=self.__pde.createCoefficient("X")       X=self.__pde.createCoefficient("X")
676            
677         d_dirac=self.__pde.createCoefficient("d_dirac")
678         y_dirac=self.__pde.createCoefficient("y_dirac")
679    
680    
681       D[0,0]=E_fpp + dt * D_fpp          
682         D[0,0]=E_fpp
683       D[0,1]=E_fps       D[0,1]=E_fps
684       D[1,0]=E_fsp + dt * D_fsp       d_dirac[0,0]=dt * D_fpp
685        
686         D[1,0]=E_fsp
687       D[1,1]=E_fss       D[1,1]=E_fss
688       D[1,2]= rho_g_surf # rho_g_surf       D[1,2]= rho_g_surf
689         d_dirac[1,0]=dt * D_fsp
690        
691       D[0,2]= -dt * H * dL_gdp * 0       D[0,2]= -dt * H * dL_gdp * 0
692       D[2,2]= 1 + dt * H       D[2,2]= 1 + dt * H
693            
# Line 724  class PorosityOneHalfModel(DualPorosity) Line 699  class PorosityOneHalfModel(DualPorosity)
699           X[0,:]=  dt * F_fp_X           X[0,:]=  dt * F_fp_X
700       X[1,:]=  dt * F_fs_X       X[1,:]=  dt * F_fs_X
701    
702       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
703       Y[1] = E_fsp *  p_f_old + E_fss * S_fg_old + 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
704       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 )
705            
706         y_dirac[0] =dt * F_fp_Y
707         y_dirac[1] =dt * F_fs_Y
708        
709       print "HHH D[0,0] = ",D[0,0]       print "HHH D[0,0] = ",D[0,0]
710       print "HHH D[0,1] = ",D[0,1]       print "HHH D[0,1] = ",D[0,1]
711       print "HHH D[0,2] = ",D[0,2]       print "HHH D[0,2] = ",D[0,2]
# Line 751  class PorosityOneHalfModel(DualPorosity) Line 729  class PorosityOneHalfModel(DualPorosity)
729           print "HHH F_mg_Y ",F_mg_Y           print "HHH F_mg_Y ",F_mg_Y
730       print "HHH H = ",H       print "HHH H = ",H
731    
732       self.__pde.setValue(A=A, D=D, X=X, Y=Y)       self.__pde.setValue(A=A, D=D, X=X, Y=Y, d_dirac=d_dirac , y_dirac=y_dirac)
733            
734       u2 = self.__pde.getSolution()       u2 = self.__pde.getSolution()
735       # this is deals with round-off errors to maintain physical meaningful values       # this is deals with round-off errors to maintain physical meaningful values
# Line 761  class PorosityOneHalfModel(DualPorosity) Line 739  class PorosityOneHalfModel(DualPorosity)
739       c_mg=clip(u2[2],minval=0)       c_mg=clip(u2[2],minval=0)
740                    
741    
      # update well production and BHP  
          wells_state_new={}  
      for I in self.wells:  
          q_I = I.getFlowRate(self.t+dt)  
              Pi_I = I.getProductivityIndex()  
              p_f_I=I.locator(p_f)  
              S_fg_I=I.locator(S_fg)  
              BHP_limit_I=I.getBHPLimit()  
742    
743               A_fw_I= self.k_w(1-S_fg_I)/self.mu_w(p_f_I)*self.rho_w(p_f_I)           q=Scalar(0., DiracDeltaFunctions(self.domain))
744               A_fg_I= self.k_g(S_fg_I)/self.mu_g(p_f_I)*self.rho_g(p_f_I)           BHP_limit=Scalar(0., DiracDeltaFunctions(self.domain))
745             water_well_mask=Scalar(0., DiracDeltaFunctions(self.domain))
746             for I in self.wells:
747             q.setTaggedValue(I.name, I.getFlowRate(self.t+dt))
748             BHP_limit.setTaggedValue(I.name, I.getBHPLimit())
749           if I.getPhase() == Well.WATER:           if I.getPhase() == Well.WATER:
750                   BHP= p_f_I - self.rho_w.rho_surf*q_I/(A_fw_I * Pi_I)                water_well_mask.setTaggedValue(I.name, 1)
751               else:      
752                   BHP= p_f_I - self.rho_g.rho_surf*q_I/(A_fg_I * Pi_I)       p_f_wells=interpolate(p_f, DiracDeltaFunctions(self.domain))
753               BHP=max(BHP, BHP_limit_I)           S_fg_wells=interpolate(S_fg, DiracDeltaFunctions(self.domain))
754           if self.verbose: print "well %s: BHP = %e"%(I.name, BHP)           A_fw_wells= self.k_w(1-S_fg_wells)/self.mu_w(p_f_wells)*self.rho_w(p_f_wells)
755                         A_fg_wells= self.k_g(S_fg_wells)/self.mu_g(p_f_wells)*self.rho_g(p_f_wells)
756               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 }          
757             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)
758                    
759             wells_state_new={ "bhp" : BHP, "q_gas": A_fg_wells * prod_index*(p_f_wells-BHP)/self.rho_g.rho_surf, "q_water": A_fw_wells * prod_index*(p_f_wells-BHP)/self.rho_w.rho_surf }
760                
761           print wells_state_new           print wells_state_new
762           return (p_f,S_fg, c_mg, wells_state_new)           return (p_f,S_fg, c_mg, wells_state_new)

Legend:
Removed from v.3523  
changed lines
  Added in v.3524

  ViewVC Help
Powered by ViewVC 1.1.26