/[escript]/trunk/escript/py_src/flows.py
ViewVC logotype

Diff of /trunk/escript/py_src/flows.py

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

revision 3074 by gross, Tue Jul 27 01:47:45 2010 UTC revision 3510 by gross, Fri May 13 06:09:49 2011 UTC
# Line 32  Some models for flow Line 32  Some models for flow
32    
33  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
34    
35  from escript import *  import escript
36  import util  import util
37  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
38  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
# Line 46  class DarcyFlow(object): Line 46  class DarcyFlow(object):
46        
47     where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,     where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
48        
49     :note: The problem is solved in a least squares formulation.     :cvar SIMPLE: simple solver
50       :cvar POST: solver using global postprocessing of flux
51       :cvar STAB: solver uses (non-symmetric) stabilization
52       :cvar SYMSTAB: solver uses symmetric stabilization
53     """     """
54         SIMPLE="SIMPLE"
55     def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False):     POST="POST"
56       STAB="STAB"
57       SYMSTAB="SYMSTAB"
58       def __init__(self, domain, useReduced=False, solver="SYMSTAB", verbose=False, w=1.):
59        """        """
60        initializes the Darcy flux problem        initializes the Darcy flux problem
61        :param domain: domain of the problem        :param domain: domain of the problem
62        :type domain: `Domain`        :type domain: `Domain`
63        :param useReduced: uses reduced oreder on flux and pressure        :param useReduced: uses reduced oreder on flux and pressure
64        :type useReduced: ``bool``        :type useReduced: ``bool``
65        :param adaptSubTolerance: switches on automatic subtolerance selection        :param solver: solver method
66        :type adaptSubTolerance: ``bool``        :type solver: in [`DarcyFlow.SIMPLE`, `DarcyFlow.POST', `DarcyFlow.STAB`, `DarcyFlow.SYMSTAB` ]
67        :param solveForFlux: if True the solver solves for the flux (do not use!)        :param verbose: if ``True`` some information on the iteration progress are printed.
68        :type solveForFlux: ``bool``          :type verbose: ``bool``
69          :param w: weighting factor for `DarcyFlow.POST` solver
70          :type w: ``float``
71          
72        """        """
73        self.domain=domain        self.domain=domain
74        self.solveForFlux=solveForFlux        self.solver=solver
75        self.useReduced=useReduced        self.useReduced=useReduced
76        self.__adaptSubTolerance=adaptSubTolerance        self.verbose=verbose
77        self.verbose=False        self.scale=1.
78                
79        self.__pde_k=LinearPDESystem(domain)        
80        self.__pde_k.setSymmetryOn()        self.__pde_v=LinearPDESystem(domain)
81        if self.useReduced: self.__pde_k.setReducedOrderOn()        self.__pde_v.setSymmetryOn()
82          if self.useReduced: self.__pde_v.setReducedOrderOn()
83        self.__pde_p=LinearSinglePDE(domain)        self.__pde_p=LinearSinglePDE(domain)
84        self.__pde_p.setSymmetryOn()        self.__pde_p.setSymmetryOn()
85        if self.useReduced: self.__pde_p.setReducedOrderOn()        if self.useReduced: self.__pde_p.setReducedOrderOn()
86          
87        self.__pde_l=LinearSinglePDE(domain)   # this is here for getSolverOptionsWeighting        if self.solver  == self.SIMPLE:
88        # self.__pde_l.setSymmetryOn()       if self.verbose: print "DarcyFlow: simple solver is used."
89        # if self.useReduced: self.__pde_l.setReducedOrderOn()           self.__pde_v.setValue(D=util.kronecker(self.domain.getDim()))
90          elif self.solver  == self.POST:
91         self.w=w
92         if util.inf(w)<0.:
93            raise ValueError,"Weighting factor must be non-negative."
94         if self.verbose: print "DarcyFlow: global postprocessing of flux is used."
95          elif self.solver  == self.STAB:
96          if self.verbose: print "DarcyFlow: (non-symmetric) stabilization is used."
97          elif  self.solver  == self.SYMSTAB:
98          if self.verbose: print "DarcyFlow: symmetric stabilization is used."
99          else:
100        raise ValueError,"unknown solver %s"%self.solver
101          self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
102          self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
103          self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
104          self.location_of_fixed_flux = escript.Vector(0, self.__pde_v.getFunctionSpaceForCoefficient("q"))
105        self.setTolerance()        self.setTolerance()
106        self.setAbsoluteTolerance()      
107        self.__f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))          
108        self.__g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))     def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
109          """
110          assigns values to model parameters
111    
112          :param f: volumetic sources/sinks
113          :type f: scalar value on the domain (e.g. `escript.Data`)
114          :param g: flux sources/sinks
115          :type g: vector values on the domain (e.g. `escript.Data`)
116          :param location_of_fixed_pressure: mask for locations where pressure is fixed
117          :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
118          :param location_of_fixed_flux:  mask for locations where flux is fixed.
119          :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
120          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
121          :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
122    
123          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
124          :note: at any point on the boundary of the domain the pressure
125                 (``location_of_fixed_pressure`` >0) or the normal component of the
126                 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
127                 is along the *x_i* axis.
128    
129          """
130          if location_of_fixed_pressure!=None:
131               self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)
132               self.__pde_p.setValue(q=self.location_of_fixed_pressure)
133          if location_of_fixed_flux!=None:
134              self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
135              self.__pde_v.setValue(q=self.location_of_fixed_flux)
136                
137                
138          # pressure  is rescaled by the factor 1/self.scale
139          if permeability!=None:
140        
141         perm=util.interpolate(permeability,self.__pde_v.getFunctionSpaceForCoefficient("A"))
142             V=util.vol(self.domain)
143             l=V**(1./self.domain.getDim())
144            
145         if perm.getRank()==0:
146            perm_inv=(1./perm)
147                self.scale=util.integrate(perm_inv)/V*l
148            perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
149            perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
150            
151            
152         elif perm.getRank()==2:
153            perm_inv=util.inverse(perm)
154                self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)/V)*l
155            perm_inv*=(1./self.scale)
156            perm=perm*self.scale
157         else:
158            raise ValueError,"illegal rank of permeability."
159            
160         self.__permeability=perm
161         self.__permeability_inv=perm_inv
162         if self.verbose: print "DarcyFlow: scaling factor for pressure is %e."%self.scale
163        
164         if self.solver  == self.SIMPLE:
165            self.__pde_p.setValue(A=self.__permeability)
166         elif self.solver  == self.POST:
167            self.__pde_p.setValue(A=self.__permeability)
168            k=util.kronecker(self.domain.getDim())
169            self.lamb = self.w*util.length(perm_inv)*l
170            self.__pde_v.setValue(D=self.__permeability_inv, A=self.lamb*self.domain.getSize()*util.outer(k,k))
171         elif self.solver  == self.STAB:
172            self.__pde_p.setValue(A=0.5*self.__permeability)
173            self.__pde_v.setValue(D=0.5*self.__permeability_inv)
174         elif  self.solver  == self.SYMSTAB:
175            self.__pde_p.setValue(A=0.5*self.__permeability)
176            self.__pde_v.setValue(D=0.5*self.__permeability_inv)
177    
178          if g != None:
179        g=util.interpolate(g, self.__pde_v.getFunctionSpaceForCoefficient("Y"))
180        if g.isEmpty():
181              g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
182        else:
183            if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g"
184        self.__g=g
185          if f !=None:
186         f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
187         if f.isEmpty():      
188              f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
189         else:
190             if f.getRank()>0: raise ValueError,"illegal rank of f."
191         self.__f=f
192     def getSolverOptionsFlux(self):     def getSolverOptionsFlux(self):
193        """        """
194        Returns the solver options used to solve the flux problems        Returns the solver options used to solve the flux problems
         
       *K^{-1} u=F*  
         
195        :return: `SolverOptions`        :return: `SolverOptions`
196        """        """
197        return self.__pde_k.getSolverOptions()        return self.__pde_v.getSolverOptions()
198                
199     def setSolverOptionsFlux(self, options=None):     def setSolverOptionsFlux(self, options=None):
200        """        """
201        Sets the solver options used to solve the flux problems        Sets the solver options used to solve the flux problems
         
       *K^{-1}u=F*  
         
202        If ``options`` is not present, the options are reset to default        If ``options`` is not present, the options are reset to default
         
203        :param options: `SolverOptions`        :param options: `SolverOptions`
       :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.  
204        """        """
205        return self.__pde_v.setSolverOptions(options)        return self.__pde_v.setSolverOptions(options)
206            
207     def getSolverOptionsPressure(self):     def getSolverOptionsPressure(self):
208        """        """
209        Returns the solver options used to solve the pressure problems        Returns the solver options used to solve the pressure problems
         
       *(Q^* K Q)p=-Q^*G*  
         
210        :return: `SolverOptions`        :return: `SolverOptions`
211        """        """
212        return self.__pde_p.getSolverOptions()        return self.__pde_p.getSolverOptions()
# Line 119  class DarcyFlow(object): Line 214  class DarcyFlow(object):
214     def setSolverOptionsPressure(self, options=None):     def setSolverOptionsPressure(self, options=None):
215        """        """
216        Sets the solver options used to solve the pressure problems        Sets the solver options used to solve the pressure problems
         
       *(Q^* K Q)p=-Q^*G*  
         
217        If ``options`` is not present, the options are reset to default        If ``options`` is not present, the options are reset to default
218                
219        :param options: `SolverOptions`        :param options: `SolverOptions`
220        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
221        """        """
222        return self.__pde_p.setSolverOptions(options)        return self.__pde_p.setSolverOptions(options)
   
    def getSolverOptionsWeighting(self):  
       """  
       Returns the solver options used to solve the pressure problems  
   
       *(D K D^*)p=-D F*  
   
       :return: `SolverOptions`  
       """  
       return self.__pde_l.getSolverOptions()  
   
    def setSolverOptionsWeighting(self, options=None):  
       """  
       Sets the solver options used to solve the pressure problems  
   
       *(D K D^*)p=-D F*  
   
       If ``options`` is not present, the options are reset to default  
   
       :param options: `SolverOptions`  
       :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.  
       """  
       return self.__pde_l.setSolverOptions(options)  
   
   
    def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):  
       """  
       assigns values to model parameters  
         
       :param f: volumetic sources/sinks  
       :type f: scalar value on the domain (e.g. `Data`)  
       :param g: flux sources/sinks  
       :type g: vector values on the domain (e.g. `Data`)  
       :param location_of_fixed_pressure: mask for locations where pressure is fixed  
       :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)  
       :param location_of_fixed_flux:  mask for locations where flux is fixed.  
       :type location_of_fixed_flux: vector values on the domain (e.g. `Data`)  
       :param permeability: permeability tensor. If scalar ``s`` is given the tensor with  
       ``s`` on the main diagonal is used.  
       :type permeability: scalar or tensor values on the domain (e.g. `Data`)  
       :note: the values of parameters which are not set by calling ``setValue`` are not altered.  
       :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)  
       or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal  
       is along the *x_i* axis.  
       """  
       if f !=None:  
      f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("X"))  
      if f.isEmpty():  
         f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))  
      else:  
         if f.getRank()>0: raise ValueError,"illegal rank of f."  
         self.__f=f  
       if g !=None:  
      g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))  
      if g.isEmpty():  
         g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))  
      else:  
         if not g.getShape()==(self.domain.getDim(),):  
            raise ValueError,"illegal shape of g"  
         self.__g=g  
       if location_of_fixed_pressure!=None:  
            self.__pde_p.setValue(q=location_of_fixed_pressure)  
            #self.__pde_l.setValue(q=location_of_fixed_pressure)  
       if location_of_fixed_flux!=None:  
            self.__pde_k.setValue(q=location_of_fixed_flux)  
               
       if permeability!=None:  
      perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))  
      if perm.getRank()==0:  
         perm_inv=(1./perm)*util.kronecker(self.domain.getDim())  
         perm=perm*util.kronecker(self.domain.getDim())  
      elif perm.getRank()==2:  
         perm_inv=util.inverse(perm)  
      else:  
         raise ValueError,"illegal rank of permeability."  
   
      self.__permeability=perm  
      self.__permeability_inv=perm_inv  
      self.__l =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))/10  
      #self.__l =(self.domain.getSize()**2*util.length(self.__permeability_inv))/10  
      if  self.solveForFlux:  
         self.__pde_k.setValue(D=self.__permeability_inv)  
      else:  
         self.__pde_k.setValue(D=self.__permeability_inv, A=self.__l*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))  
      self.__pde_p.setValue(A=self.__permeability)  
      #self.__pde_l.setValue(D=1/self.__l)  
          #self.__pde_l.setValue(A=self.__permeability)  
   
    def __applWeight(self, v, f=None):  
       # solves L p = f-Dv with p = 0  
       if self.getSolverOptionsWeighting().isVerbose() or self.verbose: print "DarcyFlux: Applying weighting operator"  
       if f == None:  
      return -util.div(v)*self.__l  
       else:  
      return (f-util.div(v))*self.__l  
       # if f == None:  
       #      self.__pde_l.setValue(Y=-util.div(v))    
       # else:  
       #      return (f-util.div(v))/self.__l  
       # return self.__pde_l.getSolution()  
         
    def __getPressure(self, v, p0, g=None):  
       # solves (G*KG)p = G^(g-v) with p = p0 where location_of_fixed_pressure>0  
       if self.getSolverOptionsPressure().isVerbose() or self.verbose: print "DarcyFlux: Pressure update"  
       if g == None:  
      self.__pde_p.setValue(X=-v, r=p0)  
       else:  
      self.__pde_p.setValue(X=g-v, r=p0)        
       p=self.__pde_p.getSolution()  
       return p  
   
    def __Aprod_v(self,dv):  
       # calculates: (a,b,c) = (K^{-1}(dv + KG * dp), L^{-1}Ddv, dp)  with (G*KG)dp = - G^*dv    
       dp=self.__getPressure(dv, p0=Data()) # dp = (G*KG)^{-1} (0-G^*dv)  
       a=util.tensor_mult(self.__permeability_inv,dv)+util.grad(dp) # a= K^{-1}u+G*dp  
       b= - self.__applWeight(dv) # b = - (D K D^*)^{-1} (0-Dv)  
       return ArithmeticTuple(a,b,-dp)  
   
    def __Msolve_PCG_v(self,r):  
       # K^{-1} u = r[0] + D^*r[1] = K^{-1}(dv + KG * dp) + D^*L^{-1}Ddv  
       if self.getSolverOptionsFlux().isVerbose() or self.verbose: print "DarcyFlux: Applying preconditioner"  
       self.__pde_k.setValue(X=r[1]*util.kronecker(self.domain), Y=r[0], r=Data())  
       # self.__pde_p.getOperator().saveMM("prec.mm")  
       return self.__pde_k.getSolution()  
   
    def __inner_PCG_v(self,v,r):  
       return util.integrate(util.inner(v,r[0])+util.div(v)*r[1])  
         
    def __Aprod_p(self,dp):  
       if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"  
       Gdp=util.grad(dp)  
       self.__pde_k.setValue(Y=-Gdp,X=Data(), r=Data())  
       du=self.__pde_k.getSolution()  
       # self.__pde_v.getOperator().saveMM("proj.mm")  
       return ArithmeticTuple(util.tensor_mult(self.__permeability,Gdp),-du)  
   
    def __getFlux(self,p, v0, f=None, g=None):  
       # solves (K^{-1}+D^*L^{-1} D) v = D^*L^{-1}f + K^{-1}g - Gp  
       if f!=None:  
      self.__pde_k.setValue(X=self.__applWeight(v0*0,self.__f)*util.kronecker(self.domain))  
       self.__pde_k.setValue(r=v0)  
       g2=util.tensor_mult(self.__permeability_inv,g)  
       if p == None:  
      self.__pde_k.setValue(Y=g2)  
       else:  
      self.__pde_k.setValue(Y=g2-util.grad(p))  
       return self.__pde_k.getSolution()    
223                
       #v=self.__getFlux(p, u0, f=self.__f, g=g2)        
    def __Msolve_PCG_p(self,r):  
       if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"  
       self.__pde_p.setValue(X=r[0]-r[1], Y=Data(), r=Data(), y=Data())  
       # self.__pde_p.getOperator().saveMM("prec.mm")  
       return self.__pde_p.getSolution()  
           
    def __inner_PCG_p(self,p,r):  
        return util.integrate(util.inner(util.grad(p), r[0]-r[1]))  
   
    def __L2(self,v):  
       return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))  
   
    def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):  
       """  
       solves the problem.  
         
       The iteration is terminated if the residual norm is less then self.getTolerance().  
   
       :param u0: initial guess for the flux. At locations in the domain marked by ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.  
       :type u0: vector value on the domain (e.g. `Data`).  
       :param p0: initial guess for the pressure. At locations in the domain marked by ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.  
       :type p0: scalar value on the domain (e.g. `Data`).  
       :param verbose: if set some information on iteration progress are printed  
       :type verbose: ``bool``  
       :return: flux and pressure  
       :rtype: ``tuple`` of `Data`.  
   
       :note: The problem is solved as a least squares form  
       *(K^[-1]+D^* (DKD^*)^[-1] D)u+G p=D^* (DKD^*)^[-1] f + K^[-1]g*  
       *G^*u+G^* K Gp=G^*g*  
   
       where *D* is the *div* operator and *(Gp)_i=p_{,i}* for the permeability *K=k_{ij}*.  
       """  
       self.verbose=verbose  
       rtol=self.getTolerance()  
       atol=self.getAbsoluteTolerance()  
       self.setSubProblemTolerance()  
       num_corrections=0  
       converged=False  
       norm_r=None  
         
       # Eliminate the hydrostatic pressure:  
       if self.verbose: print "DarcyFlux: calculate hydrostatic pressure component."  
       self.__pde_p.setValue(X=self.__g, r=p0, y=-util.inner(self.domain.getNormal(),u0))          
       p0=self.__pde_p.getSolution()  
       g2=self.__g - util.tensor_mult(self.__permeability, util.grad(p0))  
       norm_g2=util.integrate(util.inner(g2,util.tensor_mult(self.__permeability_inv,g2)))**0.5      
   
       p=p0*0  
       if self.solveForFlux:  
      v=u0.copy()  
       else:  
      v=self.__getFlux(p, u0, f=self.__f, g=g2)  
   
       while not converged and norm_g2 > 0:  
      Gp=util.grad(p)  
      KGp=util.tensor_mult(self.__permeability,Gp)  
      if self.verbose:  
         def_p=g2-(v+KGp)  
         def_v=self.__f-util.div(v)  
         print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(v))  
         print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(v)))  
         print "DarcyFlux: K^{-1}-norm of v = %e."%util.integrate(util.inner(v,util.tensor_mult(self.__permeability_inv,v)))**0.5  
         print "DarcyFlux: K^{-1}-norm of g2 = %e."%norm_g2  
         print "DarcyFlux: K-norm of grad(dp) = %e."%util.integrate(util.inner(Gp,KGp))**0.5  
      ATOL=atol+rtol*norm_g2  
      if self.verbose: print "DarcyFlux: absolute tolerance ATOL = %e."%(ATOL,)  
      if norm_r == None or norm_r>ATOL:  
         if num_corrections>max_num_corrections:  
            raise ValueError,"maximum number of correction steps reached."  
         
         if self.solveForFlux:  
            # initial residual is r=K^{-1}*(g-v-K*Gp)+D^*L^{-1}(f-Du)  
            v,r, norm_r=PCG(ArithmeticTuple(util.tensor_mult(self.__permeability_inv,g2-v)-Gp,self.__applWeight(v,self.__f),p),  
                self.__Aprod_v,  
                v,  
                self.__Msolve_PCG_v,  
                self.__inner_PCG_v,  
                atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)  
            p=r[2]  
         else:  
            # initial residual is r=G^*(g2-KGp - v)  
            p,r, norm_r=PCG(ArithmeticTuple(g2-KGp,v),  
                  self.__Aprod_p,  
                  p,  
                  self.__Msolve_PCG_p,  
                  self.__inner_PCG_p,  
                  atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)  
            v=r[1]  
         if self.verbose: print "DarcyFlux: residual norm = %e."%norm_r  
         num_corrections+=1  
      else:  
         if self.verbose: print "DarcyFlux: stopping criterium reached."  
         converged=True  
       return v,p+p0  
224     def setTolerance(self,rtol=1e-4):     def setTolerance(self,rtol=1e-4):
225        """        """
226        sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if        sets the relative tolerance ``rtol`` for the pressure for the stabelized solvers.
   
       *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*  
         
       where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`).  
