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

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

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

revision 3500 by gross, Mon Apr 11 03:09:06 2011 UTC revision 3990 by caltinay, Tue Sep 25 05:03:20 2012 UTC
# Line 1  Line 1 
1  # -*- coding: utf-8 -*-  # -*- coding: utf-8 -*-
2  ########################################################  ##############################################################################
3  #  #
4  # Copyright (c) 2003-2010 by University of Queensland  # Copyright (c) 2003-2012 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # http://www.uq.edu.au
 # http://www.uq.edu.au/esscc  
6  #  #
7  # Primary Business: Queensland, Australia  # Primary Business: Queensland, Australia
8  # Licensed under the Open Software License version 3.0  # Licensed under the Open Software License version 3.0
9  # http://www.opensource.org/licenses/osl-3.0.php  # http://www.opensource.org/licenses/osl-3.0.php
10  #  #
11  ########################################################  # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    # Development since 2012 by School of Earth Sciences
13    #
14    ##############################################################################
15    
16  __copyright__="""Copyright (c) 2003-2010 by University of Queensland  __copyright__="""Copyright (c) 2003-2012 by University of Queensland
17  Earth Systems Science Computational Center (ESSCC)  http://www.uq.edu.au
 http://www.uq.edu.au/esscc  
18  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
19  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
20  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
# Line 32  Some models for flow Line 33  Some models for flow
33    
34  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
35    
36  import escript  from . import escript
37  import util  from . import util
38  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions  from .linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
39  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES  from .pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
   
40    
41  class DarcyFlowVeryNew(object):  class DarcyFlow(object):
42     """     """
43     solves the problem     solves the problem
44        
# Line 47  class DarcyFlowVeryNew(object): Line 47  class DarcyFlowVeryNew(object):
47        
48     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,
49        
50     :note: The problem is solved in a stabelized formulation.     :cvar EVAL: direct pressure gradient evaluation for flux
51       :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*
52                   where *l* is the length scale, *K* is the inverse of the permeability tensor, and *w* is a positive weighting factor.
53       :cvar SMOOTH: global smoothing by solving the PDE *K_{ij} u_j= - p_{,j} + K_{ij} g_j*
54     """     """
55     def __init__(self, domain, useReduced=False, useVPIteration=True, *args, **kargs):     EVAL="EVAL"
56       SIMPLE="EVAL"
57       POST="POST"
58       SMOOTH="SMOOTH"
59       def __init__(self, domain, useReduced=False, solver="POST", verbose=False, w=1.):
60        """        """
61        initializes the Darcy flux problem        initializes the Darcy flux problem.
62    
63        :param domain: domain of the problem        :param domain: domain of the problem
64        :type domain: `Domain`        :type domain: `Domain`
65        :param useReduced: uses reduced oreder on flux and pressure        :param useReduced: uses reduced oreder on flux and pressure
66        :type useReduced: ``bool``        :type useReduced: ``bool``
67        :param adaptSubTolerance: switches on automatic subtolerance selection        :param solver: solver method
68        :param useVPIteration: if True altenative iteration over v and p is performed. Otherwise V and P are calculated in a single PDE.        :type solver: in [`DarcyFlow.EVAL`, `DarcyFlow.POST`, `DarcyFlow.SMOOTH` ]
69        :type useVPIteration: ``bool``            :param verbose: if ``True`` some information on the iteration progress are printed.
70          :type verbose: ``bool``
71          :param w: weighting factor for `DarcyFlow.POST` solver
72          :type w: ``float``
73          
74        """        """
75          if not solver in [DarcyFlow.EVAL, DarcyFlow.POST,  DarcyFlow.SMOOTH ] :
76              raise ValueError("unknown solver %d."%solver)
77    
78        self.domain=domain        self.domain=domain
79        self.useVPIteration=useVPIteration        self.solver=solver
80        self.useReduced=useReduced        self.useReduced=useReduced
81        self.verbose=False        self.verbose=verbose
82          self.l=None
83        if self.useVPIteration:        self.w=None
84       # this does not work yet      
85           self.__pde_k=LinearPDESystem(domain)        self.__pde_p=LinearSinglePDE(domain)
86           self.__pde_k.setSymmetryOn()        self.__pde_p.setSymmetryOn()
87           if self.useReduced: self.__pde_k.setReducedOrderOn()        if self.useReduced: self.__pde_p.setReducedOrderOn()
88      
89           self.__pde_p=LinearSinglePDE(domain)        if self.solver  == self.EVAL:
90           self.__pde_p.setSymmetryOn()           self.__pde_v=None
91           if self.useReduced: self.__pde_p.setReducedOrderOn()           if self.verbose: print("DarcyFlow: simple solver is used.")
92        else:  
93           self.__pde_k=LinearPDE(self.domain, numEquations=self.domain.getDim()+1)        elif self.solver  == self.POST:
94           self.__pde_k.setSymmetryOff()           if util.inf(w)<0.:
95                        raise ValueError("Weighting factor must be non-negative.")
96           if self.useReduced: self.__pde_k.setReducedOrderOn()           if self.verbose: print("DarcyFlow: global postprocessing of flux is used.")
97           C=self.__pde_k.createCoefficient("C")           self.__pde_v=LinearPDESystem(domain)
98           B=self.__pde_k.createCoefficient("B")           self.__pde_v.setSymmetryOn()
99           for i in range(self.domain.getDim()):           if self.useReduced: self.__pde_v.setReducedOrderOn()
100              B[i,i,self.domain.getDim()]=-1           self.w=w
101              C[self.domain.getDim(),i,i]=1           x=self.domain.getX()
102              C[i,self.domain.getDim(),i]=-0.5           self.l=min( [util.sup(x[i])-util.inf(x[i]) for i in xrange(self.domain.getDim()) ] )
103              B[self.domain.getDim(),i,i]=0.5           #self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale
104           self.__pde_k.setValue(C=C, B=B)  
105        self.__f=escript.Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))        elif self.solver  == self.SMOOTH:
106        self.__g=escript.Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))           self.__pde_v=LinearPDESystem(domain)
107        self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_k.getFunctionSpaceForCoefficient("q"))           self.__pde_v.setSymmetryOn()
108        self.location_of_fixed_flux = escript.Vector(0, self.__pde_k.getFunctionSpaceForCoefficient("q"))           if self.useReduced: self.__pde_v.setReducedOrderOn()
109        self.scale=1.           if self.verbose: print("DarcyFlow: flux smoothing is used.")
110        self.setTolerance()           self.w=0
111          
112          self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
113     def __L2(self,v):        self.__g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
114           return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))          self.__permeability_invXg=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
115     def __inner_GMRES(self,r,s):        self.__permeability_invXg_ref=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)
116           return util.integrate(util.inner(r,s))        self.ref_point_id=None
117                  self.ref_point=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)
118     def __Aprod_GMRES(self,p):        self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
119        self.__pde_k.setValue(Y=-0.5*util.grad(p), X=-p*util.kronecker(self.__pde_k.getDomain()) )        self.location_of_fixed_flux = escript.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
120        du=self.__pde_k.getSolution()        self.perm_scale=1.
121        self.__pde_p.setValue(Y=-util.div(du), X=0.5*(-du+util.tensor_mult(self.__permeability,util.grad(p))))      
122        return self.__pde_p.getSolution()          
           
    def getSolverOptionsFlux(self):  
       """  
       Returns the solver options used to solve the flux problems  
         
       *K^{-1} u=F*  
         
       :return: `SolverOptions`  
       """  
       return self.__pde_k.getSolverOptions()        
         
