/[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 3905 by gross, Tue Jun 5 08:33:41 2012 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 *  from . import escript
36  import util  from . 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
39    
40  class DarcyFlow(object):  class DarcyFlow(object):
41     """     """
# 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 EVAL: direct pressure gradient evaluation for flux
50       :cvar POST: global postprocessing of flux by solving the PDE *K_{ij} u_j + (w * K * l u_{k,k})_{,i}= - p_{,j} + K_{ij} g_j*
51                   where *l* is the length scale, *K* is the inverse of the permeability tensor, and *w* is a positive weighting factor.
52       :cvar SMOOTH: global smoothing by solving the PDE *K_{ij} u_j= - p_{,j} + K_{ij} g_j*
53     """     """
54         EVAL="EVAL"
55     def __init__(self, domain, useReduced=False, adaptSubTolerance=True, solveForFlux=False):     SIMPLE="EVAL"
56       POST="POST"
57       SMOOTH="SMOOTH"
58       def __init__(self, domain, useReduced=False, solver="POST", 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.EVAL`, `DarcyFlow.POST',  `DarcyFlow.SMOOTH' ]
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          if not solver in [DarcyFlow.EVAL, DarcyFlow.POST,  DarcyFlow.SMOOTH ] :
74              raise ValueError("unknown solver %d."%solver)
75    
76        self.domain=domain        self.domain=domain
77        self.solveForFlux=solveForFlux        self.solver=solver
78        self.useReduced=useReduced        self.useReduced=useReduced
79        self.__adaptSubTolerance=adaptSubTolerance        self.verbose=verbose
80        self.verbose=False        self.l=None
81                self.w=None
82        self.__pde_k=LinearPDESystem(domain)      
       self.__pde_k.setSymmetryOn()  
       if self.useReduced: self.__pde_k.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.EVAL:
88        # self.__pde_l.setSymmetryOn()           self.__pde_v=None
89        # if self.useReduced: self.__pde_l.setReducedOrderOn()           if self.verbose: print("DarcyFlow: simple solver is used.")
90        self.setTolerance()  
91        self.setAbsoluteTolerance()        elif self.solver  == self.POST:
92        self.__f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))           if util.inf(w)<0.:
93        self.__g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))              raise ValueError("Weighting factor must be non-negative.")
94                   if self.verbose: print("DarcyFlow: global postprocessing of flux is used.")
95             self.__pde_v=LinearPDESystem(domain)
96             self.__pde_v.setSymmetryOn()
97             if self.useReduced: self.__pde_v.setReducedOrderOn()
98             self.w=w
99             x=self.domain.getX()
100             self.l=min( [util.sup(x[i])-util.inf(x[i]) for i in xrange(self.domain.getDim()) ] )
101             #self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale
102    
103          elif self.solver  == self.SMOOTH:
104             self.__pde_v=LinearPDESystem(domain)
105             self.__pde_v.setSymmetryOn()
106             if self.useReduced: self.__pde_v.setReducedOrderOn()
107             if self.verbose: print("DarcyFlow: flux smoothing is used.")
108             self.w=0
109    
110          self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
111          self.__g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
112          self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
113          self.location_of_fixed_flux = escript.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
114          self.perm_scale=1.
115        
116            
117       def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
118          """
119          assigns values to model parameters
120    
121          :param f: volumetic sources/sinks
122          :type f: scalar value on the domain (e.g. `escript.Data`)
123          :param g: flux sources/sinks
124          :type g: vector values on the domain (e.g. `escript.Data`)
125          :param location_of_fixed_pressure: mask for locations where pressure is fixed
126          :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
127          :param location_of_fixed_flux:  mask for locations where flux is fixed.
128          :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
129          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
130          :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
131    
132          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
133          :note: at any point on the boundary of the domain the pressure
134                 (``location_of_fixed_pressure`` >0) or the normal component of the
135                 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
136                 is along the *x_i* axis.
137    
138          """
139          if location_of_fixed_pressure!=None:
140               self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)
141               self.__pde_p.setValue(q=self.location_of_fixed_pressure)
142          if location_of_fixed_flux!=None:
143              self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
144              if not self.__pde_v == None: self.__pde_v.setValue(q=self.location_of_fixed_flux)
145                
146          if permeability!=None:
147        
148             perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
149             self.perm_scale=util.Lsup(util.length(perm))
150             if self.verbose: print(("DarcyFlow: permeability scaling factor = %e."%self.perm_scale))
151             perm=perm*(1./self.perm_scale)
152            
153             if perm.getRank()==0:
154    
155                perm_inv=(1./perm)
156                perm_inv=perm_inv*util.kronecker(self.domain.getDim())
157                perm=perm*util.kronecker(self.domain.getDim())
158            
159            
160             elif perm.getRank()==2:
161                perm_inv=util.inverse(perm)
162             else:
163                raise ValueError("illegal rank of permeability.")
164            
165             self.__permeability=perm
166             self.__permeability_inv=perm_inv
167        
168             #====================
169             self.__pde_p.setValue(A=self.__permeability)
170             if self.solver  == self.EVAL:
171                  pass # no extra work required
172             elif self.solver  == self.POST:
173                  k=util.kronecker(self.domain.getDim())
174                  self.omega = self.w*util.length(perm_inv)*self.l*self.domain.getSize()
175                  #self.__pde_v.setValue(D=self.__permeability_inv, A=self.omega*util.outer(k,k))
176                  self.__pde_v.setValue(D=self.__permeability_inv, A_reduced=self.omega*util.outer(k,k))
177             elif self.solver  == self.SMOOTH:
178                self.__pde_v.setValue(D=self.__permeability_inv)
179    
180          if g != None:
181            g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
182            if g.isEmpty():
183                 g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
184            else:
185                 if not g.getShape()==(self.domain.getDim(),): raise ValueError("illegal shape of g")
186            self.__g=g
187          if f !=None:
188             f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
189             if f.isEmpty():      
190                 f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
191             else:
192                 if f.getRank()>0: raise ValueError("illegal rank of f.")
193             self.__f=f
194    
195     def getSolverOptionsFlux(self):     def getSolverOptionsFlux(self):
196        """        """
197        Returns the solver options used to solve the flux problems        Returns the solver options used to solve the flux problems
         
       *K^{-1} u=F*  
         
