/[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 3499 by gross, Fri Apr 8 00:14:35 2011 UTC revision 4154 by jfenwick, Tue Jan 22 09:30:23 2013 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-2013 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-2013 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 escriptcpp
37  import util  escore=escriptcpp
38  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions  #from . import escript
39  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES  from . import util
40    from .linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
41    from .pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
42    
43  class DarcyFlowVeryNew(object):  class DarcyFlow(object):
44     """     """
45     solves the problem     solves the problem
46        
# Line 47  class DarcyFlowVeryNew(object): Line 49  class DarcyFlowVeryNew(object):
49        
50     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,
51        
52     :note: The problem is solved in a stabelized formulation.     :cvar EVAL: direct pressure gradient evaluation for flux
53       :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*
54                   where *l* is the length scale, *K* is the inverse of the permeability tensor, and *w* is a positive weighting factor.
55       :cvar SMOOTH: global smoothing by solving the PDE *K_{ij} u_j= - p_{,j} + K_{ij} g_j*
56     """     """
57     def __init__(self, domain, useReduced=False, *args, **kargs):     EVAL="EVAL"
58       SIMPLE="EVAL"
59       POST="POST"
60       SMOOTH="SMOOTH"
61       def __init__(self, domain, useReduced=False, solver="POST", verbose=False, w=1.):
62        """        """
63        initializes the Darcy flux problem        initializes the Darcy flux problem.
64    
65        :param domain: domain of the problem        :param domain: domain of the problem
66        :type domain: `Domain`        :type domain: `Domain`
67        :param useReduced: uses reduced oreder on flux and pressure        :param useReduced: uses reduced oreder on flux and pressure
68        :type useReduced: ``bool``        :type useReduced: ``bool``
69        :param adaptSubTolerance: switches on automatic subtolerance selection        :param solver: solver method
70        :type adaptSubTolerance: ``bool``        :type solver: in [`DarcyFlow.EVAL`, `DarcyFlow.POST`, `DarcyFlow.SMOOTH` ]
71        :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.
72        :type solveForFlux: ``bool``          :type verbose: ``bool``
73        :param useVPIteration: if True altenative iteration over v and p is performed. Otherwise V and P are calculated in a single PDE.        :param w: weighting factor for `DarcyFlow.POST` solver
74        :type useVPIteration: ``bool``            :type w: ``float``
75          
76        """        """
77          if not solver in [DarcyFlow.EVAL, DarcyFlow.POST,  DarcyFlow.SMOOTH ] :
78              raise ValueError("unknown solver %d."%solver)
79    
80        self.domain=domain        self.domain=domain
81        useVPIteration=False        self.solver=solver
       self.useVPIteration=useVPIteration or True  
       #self.useVPIteration=False  
82        self.useReduced=useReduced        self.useReduced=useReduced
83        self.verbose=False        self.verbose=verbose
84          self.l=None
85        if self.useVPIteration:        self.w=None
86       # this does not work yet      
87           self.__pde_k=LinearPDESystem(domain)        self.__pde_p=LinearSinglePDE(domain)
88           self.__pde_k.setSymmetryOn()        self.__pde_p.setSymmetryOn()
89           if self.useReduced: self.__pde_k.setReducedOrderOn()        if self.useReduced: self.__pde_p.setReducedOrderOn()
90      
91           self.__pde_p=LinearSinglePDE(domain)        if self.solver  == self.EVAL:
92           self.__pde_p.setSymmetryOn()           self.__pde_v=None
93           if self.useReduced: self.__pde_p.setReducedOrderOn()           if self.verbose: print("DarcyFlow: simple solver is used.")
94        else:  
95           self.__pde_k=LinearPDE(self.domain, numEquations=self.domain.getDim()+1)        elif self.solver  == self.POST:
96           self.__pde_k.setSymmetryOff()           if util.inf(w)<0.:
97                        raise ValueError("Weighting factor must be non-negative.")
98           if self.useReduced: self.__pde_k.setReducedOrderOn()           if self.verbose: print("DarcyFlow: global postprocessing of flux is used.")
99           C=self.__pde_k.createCoefficient("C")           self.__pde_v=LinearPDESystem(domain)
100           B=self.__pde_k.createCoefficient("B")           self.__pde_v.setSymmetryOn()
101           for i in range(self.domain.getDim()):           if self.useReduced: self.__pde_v.setReducedOrderOn()
102              B[i,i,self.domain.getDim()]=-1           self.w=w
103              C[self.domain.getDim(),i,i]=1           x=self.domain.getX()
104              C[i,self.domain.getDim(),i]=-0.5           self.l=min( [util.sup(x[i])-util.inf(x[i]) for i in range(self.domain.getDim()) ] )
105              B[self.domain.getDim(),i,i]=0.5           #self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale
106           self.__pde_k.setValue(C=C, B=B)  
107        self.__f=escript.Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))        elif self.solver  == self.SMOOTH:
108        self.__g=escript.Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))           self.__pde_v=LinearPDESystem(domain)
109        self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_k.getFunctionSpaceForCoefficient("q"))           self.__pde_v.setSymmetryOn()
110        self.location_of_fixed_flux = escript.Vector(0, self.__pde_k.getFunctionSpaceForCoefficient("q"))           if self.useReduced: self.__pde_v.setReducedOrderOn()
111        self.scale=1.           if self.verbose: print("DarcyFlow: flux smoothing is used.")
112     def __L2(self,v):           self.w=0
113           return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))    
114     def __inner_GMRES(self,r,s):        self.__f=escore.Data(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
115           return util.integrate(util.inner(r,s))        self.__g=escore.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
116                  self.__permeability_invXg=escore.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
117     def __Aprod_GMRES(self,p):        self.__permeability_invXg_ref=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)
118        self.__pde_k.setValue(Y=0.5*util.grad(p), X=p*util.kronecker(self.__pde_k.getDomain()) )        self.ref_point_id=None
119        du=self.__pde_k.getSolution()        self.ref_point=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)
120        self.__pde_p.setValue(Y=util.div(du), X=0.5*(du+util.tensor_mult(self.__permeability,util.grad(p))))        self.location_of_fixed_pressure = escore.Data(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
121        return self.__pde_p.getSolution()        self.location_of_fixed_flux = escore.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
122                  self.perm_scale=1.
123     def getSolverOptionsFlux(self):      
124        """          
       Returns the solver options used to solve the flux problems  
         
       *K^{-1} u=F*  
         
       :return: `SolverOptions`  
       """  
       return self.__pde_k.getSolverOptions()        
         
125     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):
126        """        """
127        assigns values to model parameters        assigns values to model parameters
# Line 131  class DarcyFlowVeryNew(object): Line 135  class DarcyFlowVeryNew(object):
135        :param location_of_fixed_flux:  mask for locations where flux is fixed.        :param location_of_fixed_flux:  mask for locations where flux is fixed.
136        :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`)
137        :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.
138        :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`)
139    
140        :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.
141        :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 140  class DarcyFlowVeryNew(object): Line 144  class DarcyFlowVeryNew(object):
144               is along the *x_i* axis.               is along the *x_i* axis.
145    
146        """        """
147        if location_of_fixed_pressure!=None: self.location_of_fixed_pressure=location_of_fixed_pressure        if location_of_fixed_pressure!=None:
148        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")))
149                     self.ref_point_id=self.location_of_fixed_pressure.maxGlobalDataPoint()
150        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.")
151           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)
152           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,)))
153        else:             self.__pde_p.setValue(q=self.location_of_fixed_pressure)
154       if location_of_fixed_pressure!=None or location_of_fixed_flux!=None:        if location_of_fixed_flux!=None:
155          q=self.__pde_k.createCoefficient("q")            self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
156          q[self.domain.getDim()]=self.location_of_fixed_pressure            if not self.__pde_v == None:
157          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)  
158                            
       # flux is rescaled by the factor mean value(perm_inv)*length where length**self.domain.getDim()=vol(self.domain)  