123     def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):     def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
124        """        """
125        assigns values to model parameters        assigns values to model parameters
# Line 129  class DarcyFlowVeryNew(object): Line 133  class DarcyFlowVeryNew(object):
133        :param location_of_fixed_flux:  mask for locations where flux is fixed.        :param location_of_fixed_flux:  mask for locations where flux is fixed.
134        :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)        :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
135        :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.        :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
136        :type permeability: scalar or tensor values on the domain (e.g. `escript.Data`)        :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
137    
138        :note: the values of parameters which are not set by calling ``setValue`` are not altered.        :note: the values of parameters which are not set by calling ``setValue`` are not altered.
139        :note: at any point on the boundary of the domain the pressure        :note: at any point on the boundary of the domain the pressure
# Line 138  class DarcyFlowVeryNew(object): Line 142  class DarcyFlowVeryNew(object):
142               is along the *x_i* axis.               is along the *x_i* axis.
143    
144        """        """
145        if location_of_fixed_pressure!=None: self.location_of_fixed_pressure=location_of_fixed_pressure        if location_of_fixed_pressure!=None:
146        if location_of_fixed_flux!=None: self.location_of_fixed_flux=location_of_fixed_flux             self.location_of_fixed_pressure=util.wherePositive(util.interpolate(location_of_fixed_pressure, self.__pde_p.getFunctionSpaceForCoefficient("q")))
147                     self.ref_point_id=self.location_of_fixed_pressure.maxGlobalDataPoint()
148        if self.useVPIteration:             if not self.location_of_fixed_pressure.getTupleForGlobalDataPoint(*self.ref_point_id)[0] > 0: raise ValueError("pressure needs to be fixed at least one point.")
149           if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=self.location_of_fixed_pressure)             self.ref_point=self.__pde_p.getFunctionSpaceForCoefficient("q").getX().getTupleForGlobalDataPoint(*self.ref_point_id)
150           if location_of_fixed_flux!=None: self.__pde_k.setValue(q=self.location_of_fixed_flux)             if self.verbose: print(("DarcyFlow: reference point at %s."%(self.ref_point,)))
151        else:             self.__pde_p.setValue(q=self.location_of_fixed_pressure)
152       if location_of_fixed_pressure!=None or location_of_fixed_flux!=None:        if location_of_fixed_flux!=None:
153          q=self.__pde_k.createCoefficient("q")            self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
154          q[self.domain.getDim()]=self.location_of_fixed_pressure            if not self.__pde_v == None:
155          q[:self.domain.getDim()]=self.location_of_fixed_flux                self.__pde_v.setValue(q=self.location_of_fixed_flux)
         self.__pde_k.setValue(q=q)  
156                            
       # flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain)  