198        :return: `SolverOptions`        :return: `SolverOptions`
199        """        """
200        return self.__pde_k.getSolverOptions()        if self.__pde_v == None:
201              return None
202          else:
203              return self.__pde_v.getSolverOptions()
204                
205     def setSolverOptionsFlux(self, options=None):     def setSolverOptionsFlux(self, options=None):
206        """        """
207        Sets the solver options used to solve the flux problems        Sets the solver options used to solve the flux problems
         
       *K^{-1}u=F*  
         
208        If ``options`` is not present, the options are reset to default        If ``options`` is not present, the options are reset to default
         
209        :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.  
210        """        """
211        return self.__pde_v.setSolverOptions(options)        if not self.__pde_v == None:
212              self.__pde_v.setSolverOptions(options)
213            
214     def getSolverOptionsPressure(self):     def getSolverOptionsPressure(self):
215        """        """
216        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*  
         
217        :return: `SolverOptions`        :return: `SolverOptions`
218        """        """
219        return self.__pde_p.getSolverOptions()        return self.__pde_p.getSolverOptions()
# Line 119  class DarcyFlow(object): Line 221  class DarcyFlow(object):
221     def setSolverOptionsPressure(self, options=None):     def setSolverOptionsPressure(self, options=None):
222        """        """
223        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*  
         
224        If ``options`` is not present, the options are reset to default        If ``options`` is not present, the options are reset to default
225                
226        :param options: `SolverOptions`        :param options: `SolverOptions`
227        :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.
228        """        """
229        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  
230                
231        :param f: volumetic sources/sinks     def solve(self, u0, p0):
       :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()    
         
       #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):  
232        """        """
233        solves the problem.        solves the problem.
234                
       The iteration is terminated if the residual norm is less then self.getTolerance().  
   
235        :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.        :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.
236        :type u0: vector value on the domain (e.g. `Data`).        :type u0: vector value on the domain (e.g. `escript.Data`).
237        :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.        :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.
238        :type p0: scalar value on the domain (e.g. `Data`).        :type p0: scalar value on the domain (e.g. `escript.Data`).
       :param verbose: if set some information on iteration progress are printed  
       :type verbose: ``bool``  