159        if permeability!=None:        if permeability!=None:
160       perm=util.interpolate(permeability,self.__pde_k.getFunctionSpaceForCoefficient("A"))      
161           V=util.vol(self.domain)           perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
162       if perm.getRank()==0:           self.perm_scale=util.Lsup(util.length(perm))
163          perm_inv=(1./perm)           if self.verbose: print(("DarcyFlow: permeability scaling factor = %e."%self.perm_scale))
164              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."  
165                    
166       self.__permeability=perm           if perm.getRank()==0:
167       self.__permeability_inv=perm_inv  
168       if self.useVPIteration:              perm_inv=(1./perm)
169              self.__pde_k.setValue(D=0.5*self.__permeability_inv)              perm_inv=perm_inv*util.kronecker(self.domain.getDim())
170          self.__pde_p.setValue(A=0.5*self.__permeability)              perm=perm*util.kronecker(self.domain.getDim())
171            
172            
173             elif perm.getRank()==2:
174                perm_inv=util.inverse(perm)
175           else:           else:
176              D=self.__pde_k.createCoefficient("D")              raise ValueError("illegal rank of permeability.")
177              A=self.__pde_k.createCoefficient("A")          
178              D[:self.domain.getDim(),:self.domain.getDim()]=0.5*self.__permeability_inv           self.__permeability=perm
179              A[self.domain.getDim(),:,self.domain.getDim(),:]=0.5*self.__permeability           self.__permeability_inv=perm_inv
180              self.__pde_k.setValue(A=A, D=D)      
181             #====================
182             self.__pde_p.setValue(A=self.__permeability)
183             if self.solver  == self.EVAL:
184                  pass # no extra work required
185             elif self.solver  == self.POST:
186                  k=util.kronecker(self.domain.getDim())
187                  self.omega = self.w*util.length(perm_inv)*self.l*self.domain.getSize()
188                  #self.__pde_v.setValue(D=self.__permeability_inv, A=self.omega*util.outer(k,k))
189                  self.__pde_v.setValue(D=self.__permeability_inv, A_reduced=self.omega*util.outer(k,k))
190             elif self.solver  == self.SMOOTH:
191                self.__pde_v.setValue(D=self.__permeability_inv)
192    
193        if g != None:        if g != None:
194      g=util.interpolate(g, self.__pde_k.getFunctionSpaceForCoefficient("Y"))          g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
195      if g.isEmpty():          if g.isEmpty():
196            g=Vector(0,self.__pde_k.getFunctionSpaceForCoefficient("Y"))               g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
197      else:          else:
198          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")
199          self.__g=g          self.__g=g
200            self.__permeability_invXg=util.tensor_mult(self.__permeability_inv,self.__g * (1./self.perm_scale ))
201            self.__permeability_invXg_ref=util.integrate(self.__permeability_invXg)/util.vol(self.domain)
202        if f !=None:        if f !=None:
203       f=util.interpolate(f, self.__pde_k.getFunctionSpaceForCoefficient("X"))           f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
204       if f.isEmpty():           if f.isEmpty():      
205                       f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
206            f=Scalar(0,self.__pde_k.getFunctionSpaceForCoefficient("X"))           else:
207       else:               if f.getRank()>0: raise ValueError("illegal rank of f.")
208           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):  
       """  
       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:  
         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=1.e-4,  
             iter_max=100,  
             iter_restart=20, 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()))  
         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  
209    
       *|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"))  
         
210     def getSolverOptionsFlux(self):     def getSolverOptionsFlux(self):
211        """        """
212        Returns the solver options used to solve the flux problems        Returns the solver options used to solve the flux problems
         
       *K^{-1} u=F*  
         
213        :return: `SolverOptions`        :return: `SolverOptions`
214        """        """
215        return self.__pde_k.getSolverOptions()        if self.__pde_v == None:
216              return None
217          else:
218              return self.__pde_v.getSolverOptions()
219                
220     def setSolverOptionsFlux(self, options=None):     def setSolverOptionsFlux(self, options=None):
221        """        """
222        Sets the solver options used to solve the flux problems        Sets the solver options used to solve the flux problems
         
       *K^{-1}u=F*  
         
223        If ``options`` is not present, the options are reset to default        If ``options`` is not present, the options are reset to default
         
224        :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.  
225        """        """
226        return self.__pde_v.setSolverOptions(options)        if not self.__pde_v == None:
227              self.__pde_v.setSolverOptions(options)
228            
229     def getSolverOptionsPressure(self):     def getSolverOptionsPressure(self):
230        """        """
231        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*  
         
232        :return: `SolverOptions`        :return: `SolverOptions`
233        """        """
234        return self.__pde_p.getSolverOptions()        return self.__pde_p.getSolverOptions()
# Line 371  class DarcyFlow(object): Line 236  class DarcyFlow(object):
236     def setSolverOptionsPressure(self, options=None):     def setSolverOptionsPressure(self, options=None):
237        """        """
238        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*  
         
239        If ``options`` is not present, the options are reset to default        If ``options`` is not present, the options are reset to default
240                
241        :param options: `SolverOptions`        :param options: `SolverOptions`
242        :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.
243        """        """
244        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  
   
245                
246     def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):     def solve(self, u0, p0):
247        """        """
248        solves the problem.        solves the problem.
249                
       The iteration is terminated if the residual norm is less then self.getTolerance().  
   
250        :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.
251        :type u0: vector value on the domain (e.g. `escript.Data`).        :type u0: vector value on the domain (e.g. `escript.Data`).
252        :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.
253        :type p0: scalar value on the domain (e.g. `escript.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``  
254        :return: flux and pressure        :return: flux and pressure
255        :rtype: ``tuple`` of `escript.Data`.        :rtype: ``tuple`` of `escript.Data`.
256    
       :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}*.  
257        """        """
258        self.verbose=verbose        p0=util.interpolate(p0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
259        if self.useVPIteration:        if self.ref_point_id == None:
260            return self.__solveVP(u0,p0,max_iter,max_num_corrections)            p_ref=0
261        else:        else:
262            X=self.__pde_k.createCoefficient("X")            p_ref=p0.getTupleForGlobalDataPoint(*self.ref_point_id)[0]
263            Y=self.__pde_k.createCoefficient("Y")        p0_hydrostatic=p_ref+util.inner(self.__permeability_invXg_ref, self.__pde_p.getFunctionSpaceForCoefficient("q").getX() - self.ref_point)
264            Y[:self.domain.getDim()]=self.scale*util.tensor_mult(self.__permeability_inv,self.__g)        g_2=self.__g - util.tensor_mult(self.__permeability, self.__permeability_invXg_ref * self.perm_scale)
265            rtmp=self.__f * self.__l2 * self.scale        self.__pde_p.setValue(X=g_2 * 1./self.perm_scale,
266            for i in range(self.domain.getDim()): X[i,i]=rtmp                              Y=self.__f * 1./self.perm_scale,
267            X[self.domain.getDim(),:]=self.__g*self.scale                              y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux * 1./self.perm_scale ),
268            r=self.__pde_k.createCoefficient("r")                              r=p0 - p0_hydrostatic)
269            r[:self.domain.getDim()]=u0*self.scale        pp=self.__pde_p.getSolution()
270            r[self.domain.getDim()]=p0        u = self._getFlux(pp, u0)
271            self.__pde_k.setValue(X=X, Y=Y, r=r)        return u,pp + p0_hydrostatic
           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):  
       """  
       solves the problem.  