157        if permeability!=None:        if permeability!=None:
158       perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A"))      
159           V=util.vol(self.domain)           perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
160       if perm.getRank()==0:           self.perm_scale=util.Lsup(util.length(perm))
161          perm_inv=(1./perm)           if self.verbose: print(("DarcyFlow: permeability scaling factor = %e."%self.perm_scale))
162              self.scale=util.integrate(perm_inv)*V**(1./self.domain.getDim()-1.)           perm=perm*(1./self.perm_scale)
         perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))  
         perm=perm*(self.scale*util.kronecker(self.domain.getDim()))  
      elif perm.getRank()==2:  
         perm_inv=util.inverse(perm)  
             self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)*V**(2./self.domain.getDim()-1.)/self.domain.getDim())  
         perm_inv*=(1./self.scale)  
         perm=perm*self.scale  
      else:  
         raise ValueError,"illegal rank of permeability."  
163                    
164       self.__permeability=perm           if perm.getRank()==0:
165       self.__permeability_inv=perm_inv  
166       if self.useVPIteration:              perm_inv=(1./perm)
167              self.__pde_k.setValue(D=0.5*self.__permeability_inv)              perm_inv=perm_inv*util.kronecker(self.domain.getDim())
168          self.__pde_p.setValue(A=0.5*self.__permeability)              perm=perm*util.kronecker(self.domain.getDim())
169            
170            
171             elif perm.getRank()==2:
172                perm_inv=util.inverse(perm)
173           else:           else:
174              D=self.__pde_k.createCoefficient("D")              raise ValueError("illegal rank of permeability.")
175              A=self.__pde_k.createCoefficient("A")          
176              D[:self.domain.getDim(),:self.domain.getDim()]=0.5*self.__permeability_inv           self.__permeability=perm
177              A[self.domain.getDim(),:,self.domain.getDim(),:]=0.5*self.__permeability           self.__permeability_inv=perm_inv
178              self.__pde_k.setValue(A=A, D=D)      
179             #====================
180             self.__pde_p.setValue(A=self.__permeability)
181             if self.solver  == self.EVAL:
182                  pass # no extra work required
183             elif self.solver  == self.POST:
184                  k=util.kronecker(self.domain.getDim())
185                  self.omega = self.w*util.length(perm_inv)*self.l*self.domain.getSize()
186                  #self.__pde_v.setValue(D=self.__permeability_inv, A=self.omega*util.outer(k,k))
187                  self.__pde_v.setValue(D=self.__permeability_inv, A_reduced=self.omega*util.outer(k,k))
188             elif self.solver  == self.SMOOTH:
189                self.__pde_v.setValue(D=self.__permeability_inv)
190    
191        if g != None:        if g != None:
192      g=util.interpolate(g, self.__pde_k.getFunctionSpaceForCoefficient("Y"))          g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
193      if g.isEmpty():          if g.isEmpty():
194            g=Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))               g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
195      else:          else:
196          if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g"               if not g.getShape()==(self.domain.getDim(),): raise ValueError("illegal shape of g")
197          self.__g=g          self.__g=g
198            self.__permeability_invXg=util.tensor_mult(self.__permeability_inv,self.__g * (1./self.perm_scale ))
199            self.__permeability_invXg_ref=util.integrate(self.__permeability_invXg)/util.vol(self.domain)
200        if f !=None:        if f !=None:
201       f=util.interpolate(f, self.__pde_k.getFunctionSpaceForCoefficient("X"))           f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
202       if f.isEmpty():           if f.isEmpty():      
203                       f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
204            f=Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))           else:
205       else:               if f.getRank()>0: raise ValueError("illegal rank of f.")
206           if f.getRank()>0: raise ValueError,"illegal rank of f."           self.__f=f
          self.__f=f  
           
    def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10, iter_restart=20):  
       """  
       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. `escript.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. `escript.Data`).  
       :param verbose: if set some information on iteration progress are printed  
       :type verbose: ``bool``  
       :return: flux and pressure  
       :rtype: ``tuple`` of `escript.Data`.  
   
       """  
       u0_b=u0*self.location_of_fixed_flux  
       p0_b=p0*self.location_of_fixed_pressure/self.scale  
       f=self.__f-util.div(u0_b)  
       g=self.__g-u0_b - util.tensor_mult(self.__permeability,util.grad(p0_b))  
       self.verbose=verbose  
       if self.useVPIteration:  
         # get u:  
             rtol=self.getTolerance()  
         self.__pde_k.setValue(Y=0.5*util.tensor_mult(self.__permeability_inv,g),X=escript.Data())  
         du=self.__pde_k.getSolution()  
         self.__pde_p.setValue(Y=f-util.div(du), X=0.5*(g-du))  
         p=GMRES(self.__pde_p.getSolution(),  
                 self.__Aprod_GMRES,  
             p0_b*0,  
             self.__inner_GMRES,  
             atol=0,  
             rtol=rtol,  
             iter_max=max_iter,  
             iter_restart=iter_restart, verbose=self.verbose,P_R=None)  
             self.__pde_k.setValue(Y=0.5*( util.tensor_mult(self.__permeability_inv,g) + util.grad(p)) ,  
                                   X=p*util.kronecker(self.__pde_k.getDomain()))  
             #self.__pde_k.setValue(Y=0.5*util.grad(p), X=p*util.kronecker(self.__pde_k.getDomain()) )  
         u=self.__pde_k.getSolution()  
       else:  
           X=self.__pde_k.createCoefficient("X")  
           Y=self.__pde_k.createCoefficient("Y")  
           Y[:self.domain.getDim()]=0.5*util.tensor_mult(self.__permeability_inv,g)  
           Y[self.domain.getDim()]=f  
           X[self.domain.getDim(),:]=g*0.5  
           self.__pde_k.setValue(X=X, Y=Y)  
           self.__pde_k.getSolverOptions().setVerbosity(self.verbose)  
           #self.__pde_k.getSolverOptions().setPreconditioner(self.__pde_k.getSolverOptions().AMG)  
           self.__pde_k.getSolverOptions().setSolverMethod(self.__pde_k.getSolverOptions().DIRECT)  
           U=self.__pde_k.getSolution()  
           u=U[:self.domain.getDim()]  
           p=U[self.domain.getDim()]  
       # self.__pde_k.getOperator().saveMM("k.mm")  
       u=u0_b+u  
       p=(p0_b+p)*self.scale  
       if self.verbose:  
         KGp=util.tensor_mult(self.__permeability,util.grad(p)/self.scale)  
         def_p=self.__g-(u+KGp)  
         def_v=self.__f-util.div(u, self.__pde_k.getFunctionSpaceForCoefficient("X"))  
         print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(u))  
         print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))  
       return u,p  
    def setTolerance(self,rtol=1e-4):  
       """  
       sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if  
207    
       *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*  
         
       where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`).  
         
       :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  
 class DarcyFlow(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, useReduced=False, adaptSubTolerance=True, solveForFlux=False, useVPIteration=True, weighting_scale=0.1):  
       """  
       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``  
       :param solveForFlux: if True the solver solves for the flux (do not use!)  
       :type solveForFlux: ``bool``    
       :param useVPIteration: if True altenative iteration over v and p is performed. Otherwise V and P are calculated in a single PDE.  
       :type useVPIteration: ``bool``      
       """  
       self.domain=domain  
       self.useVPIteration=useVPIteration  
       self.useReduced=useReduced  
       self.weighting_scale=weighting_scale  
       if self.useVPIteration:  
          self.solveForFlux=solveForFlux  
          self.__adaptSubTolerance=adaptSubTolerance  
          self.verbose=False  
           
          self.__pde_k=LinearPDESystem(domain)  
          self.__pde_k.setSymmetryOn()  
          if self.useReduced: self.__pde_k.setReducedOrderOn()  
     
          self.__pde_p=LinearSinglePDE(domain)  
          self.__pde_p.setSymmetryOn()  
          if self.useReduced: self.__pde_p.setReducedOrderOn()  
          self.setTolerance()  
          self.setAbsoluteTolerance()  
       else:  
          self.__pde_k=LinearPDE(self.domain, numEquations=self.domain.getDim()+1)  
          self.__pde_k.setSymmetryOn()  
          if self.useReduced: self.__pde_k.setReducedOrderOn()  
          C=self.__pde_k.createCoefficient("C")  
          B=self.__pde_k.createCoefficient("B")  
          for i in range(self.domain.getDim()):  
             C[i,self.domain.getDim(),i]=1  
             B[self.domain.getDim(),i,i]=1  
          self.__pde_k.setValue(C=C, B=B)  
       self.__f=escript.Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))  
       self.__g=escript.Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))  
         
208     def getSolverOptionsFlux(self):     def getSolverOptionsFlux(self):
209        """        """
210        Returns the solver options used to solve the flux problems        Returns the solver options used to solve the flux problems
         
       *K^{-1} u=F*  
         
211        :return: `SolverOptions`        :return: `SolverOptions`
212        """        """
213        return self.__pde_k.getSolverOptions()        if self.__pde_v == None:
214              return None
215          else:
216              return self.__pde_v.getSolverOptions()
217                
218     def setSolverOptionsFlux(self, options=None):     def setSolverOptionsFlux(self, options=None):
219        """        """
220        Sets the solver options used to solve the flux problems        Sets the solver options used to solve the flux problems
         
       *K^{-1}u=F*  
         
221        If ``options`` is not present, the options are reset to default        If ``options`` is not present, the options are reset to default
         
222        :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.  
223        """        """
224        return self.__pde_v.setSolverOptions(options)        if not self.__pde_v == None:
225              self.__pde_v.setSolverOptions(options)
226            
227     def getSolverOptionsPressure(self):     def getSolverOptionsPressure(self):
228        """        """
229        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*  
         
230        :return: `SolverOptions`        :return: `SolverOptions`
231        """        """
232        return self.__pde_p.getSolverOptions()        return self.__pde_p.getSolverOptions()
# Line 371  class DarcyFlow(object): Line 234  class DarcyFlow(object):
234     def setSolverOptionsPressure(self, options=None):     def setSolverOptionsPressure(self, options=None):
235        """        """
236        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*  
         
237        If ``options`` is not present, the options are reset to default        If ``options`` is not present, the options are reset to default
238                
239        :param options: `SolverOptions`        :param options: `SolverOptions`
240        :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.
241        """        """
242        return self.__pde_p.setSolverOptions(options)        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. `escript.Data`)  
       :param g: flux sources/sinks  
       :type g: vector values on the domain (e.g. `escript.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. `escript.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. `escript.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. `escript.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 self.useVPIteration:  
          if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)  
          if location_of_fixed_flux!=None: self.__pde_k.setValue(q=location_of_fixed_flux)  
       else:  
          q=self.__pde_k.getCoefficient("q")  
          if q.isEmpty(): q=self.__pde_k.createCoefficient("q")  
          if location_of_fixed_pressure!=None: q[self.domain.getDim()]=location_of_fixed_pressure  
          if location_of_fixed_flux!=None: q[:self.domain.getDim()]=location_of_fixed_flux  
          self.__pde_k.setValue(q=q)  
               
       # flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain)  
       if permeability!=None:  
      perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A"))  
          V=util.vol(self.domain)  
      if perm.getRank()==0:  
         perm_inv=(1./perm)  
             if self.useVPIteration:  
               self.scale=1.  
             else:  
               self.scale=util.integrate(perm_inv)*V**(1./self.domain.getDim()-1.)  
   
         perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))  
         perm=perm*(self.scale*util.kronecker(self.domain.getDim()))  
      elif perm.getRank()==2:  
         perm_inv=util.inverse(perm)  
             if self.useVPIteration:  
               self.scale=1.  
             else:  
               self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)*V**(2./self.domain.getDim()-1.)/self.domain.getDim())  
           perm_inv*=(1./self.scale)  
           perm=perm*self.scale  
      else:  
         raise ValueError,"illegal rank of permeability."  
   
      self.__permeability=perm  
      self.__permeability_inv=perm_inv  
   
      self.__l2 =(util.longestEdge(self.domain)**2*util.length(self.__permeability_inv))*self.weighting_scale  
          if self.useVPIteration:  
         if  self.solveForFlux:  
            self.__pde_k.setValue(D=self.__permeability_inv)  
         else:  
            self.__pde_k.setValue(D=self.__permeability_inv, A=self.__l2*util.outer(util.kronecker(self.domain),util.kronecker(self.domain)))  
         self.__pde_p.setValue(A=self.__permeability)  
          else:  
             D=self.__pde_k.createCoefficient("D")  
             A=self.__pde_k.createCoefficient("A")  
             D[:self.domain.getDim(),:self.domain.getDim()]=self.__permeability_inv  
             for i in range(self.domain.getDim()):  
                for j in range(self.domain.getDim()):  
                  A[i,i,j,j]=self.__l2  
             A[self.domain.getDim(),:,self.domain.getDim(),:]=self.__permeability  
             self.__pde_k.setValue(A=A, D=D)  
       if g !=None:  
      g=util.interpolate(g, self.__pde_k.getFunctionSpaceForCoefficient("Y"))  
      if g.isEmpty():  
           g=Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))  
      else:  
         if not g.getShape()==(self.domain.getDim(),):  
               raise ValueError,"illegal shape of g"  
         self.__g=g  
       elif permeability!=None:  
              X  
       if f !=None:  
      f=util.interpolate(f, self.__pde_k.getFunctionSpaceForCoefficient("X"))  
      if f.isEmpty():  
           f=Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))  
      else:  
          if f.getRank()>0: raise ValueError,"illegal rank of f."  
          self.__f=f  
   
         
    def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):  
       """  
       solves the problem.  