239        :return: flux and pressure        :return: flux and pressure
240        :rtype: ``tuple`` of `Data`.        :rtype: ``tuple`` of `escript.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)  
241    
       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  
    def setTolerance(self,rtol=1e-4):  
242        """        """
243        sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if        self.__pde_p.setValue(X=self.__g * 1./self.perm_scale,
244                                Y=self.__f * 1./self.perm_scale,
245        *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*                              y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux * 1./self.perm_scale ),
246                                      r=p0)
247        where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`).        p=self.__pde_p.getSolution()
248                u = self.getFlux(p, u0)
249        :param rtol: relative tolerance for the pressure        return u,p
       :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  
250                
251        *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*     def getFlux(self,p, u0=None):
   
   
       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.5),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.)  
      self.getSolverOptionsWeighting().setTolerance(sub_tol)  
      self.getSolverOptionsWeighting().setAbsoluteTolerance(0.)  
      if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol  
   
   
 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  
       
     *(I+D^*D)u=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_v.setSolverOptions(options)  
     def getSolverOptionsPressure(self):  
     """  
     Returns the solver options used to solve the pressure problems  
       
     *(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):  
252          """          """
253          assigns values to model parameters          returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
   
         :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()):  
         """  
         returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``  
254          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
255          Note that ``g`` and ``f`` are used, see `setValue`.          Notice that ``g`` is used, see `setValue`.
256    
257          :param p: pressure.          :param p: pressure.
258          :type p: scalar value on the domain (e.g. `Data`).          :type p: scalar value on the domain (e.g. `escript.Data`).
259          :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``.
260          :type fixed_flux: vector values on the domain (e.g. `Data`).          :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
261          :return: flux          :return: flux
262          :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}*  
263          """          """
264      self.setSubProblemTolerance()          if self.solver  == self.EVAL:
265          g=self.__g             u = self.__g-self.perm_scale * util.tensor_mult(self.__permeability,util.grad(p))
266          f=self.__f          elif self.solver  == self.POST or self.solver  == self.SMOOTH:
267          self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)              self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g * 1./self.perm_scale)-util.grad(p))
268          if p == None:              if u0 == None:
269             self.__pde_v.setValue(Y=g)                 self.__pde_v.setValue(r=escript.Data())
270          else:              else:
271             self.__pde_v.setValue(Y=g-self.__Q(p))                 if not isinstance(u0, escript.Data) : u0 = escript.Vector(u0, escript.Solution(self.domain))
272          return self.__pde_v.getSolution()                 self.__pde_v.setValue(r=1./self.perm_scale * u0)
273                   u= self.__pde_v.getSolution() * self.perm_scale
274            return u
275          
276  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
277       """       """
278       solves       solves
# Line 778  class StokesProblemCartesian(Homogeneous Line 291  class StokesProblemCartesian(Homogeneous
291              sp.setTolerance()              sp.setTolerance()
292              sp.initialize(...)              sp.initialize(...)
293              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
294                sp.setStokesEquation(...) # new values for some parameters
295                v1,p1=sp.solve(v,p)
296       """       """
297       def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
298           """           """
# Line 792  class StokesProblemCartesian(Homogeneous Line 307  class StokesProblemCartesian(Homogeneous
307           """           """
308           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
309           self.domain=domain           self.domain=domain
310           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())
311           self.__pde_u.setSymmetryOn()           self.__pde_v.setSymmetryOn()
312            
313           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
314           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
# Line 801  class StokesProblemCartesian(Homogeneous Line 316  class StokesProblemCartesian(Homogeneous
316    
317           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
318           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
319       self.__pde_proj.setValue(D=1)           self.__pde_proj.setValue(D=1)
320           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
321    
322       def getSolverOptionsVelocity(self):       def getSolverOptionsVelocity(self):
# Line 810  class StokesProblemCartesian(Homogeneous Line 325  class StokesProblemCartesian(Homogeneous
325            
326       :rtype: `SolverOptions`       :rtype: `SolverOptions`
327       """       """
328       return self.__pde_u.getSolverOptions()           return self.__pde_v.getSolverOptions()
329       def setSolverOptionsVelocity(self, options=None):       def setSolverOptionsVelocity(self, options=None):
330           """           """
331       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 333  class StokesProblemCartesian(Homogeneous
333       :param options: new solver  options       :param options: new solver  options
334       :type options: `SolverOptions`       :type options: `SolverOptions`
335       """       """
336           self.__pde_u.setSolverOptions(options)           self.__pde_v.setSolverOptions(options)
337       def getSolverOptionsPressure(self):       def getSolverOptionsPressure(self):
338           """           """
339       returns the solver options used  solve the equation for pressure.       returns the solver options used  solve the equation for pressure.
340       :rtype: `SolverOptions`       :rtype: `SolverOptions`
341       """       """
342       return self.__pde_prec.getSolverOptions()           return self.__pde_prec.getSolverOptions()
343       def setSolverOptionsPressure(self, options=None):       def setSolverOptionsPressure(self, options=None):
344           """           """
345       set the solver options for solving the equation for pressure.       set the solver options for solving the equation for pressure.
346       :param options: new solver  options       :param options: new solver  options
347       :type options: `SolverOptions`       :type options: `SolverOptions`
348       """       """
349       self.__pde_prec.setSolverOptions(options)           self.__pde_prec.setSolverOptions(options)
350    
351       def setSolverOptionsDiv(self, options=None):       def setSolverOptionsDiv(self, options=None):
352           """           """
# Line 841  class StokesProblemCartesian(Homogeneous Line 356  class StokesProblemCartesian(Homogeneous
356       :param options: new solver options       :param options: new solver options
357       :type options: `SolverOptions`       :type options: `SolverOptions`
358       """       """
359       self.__pde_proj.setSolverOptions(options)           self.__pde_proj.setSolverOptions(options)
360       def getSolverOptionsDiv(self):       def getSolverOptionsDiv(self):
361           """           """
362       returns the solver options for solving the equation to project the divergence of       returns the solver options for solving the equation to project the divergence of
# Line 849  class StokesProblemCartesian(Homogeneous Line 364  class StokesProblemCartesian(Homogeneous
364            
365       :rtype: `SolverOptions`       :rtype: `SolverOptions`
366       """       """
367       return self.__pde_proj.getSolverOptions()           return self.__pde_proj.getSolverOptions()
368    
369       def updateStokesEquation(self, v, p):       def updateStokesEquation(self, v, p):
370           """           """
371           updates the Stokes equation to consider dependencies from ``v`` and ``p``           updates the Stokes equation to consider dependencies from ``v`` and ``p``
372           :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.           :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values to model parameters.
373           """           """
374           pass           pass
375       def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):       def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
# Line 875  class StokesProblemCartesian(Homogeneous Line 390  class StokesProblemCartesian(Homogeneous
390          if eta !=None:          if eta !=None:
391              k=util.kronecker(self.domain.getDim())              k=util.kronecker(self.domain.getDim())
392              kk=util.outer(k,k)              kk=util.outer(k,k)
393              self.eta=util.interpolate(eta, Function(self.domain))              self.eta=util.interpolate(eta, escript.Function(self.domain))
394          self.__pde_prec.setValue(D=1/self.eta)              self.__pde_prec.setValue(D=1/self.eta)
395              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)))
396          if restoration_factor!=None:          if restoration_factor!=None:
397              n=self.domain.getNormal()              n=self.domain.getNormal()
398              self.__pde_u.setValue(d=restoration_factor*util.outer(n,n))              self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
399          if fixed_u_mask!=None:          if fixed_u_mask!=None:
400              self.__pde_u.setValue(q=fixed_u_mask)              self.__pde_v.setValue(q=fixed_u_mask)
401          if f!=None: self.__f=f          if f!=None: self.__f=f
402          if surface_stress!=None: self.__surface_stress=surface_stress          if surface_stress!=None: self.__surface_stress=surface_stress
403          if stress!=None: self.__stress=stress          if stress!=None: self.__stress=stress
404    
405       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):
406          """          """
407          assigns values to the model parameters          assigns values to the model parameters
408    
# Line 900  class StokesProblemCartesian(Homogeneous Line 415  class StokesProblemCartesian(Homogeneous
415          :param surface_stress: normal surface stress          :param surface_stress: normal surface stress
416          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
417          :param stress: initial stress          :param stress: initial stress
418      :type stress: `Tensor` object on `FunctionSpace` `Function` or similar          :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
419          """          """
420          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
421    
# Line 913  class StokesProblemCartesian(Homogeneous Line 428  class StokesProblemCartesian(Homogeneous
428           :rtype: ``float``           :rtype: ``float``
429           """           """
430           self.__pde_proj.setValue(Y=-util.div(v))           self.__pde_proj.setValue(Y=-util.div(v))
431       self.getSolverOptionsDiv().setTolerance(tol)           self.getSolverOptionsDiv().setTolerance(tol)
432       self.getSolverOptionsDiv().setAbsoluteTolerance(0.)           self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
433           out=self.__pde_proj.getSolution()           out=self.__pde_proj.getSolution()
434           return out           return out
435    
# Line 927  class StokesProblemCartesian(Homogeneous Line 442  class StokesProblemCartesian(Homogeneous
442           :return: inner product of element p and Bv=-div(v)           :return: inner product of element p and Bv=-div(v)
443           :rtype: ``float``           :rtype: ``float``
444           """           """
445           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)))
446    
447       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
448           """           """
# Line 938  class StokesProblemCartesian(Homogeneous Line 453  class StokesProblemCartesian(Homogeneous
453           :return: inner product of p0 and p1           :return: inner product of p0 and p1
454           :rtype: ``float``           :rtype: ``float``
455           """           """
456           s0=util.interpolate(p0,Function(self.domain))           s0=util.interpolate(p0, escript.Function(self.domain))
457           s1=util.interpolate(p1,Function(self.domain))           s1=util.interpolate(p1, escript.Function(self.domain))
458           return util.integrate(s0*s1)           return util.integrate(s0*s1)
459    
460       def norm_v(self,v):       def norm_v(self,v):
# Line 955  class StokesProblemCartesian(Homogeneous Line 470  class StokesProblemCartesian(Homogeneous
470    
471       def getDV(self, p, v, tol):       def getDV(self, p, v, tol):
472           """           """
473           return the value for v for a given p (overwrite)           return the value for v for a given p
474    
475           :param p: a pressure           :param p: a pressure
476           :param v: a initial guess for the value v to return.           :param v: a initial guess for the value v to return.
477           :return: dv given as *Adv=(f-Av-B^*p)*           :return: dv given as *Adv=(f-Av-B^*p)*
478           """           """
479           self.updateStokesEquation(v,p)           self.updateStokesEquation(v,p)
480           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
481       self.getSolverOptionsVelocity().setTolerance(tol)           self.getSolverOptionsVelocity().setTolerance(tol)
482       self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)           self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
483           if self.__stress.isEmpty():           if self.__stress.isEmpty():
484              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)))
485           else:           else:
486              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)))
487           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
488           return  out           return  out
489    
490       def norm_Bv(self,Bv):       def norm_Bv(self,Bv):
# Line 979  class StokesProblemCartesian(Homogeneous Line 494  class StokesProblemCartesian(Homogeneous
494          :rtype: equal to the type of p          :rtype: equal to the type of p
495          :note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
496          """          """
497          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))
498    
499       def solve_AinvBt(self,p, tol):       def solve_AinvBt(self,p, tol):
500           """           """
# Line 989  class StokesProblemCartesian(Homogeneous Line 504  class StokesProblemCartesian(Homogeneous
504           :return: the solution of *Av=B^*p*           :return: the solution of *Av=B^*p*
505           :note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
506           """           """
507           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))
508           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
509           return  out           return  out
510    
511       def solve_prec(self,Bv, tol):       def solve_prec(self,Bv, tol):
# Line 1003  class StokesProblemCartesian(Homogeneous Line 518  class StokesProblemCartesian(Homogeneous
518           :note: boundary conditions on p are zero.           :note: boundary conditions on p are zero.
519           """           """
520           self.__pde_prec.setValue(Y=Bv)           self.__pde_prec.setValue(Y=Bv)
521       self.getSolverOptionsPressure().setTolerance(tol)           self.getSolverOptionsPressure().setTolerance(tol)
522       self.getSolverOptionsPressure().setAbsoluteTolerance(0.)           self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
523           out=self.__pde_prec.getSolution()           out=self.__pde_prec.getSolution()
524           return out           return out

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

  ViewVC Help
Powered by ViewVC 1.1.26