227                
228        :param rtol: relative tolerance for the pressure        :param rtol: relative tolerance for the pressure
229        :type rtol: non-negative ``float``        :type rtol: non-negative ``float``
# Line 386  class DarcyFlow(object): Line 231  class DarcyFlow(object):
231        if rtol<0:        if rtol<0:
232       raise ValueError,"Relative tolerance needs to be non-negative."       raise ValueError,"Relative tolerance needs to be non-negative."
233        self.__rtol=rtol        self.__rtol=rtol
234          
235     def getTolerance(self):     def getTolerance(self):
236        """        """
237        returns the relative tolerance        returns the relative tolerance
# Line 393  class DarcyFlow(object): Line 239  class DarcyFlow(object):
239        :rtype: ``float``        :rtype: ``float``
240        """        """
241        return self.__rtol        return self.__rtol
242          
243     def setAbsoluteTolerance(self,atol=0.):     def solve(self,u0,p0, max_iter=100, iter_restart=20):
244        """        """
245        sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if        solves the problem.
246                
247        *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*        The iteration is terminated if the residual norm is less then self.getTolerance().
   
   
       where ``rtol`` is an absolut tolerance (see `setTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.  
248    
249        :param atol: absolute tolerance for the pressure        :param u0: initial guess for the flux. At locations in the domain marked by ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.
250        :type atol: non-negative ``float``        :type u0: vector value on the domain (e.g. `escript.Data`).
251        """        :param p0: initial guess for the pressure. At locations in the domain marked by ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.
252        if atol<0:        :type p0: scalar value on the domain (e.g. `escript.Data`).
253       raise ValueError,"Absolute tolerance needs to be non-negative."        :param max_iter: maximum number of (outer) iteration steps for the stabilization solvers,
254        self.__atol=atol        :type max_iter: ``int``
255     def getAbsoluteTolerance(self):        :param iter_restart: number of steps after which the iteration is restarted. The larger ``iter_restart`` the larger the required memory.
256        """                             A small value for ``iter_restart`` may require a large number of iteration steps or may even lead to a failure
257        returns the absolute tolerance                             of the iteration. ``iter_restart`` is relevant for the stabilization solvers only.
258        :return: current absolute tolerance        :type iter_restart: ``int``
259        :rtype: ``float``        :return: flux and pressure
260        """        :rtype: ``tuple`` of `escript.Data`.
       return self.__atol  
    def getSubProblemTolerance(self):  
       """  
       Returns a suitable subtolerance  
       :type: ``float``  
       """  
       return max(util.EPSILON**(0.5),self.getTolerance()**2)  
261    
    def setSubProblemTolerance(self):  
262        """        """
263        Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.        # rescale initial guess:
264        """        p0=p0/self.scale
265        if self.__adaptSubTolerance:        if self.solver  == self.SIMPLE or self.solver  == self.POST :
266       sub_tol=self.getSubProblemTolerance()          self.__pde_p.setValue(X=self.__g ,
267       self.getSolverOptionsFlux().setTolerance(sub_tol)                                Y=self.__f,
268       self.getSolverOptionsFlux().setAbsoluteTolerance(0.)                                y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux),
269       self.getSolverOptionsPressure().setTolerance(sub_tol)                                r=p0)
270       self.getSolverOptionsPressure().setAbsoluteTolerance(0.)          p=self.__pde_p.getSolution()
271       self.getSolverOptionsWeighting().setTolerance(sub_tol)          u = self.getFlux(p, u0)
272       self.getSolverOptionsWeighting().setAbsoluteTolerance(0.)        elif  self.solver  == self.STAB:
273       if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol      u,p = self.__solve_STAB(u0,p0, max_iter, iter_restart)
274          elif  self.solver  == self.SYMSTAB:
275        u,p = self.__solve_SYMSTAB(u0,p0, max_iter, iter_restart)
 class DarcyFlowOld(object):  
     """  
     solves the problem  
   
     *u_i+k_{ij}*p_{,j} = g_i*  
     *u_{i,i} = f*  
   
     where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,  
   
     :note: The problem is solved in a least squares formulation.  
     """  
   
     def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):  
         """  
         initializes the Darcy flux problem  
         :param domain: domain of the problem  
         :type domain: `Domain`  
     :param useReduced: uses reduced oreder on flux and pressure  
     :type useReduced: ``bool``  
     :param adaptSubTolerance: switches on automatic subtolerance selection  
     :type adaptSubTolerance: ``bool``    
         """  
         self.domain=domain  
         if weight == None:  
            s=self.domain.getSize()  
            self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2  
            # self.__l=(3.*util.longestEdge(self.domain))**2  
            #self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2  
         else:  
            self.__l=weight  
         self.__pde_v=LinearPDESystem(domain)  
         if useReduced: self.__pde_v.setReducedOrderOn()  
         self.__pde_v.setSymmetryOn()  
         self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))  
         self.__pde_p=LinearSinglePDE(domain)  
         self.__pde_p.setSymmetryOn()  
         if useReduced: self.__pde_p.setReducedOrderOn()  
         self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))  
         self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))  
         self.setTolerance()  
         self.setAbsoluteTolerance()  
     self.__adaptSubTolerance=adaptSubTolerance  
     self.verbose=False  
     def getSolverOptionsFlux(self):  
     """  
     Returns the solver options used to solve the flux problems  
       
     *(I+D^*D)u=F*  
       
     :return: `SolverOptions`  
     """  
     return self.__pde_v.getSolverOptions()  
     def setSolverOptionsFlux(self, options=None):  
     """  
     Sets the solver options used to solve the flux problems  
276            
277      *(I+D^*D)u=F*        if self.verbose:
278                KGp=util.tensor_mult(self.__permeability,util.grad(p))
279      If ``options`` is not present, the options are reset to default          def_p=self.__g-(u+KGp)
280      :param options: `SolverOptions`          def_v=self.__f-util.div(u, self.__pde_v.getFunctionSpaceForCoefficient("X"))
281      :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.          print "DarcyFlux: |g-u-K*grad(p)|_2 = %e (|u|_2 = %e)."%(self.__L2(def_p),self.__L2(u))
282      """          print "DarcyFlux: |f-div(u)|_2 = %e (|grad(u)|_2 = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))
283      return self.__pde_v.setSolverOptions(options)        #rescale result
284      def getSolverOptionsPressure(self):        p=p*self.scale
285      """        return u,p
286      Returns the solver options used to solve the pressure problems        
287           def getFlux(self,p, u0=None):
     *(Q^*Q)p=Q^*G*  
       
     :return: `SolverOptions`  
     """  
     return self.__pde_p.getSolverOptions()  
     def setSolverOptionsPressure(self, options=None):  
     """  
     Sets the solver options used to solve the pressure problems  
       
     *(Q^*Q)p=Q^*G*  
       
     If ``options`` is not present, the options are reset to default  
     :param options: `SolverOptions`  
     :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.  
     """  
     return self.__pde_p.setSolverOptions(options)  
   
     def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):  
         """  
         assigns values to model parameters  
   
         :param f: volumetic sources/sinks  
         :type f: scalar value on the domain (e.g. `Data`)  
         :param g: flux sources/sinks  
         :type g: vector values on the domain (e.g. `Data`)  
         :param location_of_fixed_pressure: mask for locations where pressure is fixed  
         :type location_of_fixed_pressure: scalar value on the domain (e.g. `Data`)  
         :param location_of_fixed_flux:  mask for locations where flux is fixed.  
         :type location_of_fixed_flux: vector values on the domain (e.g. `Data`)  
         :param permeability: permeability tensor. If scalar ``s`` is given the tensor with  
                              ``s`` on the main diagonal is used. If vector ``v`` is given the tensor with  
                              ``v`` on the main diagonal is used.  
         :type permeability: scalar, vector or tensor values on the domain (e.g. `Data`)  
   
         :note: the values of parameters which are not set by calling ``setValue`` are not altered.  
         :note: at any point on the boundary of the domain the pressure (``location_of_fixed_pressure`` >0)  
                or the normal component of the flux (``location_of_fixed_flux[i]>0`` if direction of the normal  
                is along the *x_i* axis.  
         """  
         if f !=None:  
            f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))  
            if f.isEmpty():  
                f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))  
            else:  
                if f.getRank()>0: raise ValueError,"illegal rank of f."  
            self.__f=f  
         if g !=None:  
            g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))  
            if g.isEmpty():  
              g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))  
            else:  
              if not g.getShape()==(self.domain.getDim(),):  
                raise ValueError,"illegal shape of g"  
            self.__g=g  
   
         if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)  
         if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)  
   
         if permeability!=None:  
            perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))  
            if perm.getRank()==0:  
                perm=perm*util.kronecker(self.domain.getDim())  
            elif perm.getRank()==1:  
                perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm  
                for i in range(self.domain.getDim()): perm[i,i]=perm2[i]  
            elif perm.getRank()==2:  
               pass  
            else:  
               raise ValueError,"illegal rank of permeability."  
            self.__permeability=perm  
            self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))  
   
     def setTolerance(self,rtol=1e-4):  
         """  
         sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if  
   
         *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*  
   
         where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.  
   
         :param rtol: relative tolerance for the pressure  
         :type rtol: non-negative ``float``  
         """  
         if rtol<0:  
             raise ValueError,"Relative tolerance needs to be non-negative."  
         self.__rtol=rtol  
     def getTolerance(self):  
         """  
         returns the relative tolerance  
   
         :return: current relative tolerance  
         :rtype: ``float``  
         """  
         return self.__rtol  
   
     def setAbsoluteTolerance(self,atol=0.):  
         """  
         sets the absolute tolerance ``atol`` used to terminate the solution process. The iteration is terminated if  
   
         *|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) )*  
   
         where ``rtol`` is an absolut tolerance (see `setTolerance`), *|f|^2 = integrate(length(f)^2)* and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.  
   
         :param atol: absolute tolerance for the pressure  
         :type atol: non-negative ``float``  
         """  
         if atol<0:  
             raise ValueError,"Absolute tolerance needs to be non-negative."  
         self.__atol=atol  
     def getAbsoluteTolerance(self):  
        """  
        returns the absolute tolerance  
         
        :return: current absolute tolerance  
        :rtype: ``float``  
        """  
        return self.__atol  
     def getSubProblemTolerance(self):  
     """  
     Returns a suitable subtolerance  
     @type: ``float``  
     """  
     return max(util.EPSILON**(0.75),self.getTolerance()**2)  
     def setSubProblemTolerance(self):  
          """  
          Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.  
          """  
      if self.__adaptSubTolerance:  
          sub_tol=self.getSubProblemTolerance()  
              self.getSolverOptionsFlux().setTolerance(sub_tol)  
          self.getSolverOptionsFlux().setAbsoluteTolerance(0.)  
          self.getSolverOptionsPressure().setTolerance(sub_tol)  
          self.getSolverOptionsPressure().setAbsoluteTolerance(0.)  
          if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol  
   
     def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):  
          """  
          solves the problem.  
   
          The iteration is terminated if the residual norm is less then self.getTolerance().  
   
          :param u0: initial guess for the flux. At locations in the domain marked by ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.  
          :type u0: vector value on the domain (e.g. `Data`).  
          :param p0: initial guess for the pressure. At locations in the domain marked by ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.  
          :type p0: scalar value on the domain (e.g. `Data`).  
          :param verbose: if set some information on iteration progress are printed  
          :type verbose: ``bool``  
          :return: flux and pressure  
          :rtype: ``tuple`` of `Data`.  
   
          :note: The problem is solved as a least squares form  
   
          *(I+D^*D)u+Qp=D^*f+g*  
          *Q^*u+Q^*Qp=Q^*g*  
   
          where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}* for the permeability *k_{ij}*.  
          We eliminate the flux form the problem by setting  
   
          *u=(I+D^*D)^{-1}(D^*f-g-Qp)* with u=u0 on location_of_fixed_flux  
   
          form the first equation. Inserted into the second equation we get  
   
          *Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))* with p=p0  on location_of_fixed_pressure  
   
          which is solved using the PCG method (precondition is *Q^*Q*). In each iteration step  
          PDEs with operator *I+D^*D* and with *Q^*Q* needs to be solved using a sub iteration scheme.  
          """  
          self.verbose=verbose  
          rtol=self.getTolerance()  
          atol=self.getAbsoluteTolerance()  
      self.setSubProblemTolerance()  
          num_corrections=0  
          converged=False  
          p=p0  
          norm_r=None  
          while not converged:  
                v=self.getFlux(p, fixed_flux=u0)  
                Qp=self.__Q(p)  
                norm_v=self.__L2(v)  
                norm_Qp=self.__L2(Qp)  
                if norm_v == 0.:  
                   if norm_Qp == 0.:  
                      return v,p  
                   else:  
                     fac=norm_Qp  
                else:  
                   if norm_Qp == 0.:  
                     fac=norm_v  
                   else:  
                     fac=2./(1./norm_v+1./norm_Qp)  
                ATOL=(atol+rtol*fac)  
                if self.verbose:  
                     print "DarcyFlux: L2 norm of v = %e."%norm_v  
                     print "DarcyFlux: L2 norm of k.util.grad(p) = %e."%norm_Qp  
                     print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-Qp)**2)**(0.5),)  
                     print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)  
                     print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL  
                if norm_r == None or norm_r>ATOL:  
                    if num_corrections>max_num_corrections:  
                          raise ValueError,"maximum number of correction steps reached."  
                    p,r, norm_r=PCG(self.__g-util.interpolate(v,Function(self.domain))-Qp,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.5*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)  
                    num_corrections+=1  
                else:  
                    converged=True  
          return v,p  
     def __L2(self,v):  
          return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))  
   
     def __Q(self,p):  
           return util.tensor_mult(self.__permeability,util.grad(p))  
   
     def __Aprod(self,dp):  
           if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"  
           Qdp=self.__Q(dp)  
           self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())  
           du=self.__pde_v.getSolution()  
           # self.__pde_v.getOperator().saveMM("proj.mm")  
           return Qdp+du  
     def __inner_GMRES(self,r,s):  
          return util.integrate(util.inner(r,s))  
   
     def __inner_PCG(self,p,r):  
          return util.integrate(util.inner(self.__Q(p), r))  
   
     def __Msolve_PCG(self,r):  
       if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"  
           self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())  
           # self.__pde_p.getOperator().saveMM("prec.mm")  
           return self.__pde_p.getSolution()  
   
     def getFlux(self,p=None, fixed_flux=Data()):  
288          """          """
289          returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``          returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
290          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
291          Note that ``g`` and ``f`` are used, see `setValue`.          Notice that ``g`` and ``f`` are used, see `setValue`.
292    
293          :param p: pressure.          :param p: pressure.
294          :type p: scalar value on the domain (e.g. `Data`).          :type p: scalar value on the domain (e.g. `escript.Data`).
295          :param fixed_flux: flux on the locations of the domain marked be ``location_of_fixed_flux``.          :param u0: flux on the locations of the domain marked be ``location_of_fixed_flux``.
296          :type fixed_flux: vector values on the domain (e.g. `Data`).          :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
297          :return: flux          :return: flux
298          :rtype: `Data`          :rtype: `escript.Data`
         :note: the method uses the least squares solution *u=(I+D^*D)^{-1}(D^*f-g-Qp)* where *D* is the *div* operator and *(Qp)_i=k_{ij}p_{,j}*  
                for the permeability *k_{ij}*  
299          """          """
300      self.setSubProblemTolerance()          if self.solver  == self.SIMPLE or self.solver  == self.POST  :
301          g=self.__g              KGp=util.tensor_mult(self.__permeability,util.grad(p))
302          f=self.__f              self.__pde_v.setValue(Y=self.__g-KGp, X=escript.Data())
303          self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)              if u0 == None:
304          if p == None:             self.__pde_v.setValue(r=escript.Data())
305             self.__pde_v.setValue(Y=g)          else:
306          else:             self.__pde_v.setValue(r=u0)
307             self.__pde_v.setValue(Y=g-self.__Q(p))              u= self.__pde_v.getSolution()
308          return self.__pde_v.getSolution()      elif self.solver  == self.POST:
309                self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g)-util.grad(p),
310                                      X=self.lamb * self.__f * util.kronecker(self.domain.getDim()))
311                if u0 == None:
312               self.__pde_v.setValue(r=escript.Data())
313            else:
314               self.__pde_v.setValue(r=u0)
315                u= self.__pde_v.getSolution()
316        elif self.solver  == self.STAB:
317             gp=util.grad(p)
318             self.__pde_v.setValue(Y=0.5*(util.tensor_mult(self.__permeability_inv,self.__g)+gp),
319                                   X= p * util.kronecker(self.domain.getDim()),
320                                   y= - p * self.domain.getNormal())                          
321             if u0 == None:
322               self.__pde_v.setValue(r=escript.Data())
323             else:
324               self.__pde_v.setValue(r=u0)
325             u= self.__pde_v.getSolution()
326        elif  self.solver  == self.SYMSTAB:
327             gp=util.grad(p)
328             self.__pde_v.setValue(Y=0.5*(util.tensor_mult(self.__permeability_inv,self.__g)-gp),
329                                   X= escript.Data() ,
330                                   y= escript.Data() )                          
331             if u0 == None:
332               self.__pde_v.setValue(r=escript.Data())
333             else:
334               self.__pde_v.setValue(r=u0)
335             u= self.__pde_v.getSolution()
336        return u
337          
338        
339       def __solve_STAB(self, u0, p0, max_iter, iter_restart):
340              # p0 is used as an initial guess
341          u=self.getFlux(p0, u0)  
342              self.__pde_p.setValue( Y=self.__f-util.div(u),
343                                     X=0.5*(self.__g - u - util.tensor_mult(self.__permeability,util.grad(p0)) ),
344                                     y= escript.Data(),
345                                     r=escript.Data())
346    
347          dp=self.__pde_p.getSolution()
348          p=GMRES(dp,
349                  self.__STAB_Aprod,
350              p0,
351              self.__inner,
352              atol=self.__norm(p0+dp)*self.getTolerance() ,
353              rtol=0.,
354              iter_max=max_iter,
355              iter_restart=iter_restart,
356              verbose=self.verbose,P_R=None)
357                
358              u=self.getFlux(p, u0)
359              return u,p
360    
361       def __solve_SYMSTAB(self, u0, p0, max_iter, iter_restart):
362              # p0 is used as an initial guess
363          u=self.getFlux(p0, u0)  
364              self.__pde_p.setValue( Y= self.__f,
365                                     X=  0.5*(self.__g + u - util.tensor_mult(self.__permeability,util.grad(p0)) ),
366                                     y=  -  util.inner(self.domain.getNormal(), u),
367                                     r=escript.Data())
368          dp=self.__pde_p.getSolution()
369          
370          print dp
371              print p0+dp
372              
373          p=GMRES(dp,
374                  self.__SYMSTAB_Aprod,
375              p0,
376              self.__inner,
377              atol=self.__norm(p0+dp)*self.getTolerance() ,
378              rtol=0.,
379              iter_max=max_iter,
380              iter_restart=iter_restart,
381              verbose=self.verbose,P_R=None)
382                
383              u=self.getFlux(p, u0)
384              return u,p
385    
386       def __L2(self,v):
387             return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))      
388      
389       def __norm(self,r):
390             return util.sqrt(self.__inner(r,r))
391            
392       def __inner(self,r,s):
393             return util.integrate(util.inner(r,s), escript.Function(self.domain))
394            
395       def __STAB_Aprod(self,p):
396          gp=util.grad(p)
397          self.__pde_v.setValue(Y=-0.5*gp,
398                                X=-p*util.kronecker(self.__pde_v.getDomain()),
399                                y= p * self.domain.getNormal(),  
400                                r=escript.Data())
401          u = -self.__pde_v.getSolution()
402          self.__pde_p.setValue(Y=util.div(u),
403                                X=0.5*(u+util.tensor_mult(self.__permeability,gp)),
404                                y=escript.Data(),
405                                r=escript.Data())
406        
407          return  self.__pde_p.getSolution()
408      
409       def __SYMSTAB_Aprod(self,p):
410          gp=util.grad(p)
411          self.__pde_v.setValue(Y=0.5*gp ,
412                                X=escript.Data(),
413                                y=escript.Data(),  
414                                r=escript.Data())
415          u = -self.__pde_v.getSolution()
416          self.__pde_p.setValue(Y=escript.Data(),
417                                X=0.5*(-u+util.tensor_mult(self.__permeability,gp)),
418                                y=escript.Data(),
419                                r=escript.Data())
420        
421          return  self.__pde_p.getSolution()
422          
423    
424  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
425       """       """
# Line 792  class StokesProblemCartesian(Homogeneous Line 453  class StokesProblemCartesian(Homogeneous
453           """           """
454           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
455           self.domain=domain           self.domain=domain
456           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
457           self.__pde_u.setSymmetryOn()           self.__pde_v.setSymmetryOn()
458            
459           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
460           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
# Line 810  class StokesProblemCartesian(Homogeneous Line 471  class StokesProblemCartesian(Homogeneous
471            
472       :rtype: `SolverOptions`       :rtype: `SolverOptions`
473       """       """
474       return self.__pde_u.getSolverOptions()       return self.__pde_v.getSolverOptions()
475       def setSolverOptionsVelocity(self, options=None):       def setSolverOptionsVelocity(self, options=None):
476           """           """
477       set the solver options for solving the equation for velocity.       set the solver options for solving the equation for velocity.
# Line 818  class StokesProblemCartesian(Homogeneous Line 479  class StokesProblemCartesian(Homogeneous
479       :param options: new solver  options       :param options: new solver  options
480       :type options: `SolverOptions`       :type options: `SolverOptions`
481       """       """
482           self.__pde_u.setSolverOptions(options)           self.__pde_v.setSolverOptions(options)
483       def getSolverOptionsPressure(self):       def getSolverOptionsPressure(self):
484           """           """
485       returns the solver options used  solve the equation for pressure.       returns the solver options used  solve the equation for pressure.
# Line 875  class StokesProblemCartesian(Homogeneous Line 536  class StokesProblemCartesian(Homogeneous
536          if eta !=None:          if eta !=None:
537              k=util.kronecker(self.domain.getDim())              k=util.kronecker(self.domain.getDim())
538              kk=util.outer(k,k)              kk=util.outer(k,k)
539              self.eta=util.interpolate(eta, Function(self.domain))              self.eta=util.interpolate(eta, escript.Function(self.domain))
540          self.__pde_prec.setValue(D=1/self.eta)          self.__pde_prec.setValue(D=1/self.eta)
541              self.__pde_u.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))              self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
542          if restoration_factor!=None:          if restoration_factor!=None:
543              n=self.domain.getNormal()              n=self.domain.getNormal()
544              self.__pde_u.setValue(d=restoration_factor*util.outer(n,n))              self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
545          if fixed_u_mask!=None:          if fixed_u_mask!=None:
546              self.__pde_u.setValue(q=fixed_u_mask)              self.__pde_v.setValue(q=fixed_u_mask)
547          if f!=None: self.__f=f          if f!=None: self.__f=f
548          if surface_stress!=None: self.__surface_stress=surface_stress          if surface_stress!=None: self.__surface_stress=surface_stress
549          if stress!=None: self.__stress=stress          if stress!=None: self.__stress=stress
550    
551       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data(), restoration_factor=0):       def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
552          """          """
553          assigns values to the model parameters          assigns values to the model parameters
554    
# Line 927  class StokesProblemCartesian(Homogeneous Line 588  class StokesProblemCartesian(Homogeneous
588           :return: inner product of element p and Bv=-div(v)           :return: inner product of element p and Bv=-div(v)
589           :rtype: ``float``           :rtype: ``float``
590           """           """
591           return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))           return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
592    
593       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
594           """           """
# Line 938  class StokesProblemCartesian(Homogeneous Line 599  class StokesProblemCartesian(Homogeneous
599           :return: inner product of p0 and p1           :return: inner product of p0 and p1
600           :rtype: ``float``           :rtype: ``float``
601           """           """
602           s0=util.interpolate(p0,Function(self.domain))           s0=util.interpolate(p0, escript.Function(self.domain))
603           s1=util.interpolate(p1,Function(self.domain))           s1=util.interpolate(p1, escript.Function(self.domain))
604           return util.integrate(s0*s1)           return util.integrate(s0*s1)
605    
606       def norm_v(self,v):       def norm_v(self,v):
# Line 962  class StokesProblemCartesian(Homogeneous Line 623  class StokesProblemCartesian(Homogeneous
623           :return: dv given as *Adv=(f-Av-B^*p)*           :return: dv given as *Adv=(f-Av-B^*p)*
624           """           """
625           self.updateStokesEquation(v,p)           self.updateStokesEquation(v,p)
626           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
627       self.getSolverOptionsVelocity().setTolerance(tol)       self.getSolverOptionsVelocity().setTolerance(tol)
628       self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)       self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
629           if self.__stress.isEmpty():           if self.__stress.isEmpty():
630              self.__pde_u.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))              self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
631           else:           else:
632              self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))              self.__pde_v.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
633           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
634           return  out           return  out
635    
636       def norm_Bv(self,Bv):       def norm_Bv(self,Bv):
# Line 979  class StokesProblemCartesian(Homogeneous Line 640  class StokesProblemCartesian(Homogeneous
640          :rtype: equal to the type of p          :rtype: equal to the type of p
641          :note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
642          """          """
643          return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))          return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
644    
645       def solve_AinvBt(self,p, tol):       def solve_AinvBt(self,p, tol):
646           """           """
# Line 989  class StokesProblemCartesian(Homogeneous Line 650  class StokesProblemCartesian(Homogeneous
650           :return: the solution of *Av=B^*p*           :return: the solution of *Av=B^*p*
651           :note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
652           """           """
653           self.__pde_u.setValue(Y=Data(), y=Data(), X=-p*util.kronecker(self.domain))           self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
654           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
655           return  out           return  out
656    
657       def solve_prec(self,Bv, tol):       def solve_prec(self,Bv, tol):

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

  ViewVC Help
Powered by ViewVC 1.1.26