243                
244        The iteration is terminated if the residual norm is less then self.getTolerance().     def solve(self, u0, p0):
   
       :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. `escript.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. `escript.Data`).  
       :param verbose: if set some information on iteration progress are printed  
       :type verbose: ``bool``  
       :return: flux and pressure  
       :rtype: ``tuple`` of `escript.Data`.  
   
       :note: The problem is solved as a least squares form  
              *(K^[-1]+D^* l2 D)u+G p=D^* l2 * 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  
       if self.useVPIteration:  
           return self.__solveVP(u0,p0,max_iter,max_num_corrections)  
       else:  
           X=self.__pde_k.createCoefficient("X")  
           Y=self.__pde_k.createCoefficient("Y")  
           Y[:self.domain.getDim()]=self.scale*util.tensor_mult(self.__permeability_inv,self.__g)  
           rtmp=self.__f * self.__l2 * self.scale  
           for i in range(self.domain.getDim()): X[i,i]=rtmp  
           X[self.domain.getDim(),:]=self.__g*self.scale  
           r=self.__pde_k.createCoefficient("r")  
           r[:self.domain.getDim()]=u0*self.scale  
           r[self.domain.getDim()]=p0  
           self.__pde_k.setValue(X=X, Y=Y, r=r)  
           self.__pde_k.getSolverOptions().setVerbosity(self.verbose)  
           #self.__pde_k.getSolverOptions().setPreconditioner(self.__pde_k.getSolverOptions().AMG)  
           self.__pde_k.getSolverOptions().setSolverMethod(self.__pde_k.getSolverOptions().DIRECT)  
           U=self.__pde_k.getSolution()  
           # self.__pde_k.getOperator().saveMM("k.mm")  
           u=U[:self.domain.getDim()]*(1./self.scale)  
           p=U[self.domain.getDim()]  
           if self.verbose:  
         KGp=util.tensor_mult(self.__permeability,util.grad(p)/self.scale)  
         def_p=self.__g-(u+KGp)  
         def_v=self.__f-util.div(u, self.__pde_k.getFunctionSpaceForCoefficient("X"))  
         print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(u))  
         print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))  
           return u,p  
   
    def __solveVP(self,u0,p0, max_iter=100, max_num_corrections=10):  
245        """        """
246        solves the problem.        solves the problem.
247                
       The iteration is terminated if the residual norm is less than self.getTolerance().  
   
248        :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.
249        :type u0: vector value on the domain (e.g. `escript.Data`).        :type u0: vector value on the domain (e.g. `escript.Data`).
250        :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.
# Line 538  class DarcyFlow(object): Line 252  class DarcyFlow(object):
252        :return: flux and pressure        :return: flux and pressure
253        :rtype: ``tuple`` of `escript.Data`.        :rtype: ``tuple`` of `escript.Data`.
254    
       :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}*.  
255        """        """
256        rtol=self.getTolerance()        p0=util.interpolate(p0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
257        atol=self.getAbsoluteTolerance()        if self.ref_point_id == None:
258        self.setSubProblemTolerance()            p_ref=0
       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  
   
    def __applWeight(self, v, f=None):  
       # solves L p = f-Dv with p = 0  
       if self.verbose: print "DarcyFlux: Applying weighting operator"  
       if f == None:  
      return -util.div(v)*self.__l2  
       else:  
      return (f-util.div(v))*self.__l2  
    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)  
259        else:        else:
260       self.__pde_p.setValue(X=g-v, r=p0)                  p_ref=p0.getTupleForGlobalDataPoint(*self.ref_point_id)[0]
261        p=self.__pde_p.getSolution()        p0_hydrostatic=p_ref+util.inner(self.__permeability_invXg_ref, self.__pde_p.getFunctionSpaceForCoefficient("q").getX() - self.ref_point)
262        return p        g_2=self.__g - util.tensor_mult(self.__permeability, self.__permeability_invXg_ref * self.perm_scale)
263          self.__pde_p.setValue(X=g_2 * 1./self.perm_scale,
264     def __Aprod_v(self,dv):                              Y=self.__f * 1./self.perm_scale,
265        # calculates: (a,b,c) = (K^{-1}(dv + KG * dp), L^{-1}Ddv, dp)  with (G*KG)dp = - G^*dv                                y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux * 1./self.perm_scale ),
266        dp=self.__getPressure(dv, p0=escript.Data()) # dp = (G*KG)^{-1} (0-G^*dv)                              r=p0 - p0_hydrostatic)
267        a=util.tensor_mult(self.__permeability_inv,dv)+util.grad(dp) # a= K^{-1}u+G*dp        pp=self.__pde_p.getSolution()
268        b= - self.__applWeight(dv) # b = - (D K D^*)^{-1} (0-Dv)        u = self._getFlux(pp, u0)
269        return ArithmeticTuple(a,b,-dp)        return u,pp + p0_hydrostatic
   
    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=escript.Data())  
       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])  
270                
271     def __Aprod_p(self,dp):     def getFlux(self,p, u0=None):
       if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"  
       Gdp=util.grad(dp)  
       self.__pde_k.setValue(Y=-Gdp,X=escript.Data(), r=escript.Data())  
       du=self.__pde_k.getSolution()  
       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=escript.Data(), r=escript.Data(), y=escript.Data())  
       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,escript.Function(self.domain)))**2))  
   
    def __L2_r(self,v):  
       return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.ReducedFunction(self.domain)))**2))  
   
    def setTolerance(self,rtol=1e-4):  
       """  
       sets the relative tolerance ``rtol`` used to terminate the solution process. The iteration is terminated if  
   
       *|g-v-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*  
         
       where ``atol`` is an absolut tolerance (see `setAbsoluteTolerance`).  
         
       :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-K gard(p)|_PCG <= atol + rtol * |K^{1/2}g2|_0*  
   
   
       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.)  
      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):  
272          """          """
273          initializes the Darcy flux problem          returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
         :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.__l2=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2  
            # self.__l2=(3.*util.longestEdge(self.domain))**2  
            #self.__l2=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2  
         else:  
            self.__l2=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.__l2*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=escript.Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))  
         self.__g=escript.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):  
         """  
         assigns values to model parameters  
   
         :param f: volumetic sources/sinks  
         :type f: scalar value on the domain (e.g. `escript.Data`)  
         :param g: flux sources/sinks  
         :type g: vector values on the domain (e.g. `escript.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. `escript.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. `escript.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. `escript.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=escript.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=escript.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. `escript.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. `escript.Data`).  
          :param verbose: if set some information on iteration progress are printed  
          :type verbose: ``bool``  
          :return: flux and pressure  
          :rtype: ``tuple`` of `escript.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,escript.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,escript.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,escript.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=escript.Data(), r=escript.Data())  
           du=self.__pde_v.getSolution()  
           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=escript.Data(), r=escript.Data())  
           return self.__pde_p.getSolution()  
   
     def getFlux(self,p=None, fixed_flux=escript.Data()):  
         """  
         returns the flux for a given pressure ``p`` where the flux is equal to ``fixed_flux``  
274          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
275          Note that ``g`` and ``f`` are used, see `setValue`.          Notice that ``g`` is used, see `setValue`.
276    
277          :param p: pressure.          :param p: pressure.
278          :type p: scalar value on the domain (e.g. `escript.Data`).          :type p: scalar value on the domain (e.g. `escript.Data`).
279          :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``.
280          :type fixed_flux: vector values on the domain (e.g. `escript.Data`).          :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
281          :return: flux          :return: flux
282          :rtype: `escript.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}*  
283          """          """
284      self.setSubProblemTolerance()          p=util.interpolate(p, self.__pde_p.getFunctionSpaceForCoefficient("q"))
285          g=self.__g          if self.ref_point_id == None:
286          f=self.__f              p_ref=0
         self.__pde_v.setValue(X=self.__l2*f*util.kronecker(self.domain), r=fixed_flux)  
         if p == None:  
            self.__pde_v.setValue(Y=g)  