272                
273        The iteration is terminated if the residual norm is less than self.getTolerance().     def getFlux(self,p, u0=None):
   
       :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`).  
       :return: flux and pressure  
       :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}*.  
       """  
       rtol=self.getTolerance()  
       atol=self.getAbsoluteTolerance()  
       self.setSubProblemTolerance()  
       num_corrections=0  
       converged=False  
       norm_r=None  
         
       # Eliminate the hydrostatic pressure:  
       if self.verbose: print "DarcyFlux: calculate hydrostatic pressure component."  
       self.__pde_p.setValue(X=self.__g, r=p0, y=-util.inner(self.domain.getNormal(),u0))          
       p0=self.__pde_p.getSolution()  
       g2=self.__g - util.tensor_mult(self.__permeability, util.grad(p0))  
       norm_g2=util.integrate(util.inner(g2,util.tensor_mult(self.__permeability_inv,g2)))**0.5      
   
       p=p0*0  
       if self.solveForFlux:  
      v=u0.copy()  
       else:  
      v=self.__getFlux(p, u0, f=self.__f, g=g2)  
   
       while not converged and norm_g2 > 0:  
      Gp=util.grad(p)  
      KGp=util.tensor_mult(self.__permeability,Gp)  
      if self.verbose:  
         def_p=g2-(v+KGp)  
         def_v=self.__f-util.div(v)  
         print "DarcyFlux: L2: g-v-K*grad(p) = %e (v = %e)."%(self.__L2(def_p),self.__L2(v))  
         print "DarcyFlux: L2: f-div(v) = %e (grad(v) = %e)."%(self.__L2(def_v),self.__L2(util.grad(v)))  
         print "DarcyFlux: K^{-1}-norm of v = %e."%util.integrate(util.inner(v,util.tensor_mult(self.__permeability_inv,v)))**0.5  
         print "DarcyFlux: K^{-1}-norm of g2 = %e."%norm_g2  
         print "DarcyFlux: K-norm of grad(dp) = %e."%util.integrate(util.inner(Gp,KGp))**0.5  
      ATOL=atol+rtol*norm_g2  
      if self.verbose: print "DarcyFlux: absolute tolerance ATOL = %e."%(ATOL,)  
      if norm_r == None or norm_r>ATOL:  
         if num_corrections>max_num_corrections:  
            raise ValueError,"maximum number of correction steps reached."  
         
         if self.solveForFlux:  
            # initial residual is r=K^{-1}*(g-v-K*Gp)+D^*L^{-1}(f-Du)  
            v,r, norm_r=PCG(ArithmeticTuple(util.tensor_mult(self.__permeability_inv,g2-v)-Gp,self.__applWeight(v,self.__f),p),  
                self.__Aprod_v,  
                v,  
                self.__Msolve_PCG_v,  
                self.__inner_PCG_v,  
                atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)  
            p=r[2]  
         else:  
            # initial residual is r=G^*(g2-KGp - v)  
            p,r, norm_r=PCG(ArithmeticTuple(g2-KGp,v),  
                  self.__Aprod_p,  
                  p,  
                  self.__Msolve_PCG_p,  
                  self.__inner_PCG_p,  
                  atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)  
            v=r[1]  
         if self.verbose: print "DarcyFlux: residual norm = %e."%norm_r  
         num_corrections+=1  
      else:  
         if self.verbose: print "DarcyFlux: stopping criterium reached."  
         converged=True  
       return v,p+p0  
   
    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)  
       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=escript.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=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])  
         
    def __Aprod_p(self,dp):  
       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):  
274          """          """
275          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``  
276          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
277          Note that ``g`` and ``f`` are used, see `setValue`.          Notice that ``g`` is used, see `setValue`.
278    
279          :param p: pressure.          :param p: pressure.
280          :type p: scalar value on the domain (e.g. `escript.Data`).          :type p: scalar value on the domain (e.g. `escript.Data`).
281          :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``.
282          :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``
283          :return: flux          :return: flux
284          :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}*  
285          """          """
286      self.setSubProblemTolerance()          p=util.interpolate(p, self.__pde_p.getFunctionSpaceForCoefficient("q"))
287          g=self.__g          if self.ref_point_id == None:
288          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)  
289          else:          else:
290             self.__pde_v.setValue(Y=g-self.__Q(p))              p_ref=p.getTupleForGlobalDataPoint(*self.ref_point_id)[0]
291          return self.__pde_v.getSolution()          p_hydrostatic=p_ref+util.inner(self.__permeability_invXg_ref, self.__pde_p.getFunctionSpaceForCoefficient("q").getX() - self.ref_point)
292            return self._getFlux(p-p_hydrostatic, u0)
293    
294       def _getFlux(self, pp, u0=None):
295            """
296            returns the flux for a given pressure ``pp`` where the flux is equal to
297            ``u0`` on locations where ``location_of_fixed_flux`` is positive (see
298            `setValue`). Notice that ``g`` is used, see `setValue`.
299    
300            :param pp: pressure.
301            :type pp: scalar value on the domain (i.e. `escript.Data`).
302            :param u0: flux on the locations of the domain marked in ``location_of_fixed_flux``.
303            :type u0: vector values on the domain (i.e. `escript.Data`) or ``None``
304            :return: flux
305            :rtype: `escript.Data`
306            """
307            if self.solver  == self.EVAL:
308               u = self.__g - util.tensor_mult(self.__permeability, self.perm_scale * (util.grad(pp) + self.__permeability_invXg_ref))
309            elif self.solver  == self.POST or self.solver  == self.SMOOTH:
310                self.__pde_v.setValue(Y= self.__permeability_invXg - (util.grad(pp) + self.__permeability_invXg_ref))
311                print
312                if u0 == None:
313                   self.__pde_v.setValue(r=escore.Data())
314                else:
315                   if not isinstance(u0, escore.Data) : u0 = escore.Vector(u0, escore.Solution(self.domain))
316                   self.__pde_v.setValue(r=1./self.perm_scale * u0)
317                u= self.__pde_v.getSolution() * self.perm_scale
318            return u
319          
320  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
321       """       """
322       solves       solves
# Line 1074  class StokesProblemCartesian(Homogeneous Line 335  class StokesProblemCartesian(Homogeneous
335              sp.setTolerance()              sp.setTolerance()
336              sp.initialize(...)              sp.initialize(...)
337              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
338                sp.setStokesEquation(...) # new values for some parameters
339                v1,p1=sp.solve(v,p)
340       """       """
341       def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
342           """           """
# Line 1088  class StokesProblemCartesian(Homogeneous Line 351  class StokesProblemCartesian(Homogeneous
351           """           """
352           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
353           self.domain=domain           self.domain=domain
354           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())
355           self.__pde_u.setSymmetryOn()           self.__pde_v.setSymmetryOn()
356            
357           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
358           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
# Line 1097  class StokesProblemCartesian(Homogeneous Line 360  class StokesProblemCartesian(Homogeneous
360    
361           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
362           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
363       self.__pde_proj.setValue(D=1)           self.__pde_proj.setValue(D=1)
364           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
365    
366       def getSolverOptionsVelocity(self):       def getSolverOptionsVelocity(self):
# Line 1106  class StokesProblemCartesian(Homogeneous Line 369  class StokesProblemCartesian(Homogeneous
369            
370       :rtype: `SolverOptions`       :rtype: `SolverOptions`
371       """       """
372       return self.__pde_u.getSolverOptions()           return self.__pde_v.getSolverOptions()
373       def setSolverOptionsVelocity(self, options=None):       def setSolverOptionsVelocity(self, options=None):
374           """           """
375       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 377  class StokesProblemCartesian(Homogeneous
377       :param options: new solver  options       :param options: new solver  options
378       :type options: `SolverOptions`       :type options: `SolverOptions`
379       """       """
380           self.__pde_u.setSolverOptions(options)           self.__pde_v.setSolverOptions(options)
381       def getSolverOptionsPressure(self):       def getSolverOptionsPressure(self):
382           """           """
383       returns the solver options used  solve the equation for pressure.       returns the solver options used  solve the equation for pressure.
384       :rtype: `SolverOptions`       :rtype: `SolverOptions`
385       """       """
386       return self.__pde_prec.getSolverOptions()           return self.__pde_prec.getSolverOptions()
387       def setSolverOptionsPressure(self, options=None):       def setSolverOptionsPressure(self, options=None):
388           """           """
389       set the solver options for solving the equation for pressure.       set the solver options for solving the equation for pressure.
390       :param options: new solver  options       :param options: new solver  options
391       :type options: `SolverOptions`       :type options: `SolverOptions`
392       """       """
393       self.__pde_prec.setSolverOptions(options)           self.__pde_prec.setSolverOptions(options)
394    
395       def setSolverOptionsDiv(self, options=None):       def setSolverOptionsDiv(self, options=None):
396           """           """
# Line 1137  class StokesProblemCartesian(Homogeneous Line 400  class StokesProblemCartesian(Homogeneous
400       :param options: new solver options       :param options: new solver options
401       :type options: `SolverOptions`       :type options: `SolverOptions`
402       """       """
403       self.__pde_proj.setSolverOptions(options)           self.__pde_proj.setSolverOptions(options)
404       def getSolverOptionsDiv(self):       def getSolverOptionsDiv(self):
405           """           """
406       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 408  class StokesProblemCartesian(Homogeneous
408            
409       :rtype: `SolverOptions`       :rtype: `SolverOptions`
410       """       """
411       return self.__pde_proj.getSolverOptions()           return self.__pde_proj.getSolverOptions()
412    
413       def updateStokesEquation(self, v, p):       def updateStokesEquation(self, v, p):
414           """           """
415           updates the Stokes equation to consider dependencies from ``v`` and ``p``           updates the Stokes equation to consider dependencies from ``v`` and ``p``
416           :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.
417           """           """
418           pass           pass
419       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 1171  class StokesProblemCartesian(Homogeneous Line 434  class StokesProblemCartesian(Homogeneous
434          if eta !=None:          if eta !=None:
435              k=util.kronecker(self.domain.getDim())              k=util.kronecker(self.domain.getDim())
436              kk=util.outer(k,k)              kk=util.outer(k,k)
437              self.eta=util.interpolate(eta, escript.Function(self.domain))              self.eta=util.interpolate(eta, escore.Function(self.domain))
438          self.__pde_prec.setValue(D=1/self.eta)              self.__pde_prec.setValue(D=1/self.eta)
439              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)))
440          if restoration_factor!=None:          if restoration_factor!=None:
441              n=self.domain.getNormal()              n=self.domain.getNormal()
442              self.__pde_u.setValue(d=restoration_factor*util.outer(n,n))              self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
443          if fixed_u_mask!=None:          if fixed_u_mask!=None:
444              self.__pde_u.setValue(q=fixed_u_mask)              self.__pde_v.setValue(q=fixed_u_mask)
445          if f!=None: self.__f=f          if f!=None: self.__f=f
446          if surface_stress!=None: self.__surface_stress=surface_stress          if surface_stress!=None: self.__surface_stress=surface_stress
447          if stress!=None: self.__stress=stress          if stress!=None: self.__stress=stress
448    
449       def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):       def initialize(self,f=escore.Data(),fixed_u_mask=escore.Data(),eta=1,surface_stress=escore.Data(),stress=escore.Data(), restoration_factor=0):
450          """          """
451          assigns values to the model parameters          assigns values to the model parameters
452    
# Line 1196  class StokesProblemCartesian(Homogeneous Line 459  class StokesProblemCartesian(Homogeneous
459          :param surface_stress: normal surface stress          :param surface_stress: normal surface stress
460          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
461          :param stress: initial stress          :param stress: initial stress
462      :type stress: `Tensor` object on `FunctionSpace` `Function` or similar          :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
463          """          """
464          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
465    
# Line 1209  class StokesProblemCartesian(Homogeneous Line 472  class StokesProblemCartesian(Homogeneous
472           :rtype: ``float``           :rtype: ``float``
473           """           """
474           self.__pde_proj.setValue(Y=-util.div(v))           self.__pde_proj.setValue(Y=-util.div(v))
475       self.getSolverOptionsDiv().setTolerance(tol)           self.getSolverOptionsDiv().setTolerance(tol)
476       self.getSolverOptionsDiv().setAbsoluteTolerance(0.)           self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
477           out=self.__pde_proj.getSolution()           out=self.__pde_proj.getSolution()
478           return out           return out
479    
# Line 1223  class StokesProblemCartesian(Homogeneous Line 486  class StokesProblemCartesian(Homogeneous
486           :return: inner product of element p and Bv=-div(v)           :return: inner product of element p and Bv=-div(v)
487           :rtype: ``float``           :rtype: ``float``
488           """           """
489           return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))           return util.integrate(util.interpolate(p,escore.Function(self.domain))*util.interpolate(Bv, escore.Function(self.domain)))
490    
491       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
492           """           """
# Line 1234  class StokesProblemCartesian(Homogeneous Line 497  class StokesProblemCartesian(Homogeneous
497           :return: inner product of p0 and p1           :return: inner product of p0 and p1
498           :rtype: ``float``           :rtype: ``float``
499           """           """
500           s0=util.interpolate(p0, escript.Function(self.domain))           s0=util.interpolate(p0, escore.Function(self.domain))
501           s1=util.interpolate(p1, escript.Function(self.domain))           s1=util.interpolate(p1, escore.Function(self.domain))
502           return util.integrate(s0*s1)           return util.integrate(s0*s1)
503    
504       def norm_v(self,v):       def norm_v(self,v):
# Line 1251  class StokesProblemCartesian(Homogeneous Line 514  class StokesProblemCartesian(Homogeneous
514    
515       def getDV(self, p, v, tol):       def getDV(self, p, v, tol):
516           """           """
517           return the value for v for a given p (overwrite)           return the value for v for a given p
518    
519           :param p: a pressure           :param p: a pressure
520           :param v: a initial guess for the value v to return.           :param v: a initial guess for the value v to return.
521           :return: dv given as *Adv=(f-Av-B^*p)*           :return: dv given as *Adv=(f-Av-B^*p)*
522           """           """
523           self.updateStokesEquation(v,p)           self.updateStokesEquation(v,p)
524           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress)           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
525       self.getSolverOptionsVelocity().setTolerance(tol)           self.getSolverOptionsVelocity().setTolerance(tol)
526       self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)           self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
527           if self.__stress.isEmpty():           if self.__stress.isEmpty():
528              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)))
529           else:           else:
530              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)))
531           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
532           return  out           return  out
533    
534       def norm_Bv(self,Bv):       def norm_Bv(self,Bv):
# Line 1275  class StokesProblemCartesian(Homogeneous Line 538  class StokesProblemCartesian(Homogeneous
538          :rtype: equal to the type of p          :rtype: equal to the type of p
539          :note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
540          """          """
541          return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))          return util.sqrt(util.integrate(util.interpolate(Bv, escore.Function(self.domain))**2))
542    
543       def solve_AinvBt(self,p, tol):       def solve_AinvBt(self,p, tol):
544           """           """
# Line 1285  class StokesProblemCartesian(Homogeneous Line 548  class StokesProblemCartesian(Homogeneous
548           :return: the solution of *Av=B^*p*           :return: the solution of *Av=B^*p*
549           :note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
550           """           """
551           self.__pde_u.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))           self.__pde_v.setValue(Y=escore.Data(), y=escore.Data(), X=-p*util.kronecker(self.domain))
552           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
553           return  out           return  out
554    
555       def solve_prec(self,Bv, tol):       def solve_prec(self,Bv, tol):
556           """           """
557           applies preconditioner for for *BA^{-1}B^** to *Bv*           applies preconditioner for for *BA^{-1}B^** to *Bv*
558           with accuracy `self.getSubProblemTolerance()`           with accuracy ``self.getSubProblemTolerance()``
559    
560           :param Bv: velocity increment           :param Bv: velocity increment
561           :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^ * )*
562           :note: boundary conditions on p are zero.           :note: boundary conditions on p are zero.
563           """           """
564           self.__pde_prec.setValue(Y=Bv)           self.__pde_prec.setValue(Y=Bv)
565       self.getSolverOptionsPressure().setTolerance(tol)           self.getSolverOptionsPressure().setTolerance(tol)
566       self.getSolverOptionsPressure().setAbsoluteTolerance(0.)           self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
567           out=self.__pde_prec.getSolution()           out=self.__pde_prec.getSolution()
568           return out           return out

Legend:
Removed from v.3499  
changed lines
  Added in v.4154

  ViewVC Help
Powered by ViewVC 1.1.26