287          else:          else:
288             self.__pde_v.setValue(Y=g-self.__Q(p))              p_ref=p.getTupleForGlobalDataPoint(*self.ref_point_id)[0]
289          return self.__pde_v.getSolution()          p_hydrostatic=p_ref+util.inner(self.__permeability_invXg_ref, self.__pde_p.getFunctionSpaceForCoefficient("q").getX() - self.ref_point)
290            return self._getFlux(p-p_hydrostatic, u0)
291    
292       def _getFlux(self, pp, u0=None):
293            """
294            returns the flux for a given pressure ``pp`` where the flux is equal to
295            ``u0`` on locations where ``location_of_fixed_flux`` is positive (see
296            `setValue`). Notice that ``g`` is used, see `setValue`.
297    
298            :param pp: pressure.
299            :type pp: scalar value on the domain (i.e. `escript.Data`).
300            :param u0: flux on the locations of the domain marked in ``location_of_fixed_flux``.
301            :type u0: vector values on the domain (i.e. `escript.Data`) or ``None``
302            :return: flux
303            :rtype: `escript.Data`
304            """
305            if self.solver  == self.EVAL:
306               u = self.__g - util.tensor_mult(self.__permeability, self.perm_scale * (util.grad(pp) + self.__permeability_invXg_ref))
307            elif self.solver  == self.POST or self.solver  == self.SMOOTH:
308                self.__pde_v.setValue(Y= self.__permeability_invXg - (util.grad(pp) + self.__permeability_invXg_ref))
309                print
310                if u0 == None:
311                   self.__pde_v.setValue(r=escript.Data())
312                else:
313                   if not isinstance(u0, escript.Data) : u0 = escript.Vector(u0, escript.Solution(self.domain))
314                   self.__pde_v.setValue(r=1./self.perm_scale * u0)
315                u= self.__pde_v.getSolution() * self.perm_scale
316            return u
317          
318  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
319       """       """
320       solves       solves
# Line 1074  class StokesProblemCartesian(Homogeneous Line 333  class StokesProblemCartesian(Homogeneous
333              sp.setTolerance()              sp.setTolerance()
334              sp.initialize(...)              sp.initialize(...)
335              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
336                sp.setStokesEquation(...) # new values for some parameters
337                v1,p1=sp.solve(v,p)
338       """       """
339       def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
340           """           """
# Line 1088  class StokesProblemCartesian(Homogeneous Line 349  class StokesProblemCartesian(Homogeneous
349           """           """
350           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
351           self.domain=domain           self.domain=domain
352           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())
353           self.__pde_u.setSymmetryOn()           self.__pde_v.setSymmetryOn()
354            
355           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
356           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
# Line 1097  class StokesProblemCartesian(Homogeneous Line 358  class StokesProblemCartesian(Homogeneous
358    
359           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
360           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
361       self.__pde_proj.setValue(D=1)           self.__pde_proj.setValue(D=1)
362           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
363    
364       def getSolverOptionsVelocity(self):       def getSolverOptionsVelocity(self):
# Line 1106  class StokesProblemCartesian(Homogeneous Line 367  class StokesProblemCartesian(Homogeneous
367            
368       :rtype: `SolverOptions`       :rtype: `SolverOptions`
369       """       """
370       return self.__pde_u.getSolverOptions()           return self.__pde_v.getSolverOptions()
371       def setSolverOptionsVelocity(self, options=None):       def setSolverOptionsVelocity(self, options=None):
372           """           """
373       set the solver options for solving the equation for velocity.       set the solver options for solving the equation for velocity.
# Line 1114  class StokesProblemCartesian(Homogeneous Line 375  class StokesProblemCartesian(Homogeneous
375       :param options: new solver  options       :param options: new solver  options
376       :type options: `SolverOptions`       :type options: `SolverOptions`
377       """       """
378           self.__pde_u.setSolverOptions(options)           self.__pde_v.setSolverOptions(options)
379       def getSolverOptionsPressure(self):       def getSolverOptionsPressure(self):
380           """           """
381       returns the solver options used  solve the equation for pressure.       returns the solver options used  solve the equation for pressure.
382       :rtype: `SolverOptions`       :rtype: `SolverOptions`
383       """       """
384       return self.__pde_prec.getSolverOptions()           return self.__pde_prec.getSolverOptions()
385       def setSolverOptionsPressure(self, options=None):       def setSolverOptionsPressure(self, options=None):
386           """           """
387       set the solver options for solving the equation for pressure.       set the solver options for solving the equation for pressure.
388       :param options: new solver  options       :param options: new solver  options
389       :type options: `SolverOptions`       :type options: `SolverOptions`
390       """       """
391       self.__pde_prec.setSolverOptions(options)           self.__pde_prec.setSolverOptions(options)
392    
393       def setSolverOptionsDiv(self, options=None):       def setSolverOptionsDiv(self, options=None):
394           """           """
# Line 1137  class StokesProblemCartesian(Homogeneous Line 398  class StokesProblemCartesian(Homogeneous
398       :param options: new solver options       :param options: new solver options
399       :type options: `SolverOptions`       :type options: `SolverOptions`
400       """       """
401       self.__pde_proj.setSolverOptions(options)           self.__pde_proj.setSolverOptions(options)
402       def getSolverOptionsDiv(self):       def getSolverOptionsDiv(self):
403           """           """
404       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 1145  class StokesProblemCartesian(Homogeneous Line 406  class StokesProblemCartesian(Homogeneous
406            
407       :rtype: `SolverOptions`       :rtype: `SolverOptions`
408       """       """
409       return self.__pde_proj.getSolverOptions()           return self.__pde_proj.getSolverOptions()
410    
411       def updateStokesEquation(self, v, p):       def updateStokesEquation(self, v, p):
412           """           """
413           updates the Stokes equation to consider dependencies from ``v`` and ``p``           updates the Stokes equation to consider dependencies from ``v`` and ``p``
414           :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.
415           """           """
416           pass           pass
417       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 1172  class StokesProblemCartesian(Homogeneous Line 433  class StokesProblemCartesian(Homogeneous
433              k=util.kronecker(self.domain.getDim())              k=util.kronecker(self.domain.getDim())
434              kk=util.outer(k,k)              kk=util.outer(k,k)
435              self.eta=util.interpolate(eta, escript.Function(self.domain))              self.eta=util.interpolate(eta, escript.Function(self.domain))
436          self.__pde_prec.setValue(D=1/self.eta)              self.__pde_prec.setValue(D=1/self.eta)
437              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)))
438          if restoration_factor!=None:          if restoration_factor!=None:
439              n=self.domain.getNormal()              n=self.domain.getNormal()
440              self.__pde_u.setValue(d=restoration_factor*util.outer(n,n))              self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
441          if fixed_u_mask!=None:          if fixed_u_mask!=None:
442              self.__pde_u.setValue(q=fixed_u_mask)              self.__pde_v.setValue(q=fixed_u_mask)
443          if f!=None: self.__f=f          if f!=None: self.__f=f
444          if surface_stress!=None: self.__surface_stress=surface_stress          if surface_stress!=None: self.__surface_stress=surface_stress
445          if stress!=None: self.__stress=stress          if stress!=None: self.__stress=stress
# Line 1196  class StokesProblemCartesian(Homogeneous Line 457  class StokesProblemCartesian(Homogeneous
457          :param surface_stress: normal surface stress          :param surface_stress: normal surface stress
458          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
459          :param stress: initial stress          :param stress: initial stress
460      :type stress: `Tensor` object on `FunctionSpace` `Function` or similar          :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
461          """          """
462          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
463    
# Line 1209  class StokesProblemCartesian(Homogeneous Line 470  class StokesProblemCartesian(Homogeneous
470           :rtype: ``float``           :rtype: ``float``
471           """           """
472           self.__pde_proj.setValue(Y=-util.div(v))           self.__pde_proj.setValue(Y=-util.div(v))
473       self.getSolverOptionsDiv().setTolerance(tol)           self.getSolverOptionsDiv().setTolerance(tol)
474       self.getSolverOptionsDiv().setAbsoluteTolerance(0.)           self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
475           out=self.__pde_proj.getSolution()           out=self.__pde_proj.getSolution()
476           return out           return out
477    
# Line 1251  class StokesProblemCartesian(Homogeneous Line 512  class StokesProblemCartesian(Homogeneous
512    
513       def getDV(self, p, v, tol):       def getDV(self, p, v, tol):
514           """           """
515           return the value for v for a given p (overwrite)           return the value for v for a given p
516    
517           :param p: a pressure           :param p: a pressure
518           :param v: a initial guess for the value v to return.           :param v: a initial guess for the value v to return.
519           :return: dv given as *Adv=(f-Av-B^*p)*           :return: dv given as *Adv=(f-Av-B^*p)*
520           """           """
521           self.updateStokesEquation(v,p)           self.updateStokesEquation(v,p)
522           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
523       self.getSolverOptionsVelocity().setTolerance(tol)           self.getSolverOptionsVelocity().setTolerance(tol)
524       self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)           self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
525           if self.__stress.isEmpty():           if self.__stress.isEmpty():
526              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)))
527           else:           else:
528              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)))
529           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
530           return  out           return  out
531    
532       def norm_Bv(self,Bv):       def norm_Bv(self,Bv):
# Line 1285  class StokesProblemCartesian(Homogeneous Line 546  class StokesProblemCartesian(Homogeneous
546           :return: the solution of *Av=B^*p*           :return: the solution of *Av=B^*p*
547           :note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
548           """           """
549           self.__pde_u.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))           self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
550           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
551           return  out           return  out
552    
553       def solve_prec(self,Bv, tol):       def solve_prec(self,Bv, tol):
554           """           """
555           applies preconditioner for for *BA^{-1}B^** to *Bv*           applies preconditioner for for *BA^{-1}B^** to *Bv*
556           with accuracy `self.getSubProblemTolerance()`           with accuracy ``self.getSubProblemTolerance()``
557    
558           :param Bv: velocity increment           :param Bv: velocity increment
559           :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*           :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
560           :note: boundary conditions on p are zero.           :note: boundary conditions on p are zero.
561           """           """
562           self.__pde_prec.setValue(Y=Bv)           self.__pde_prec.setValue(Y=Bv)
563       self.getSolverOptionsPressure().setTolerance(tol)           self.getSolverOptionsPressure().setTolerance(tol)
564       self.getSolverOptionsPressure().setAbsoluteTolerance(0.)           self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
565           out=self.__pde_prec.getSolution()           out=self.__pde_prec.getSolution()
566           return out           return out

Legend:
Removed from v.3500  
changed lines
  Added in v.3990

  ViewVC Help
Powered by ViewVC 1.1.26