/[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 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC revision 3515 by gross, Thu May 19 08:20:57 2011 UTC
# Line 1  Line 1 
1    # -*- coding: utf-8 -*-
2  ########################################################  ########################################################
3  #  #
4  # Copyright (c) 2003-2009 by University of Queensland  # Copyright (c) 2003-2010 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
6  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
7  #  #
# Line 10  Line 11 
11  #  #
12  ########################################################  ########################################################
13    
14  __copyright__="""Copyright (c) 2003-2008 by University of Queensland  __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
16  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
17  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
# Line 21  __url__="https://launchpad.net/escript-f Line 22  __url__="https://launchpad.net/escript-f
22  """  """
23  Some models for flow  Some models for flow
24    
25  @var __author__: name of author  :var __author__: name of author
26  @var __copyright__: copyrights  :var __copyright__: copyrights
27  @var __license__: licence agreement  :var __license__: licence agreement
28  @var __url__: url entry point on documentation  :var __url__: url entry point on documentation
29  @var __version__: version  :var __version__: version
30  @var __date__: date of the version  :var __date__: date of the version
31  """  """
32    
33  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
34    
35  from escript import *  import escript
36  import util  import util
37  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
38  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39    
40  class DarcyFlow(object):  class DarcyFlow(object):
41      """     """
42      solves the problem     solves the problem
43      
44      M{u_i+k_{ij}*p_{,j} = g_i}     *u_i+k_{ij}*p_{,j} = g_i*
45      M{u_{i,i} = f}     *u_{i,i} = f*
46      
47      where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,     where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
48      
49      @note: The problem is solved in a least squares formulation.     :cvar SIMPLE: simple solver
50      """     :cvar POST: solver using global postprocessing of flux
51       :cvar STAB: solver uses (non-symmetric) stabilization
52      def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):     :cvar SYMSTAB: solver uses symmetric stabilization
53          """     """
54          initializes the Darcy flux problem     SIMPLE="SIMPLE"
55          @param domain: domain of the problem     POST="POST"
56          @type domain: L{Domain}     STAB="STAB"
57      @param useReduced: uses reduced oreder on flux and pressure     SYMSTAB="SYMSTAB"
58      @type useReduced: C{bool}     def __init__(self, domain, useReduced=False, solver="SYMSTAB", verbose=False, w=1.):
59      @param adaptSubTolerance: switches on automatic subtolerance selection        """
60      @type adaptSubTolerance: C{bool}            initializes the Darcy flux problem
61          """        :param domain: domain of the problem
62          self.domain=domain        :type domain: `Domain`
63          if weight == None:        :param useReduced: uses reduced oreder on flux and pressure
64             s=self.domain.getSize()        :type useReduced: ``bool``
65             self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2        :param solver: solver method
66             # self.__l=(3.*util.longestEdge(self.domain))**2        :type solver: in [`DarcyFlow.SIMPLE`, `DarcyFlow.POST', `DarcyFlow.STAB`, `DarcyFlow.SYMSTAB` ]
67             # self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2        :param verbose: if ``True`` some information on the iteration progress are printed.
68          else:        :type verbose: ``bool``
69             self.__l=weight        :param w: weighting factor for `DarcyFlow.POST` solver
70          self.__pde_v=LinearPDESystem(domain)        :type w: ``float``
71          if useReduced: self.__pde_v.setReducedOrderOn()        
72          self.__pde_v.setSymmetryOn()        """
73          self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))        self.domain=domain
74          self.__pde_p=LinearSinglePDE(domain)        self.solver=solver
75          self.__pde_p.setSymmetryOn()        self.useReduced=useReduced
76          if useReduced: self.__pde_p.setReducedOrderOn()        self.verbose=verbose
77          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))        self.scale=1.
78          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))        
79          self.setTolerance()        
80          self.setAbsoluteTolerance()        self.__pde_v=LinearPDESystem(domain)
81      self.__adaptSubTolerance=adaptSubTolerance        self.__pde_v.setSymmetryOn()
82      self.verbose=False        if self.useReduced: self.__pde_v.setReducedOrderOn()
83      def getSolverOptionsFlux(self):        self.__pde_p=LinearSinglePDE(domain)
84      """        self.__pde_p.setSymmetryOn()
85      Returns the solver options used to solve the flux problems        if self.useReduced: self.__pde_p.setReducedOrderOn()
86              
87      M{(I+D^*D)u=F}        if self.solver  == self.SIMPLE:
88             if self.verbose: print "DarcyFlow: simple solver is used."
89      @return: L{SolverOptions}           self.__pde_v.setValue(D=util.kronecker(self.domain.getDim()))
90      """        elif self.solver  == self.POST:
91      return self.__pde_v.getSolverOptions()       self.w=w
92      def setSolverOptionsFlux(self, options=None):       if util.inf(w)<0.:
93      """          raise ValueError,"Weighting factor must be non-negative."
94      Sets the solver options used to solve the flux problems       if self.verbose: print "DarcyFlow: global postprocessing of flux is used."
95              elif self.solver  == self.STAB:
96      M{(I+D^*D)u=F}        if self.verbose: print "DarcyFlow: (non-symmetric) stabilization is used."
97              elif  self.solver  == self.SYMSTAB:
98      If C{options} is not present, the options are reset to default        if self.verbose: print "DarcyFlow: symmetric stabilization is used."
99      @param options: L{SolverOptions}        else:
100      @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.      raise ValueError,"unknown solver %s"%self.solver
101      """        self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
102      return self.__pde_v.setSolverOptions(options)        self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
103      def getSolverOptionsPressure(self):        self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
104      """        self.location_of_fixed_flux = escript.Vector(0, self.__pde_v.getFunctionSpaceForCoefficient("q"))
105      Returns the solver options used to solve the pressure problems        self.setTolerance()
106            
107      M{(Q^*Q)p=Q^*G}          
108       def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
109          """
110          assigns values to model parameters
111    
112          :param f: volumetic sources/sinks
113          :type f: scalar value on the domain (e.g. `escript.Data`)
114          :param g: flux sources/sinks
115          :type g: vector values on the domain (e.g. `escript.Data`)
116          :param location_of_fixed_pressure: mask for locations where pressure is fixed
117          :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
118          :param location_of_fixed_flux:  mask for locations where flux is fixed.
119          :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
120          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
121          :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
122    
123          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
124          :note: at any point on the boundary of the domain the pressure
125                 (``location_of_fixed_pressure`` >0) or the normal component of the
126                 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
127                 is along the *x_i* axis.
128    
129          """
130          if location_of_fixed_pressure!=None:
131               self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)
132               self.__pde_p.setValue(q=self.location_of_fixed_pressure)
133          if location_of_fixed_flux!=None:
134              self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
135              self.__pde_v.setValue(q=self.location_of_fixed_flux)
136          
137                
138          # pressure  is rescaled by the factor 1/self.scale
139          if permeability!=None:
140            
141      @return: L{SolverOptions}       perm=util.interpolate(permeability,self.__pde_v.getFunctionSpaceForCoefficient("A"))
142      """           V=util.vol(self.domain)
143      return self.__pde_p.getSolverOptions()           l=V**(1./self.domain.getDim())
144      def setSolverOptionsPressure(self, options=None):          
145      """       if perm.getRank()==0:
146      Sets the solver options used to solve the pressure problems          perm_inv=(1./perm)
147                    self.scale=util.integrate(perm_inv)/V*l
148      M{(Q^*Q)p=Q^*G}          perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
149                perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
150      If C{options} is not present, the options are reset to default          
151      @param options: L{SolverOptions}          
152      @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.       elif perm.getRank()==2:
153      """          perm_inv=util.inverse(perm)
154      return self.__pde_p.setSolverOptions(options)              self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)/V)*l
155            perm_inv*=(1./self.scale)
156      def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):          perm=perm*self.scale
157          """       else:
158          assigns values to model parameters          raise ValueError,"illegal rank of permeability."
159            
160          @param f: volumetic sources/sinks       self.__permeability=perm
161          @type f: scalar value on the domain (e.g. L{Data})       self.__permeability_inv=perm_inv
162          @param g: flux sources/sinks       if self.verbose: print "DarcyFlow: scaling factor for pressure is %e."%self.scale
         @type g: vector values on the domain (e.g. L{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. L{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. L{Data})  
         @param permeability: permeability tensor. If scalar C{s} is given the tensor with  
                              C{s} on the main diagonal is used. If vector C{v} is given the tensor with  
                              C{v} on the main diagonal is used.  
         @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})  
   
         @note: the values of parameters which are not set by calling C{setValue} are not altered.  
         @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)  
                or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal  
                is along the M{x_i} axis.  
         """  
         if f !=None:  
            f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))  
            if f.isEmpty():  
                f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))  
            else:  
                if f.getRank()>0: raise ValueError,"illegal rank of f."  
            self.__f=f  
         if g !=None:  
            g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))  
            if g.isEmpty():  
              g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))  
            else:  
              if not g.getShape()==(self.domain.getDim(),):  
                raise ValueError,"illegal shape of g"  
            self.__g=g  
   
         if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)  
         if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)  
   
         if permeability!=None:  
            perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))  
            if perm.getRank()==0:  
                perm=perm*util.kronecker(self.domain.getDim())  
            elif perm.getRank()==1:  
                perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm  
                for i in range(self.domain.getDim()): perm[i,i]=perm2[i]  
            elif perm.getRank()==2:  
               pass  
            else:  
               raise ValueError,"illegal rank of permeability."  
            self.__permeability=perm  
            self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))  
   
     def setTolerance(self,rtol=1e-4):  
         """  
         sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if  
   
         M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }  
   
         where C{atol} is an absolut tolerance (see L{setAbsoluteTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.  
   
         @param rtol: relative tolerance for the pressure  
         @type rtol: non-negative C{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: C{float}  
         """  
         return self.__rtol  
   
     def setAbsoluteTolerance(self,atol=0.):  
         """  
         sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if  
   
         M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }  
   
         where C{rtol} is an absolut tolerance (see L{setTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.  
   
         @param atol: absolute tolerance for the pressure  
         @type atol: non-negative C{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: C{float}  
        """  
        return self.__atol  
     def getSubProblemTolerance(self):  
     """  
     Returns a suitable subtolerance  
     @type: C{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 C{location_of_fixed_flux} the value of C{u0} is kept unchanged.  
          @type u0: vector value on the domain (e.g. L{Data}).  
          @param p0: initial guess for the pressure. At locations in the domain marked by C{location_of_fixed_pressure} the value of C{p0} is kept unchanged.  
          @type p0: scalar value on the domain (e.g. L{Data}).  
          @param verbose: if set some information on iteration progress are printed  
          @type verbose: C{bool}  
          @return: flux and pressure  
          @rtype: C{tuple} of L{Data}.  
   
          @note: The problem is solved as a least squares form  
   
          M{(I+D^*D)u+Qp=D^*f+g}  
          M{Q^*u+Q^*Qp=Q^*g}  
   
          where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.  
          We eliminate the flux form the problem by setting  
   
          M{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  
   
          M{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 M{Q^*Q}). In each iteration step  
          PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.  
          """  
          self.verbose=verbose  
          rtol=self.getTolerance()  
          atol=self.getAbsoluteTolerance()  
      self.setSubProblemTolerance()  
163            
164           num_corrections=0       if self.solver  == self.SIMPLE:
165           converged=False          self.__pde_p.setValue(A=self.__permeability)
166           p=p0       elif self.solver  == self.POST:
167           norm_r=None          self.__pde_p.setValue(A=self.__permeability)
168           while not converged:          k=util.kronecker(self.domain.getDim())
169                 v=self.getFlux(p, fixed_flux=u0)          self.lamb = self.w*util.length(perm_inv)*l
170                 Qp=self.__Q(p)          self.__pde_v.setValue(D=self.__permeability_inv, A=self.lamb*self.domain.getSize()*util.outer(k,k))
171                 norm_v=self.__L2(v)       elif self.solver  == self.STAB:
172                 norm_Qp=self.__L2(Qp)          self.__pde_p.setValue(A=0.5*self.__permeability)
173                 if norm_v == 0.:          self.__pde_v.setValue(D=0.5*self.__permeability_inv)
174                    if norm_Qp == 0.:       elif  self.solver  == self.SYMSTAB:
175                       return v,p          self.__pde_p.setValue(A=0.5*self.__permeability)
176                    else:          self.__pde_v.setValue(D=0.5*self.__permeability_inv)
177                      fac=norm_Qp  
178                 else:        if g != None:
179                    if norm_Qp == 0.:      g=util.interpolate(g, self.__pde_v.getFunctionSpaceForCoefficient("Y"))
180                      fac=norm_v      if g.isEmpty():
181                    else:            g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
182                      fac=2./(1./norm_v+1./norm_Qp)      else:
183                 ATOL=(atol+rtol*fac)          if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g"
184                 if self.verbose:      self.__g=g
185                      print "DarcyFlux: L2 norm of v = %e."%norm_v        if f !=None:
186                      print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp       f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
187                      print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-Qp)**2)**(0.5),)       if f.isEmpty():      
188                      print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)            f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
189                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL       else:
190                 if norm_r == None or norm_r>ATOL:           if f.getRank()>0: raise ValueError,"illegal rank of f."
191                     if num_corrections>max_num_corrections:       self.__f=f
192                           raise ValueError,"maximum number of correction steps reached."     def getSolverOptionsFlux(self):
193                     p,r, norm_r=PCG(self.__g-util.interpolate(v,Function(self.domain))-Qp,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.5*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)        """
194                     num_corrections+=1        Returns the solver options used to solve the flux problems
195                 else:        :return: `SolverOptions`
196                     converged=True        """
197           return v,p        return self.__pde_v.getSolverOptions()
198      def __L2(self,v):        
199           return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))     def setSolverOptionsFlux(self, options=None):
200          """
201      def __Q(self,p):        Sets the solver options used to solve the flux problems
202            return util.tensor_mult(self.__permeability,util.grad(p))        If ``options`` is not present, the options are reset to default
203          :param options: `SolverOptions`
204      def __Aprod(self,dp):        """
205            if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"        return self.__pde_v.setSolverOptions(options)
206            Qdp=self.__Q(dp)      
207            self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())     def getSolverOptionsPressure(self):
208            du=self.__pde_v.getSolution()        """
209            # self.__pde_v.getOperator().saveMM("proj.mm")        Returns the solver options used to solve the pressure problems
210            return Qdp+du        :return: `SolverOptions`
211      def __inner_GMRES(self,r,s):        """
212           return util.integrate(util.inner(r,s))        return self.__pde_p.getSolverOptions()
213          
214      def __inner_PCG(self,p,r):     def setSolverOptionsPressure(self, options=None):
215           return util.integrate(util.inner(self.__Q(p), r))        """
216          Sets the solver options used to solve the pressure problems
217      def __Msolve_PCG(self,r):        If ``options`` is not present, the options are reset to default
218        if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"        
219            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())        :param options: `SolverOptions`
220            # self.__pde_p.getOperator().saveMM("prec.mm")        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
221            return self.__pde_p.getSolution()        """
222          return self.__pde_p.setSolverOptions(options)
223      def getFlux(self,p=None, fixed_flux=Data()):        
224          """     def setTolerance(self,rtol=1e-4):
225          returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}        """
226          on locations where C{location_of_fixed_flux} is positive (see L{setValue}).        sets the relative tolerance ``rtol`` for the pressure for the stabelized solvers.
227          Note that C{g} and C{f} are used, see L{setValue}.        
228          :param rtol: relative tolerance for the pressure
229          @param p: pressure.        :type rtol: non-negative ``float``
230          @type p: scalar value on the domain (e.g. L{Data}).        """
231          @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.        if rtol<0:
232          @type fixed_flux: vector values on the domain (e.g. L{Data}).       raise ValueError,"Relative tolerance needs to be non-negative."
233          @param tol: relative tolerance to be used.        self.__rtol=rtol
234          @type tol: positive C{float}.        
235          @return: flux     def getTolerance(self):
236          @rtype: L{Data}        """
237          @note: the method uses the least squares solution M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}}        returns the relative tolerance
238                 for the permeability M{k_{ij}}        :return: current relative tolerance
239          """        :rtype: ``float``
240      self.setSubProblemTolerance()        """
241          g=self.__g        return self.__rtol
242          f=self.__f        
243          self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)     def solve(self,u0,p0, max_iter=100, iter_restart=20):
244          if p == None:        """
245             self.__pde_v.setValue(Y=g)        solves the problem.
246          else:        
247             self.__pde_v.setValue(Y=g-self.__Q(p))        The iteration is terminated if the residual norm is less then self.getTolerance().
248          return self.__pde_v.getSolution()  
249          :param u0: initial guess for the flux. At locations in the domain marked by ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.
250          :type u0: vector value on the domain (e.g. `escript.Data`).
251          :param p0: initial guess for the pressure. At locations in the domain marked by ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.
252          :type p0: scalar value on the domain (e.g. `escript.Data`).
253          :param max_iter: maximum number of (outer) iteration steps for the stabilization solvers,
254          :type max_iter: ``int``
255          :param iter_restart: number of steps after which the iteration is restarted. The larger ``iter_restart`` the larger the required memory.
256                               A small value for ``iter_restart`` may require a large number of iteration steps or may even lead to a failure
257                               of the iteration. ``iter_restart`` is relevant for the stabilization solvers only.
258          :type iter_restart: ``int``
259          :return: flux and pressure
260          :rtype: ``tuple`` of `escript.Data`.
261    
262          """
263          # rescale initial guess:
264          p0=p0/self.scale
265          if self.solver  == self.SIMPLE or self.solver  == self.POST :
266            self.__pde_p.setValue(X=self.__g ,
267                                  Y=self.__f,
268                                  y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux),
269                                  r=p0)
270            p=self.__pde_p.getSolution()
271            u = self.getFlux(p, u0)
272          elif  self.solver  == self.STAB:
273        u,p = self.__solve_STAB(u0,p0, max_iter, iter_restart)
274          elif  self.solver  == self.SYMSTAB:
275        u,p = self.__solve_SYMSTAB(u0,p0, max_iter, iter_restart)
276        
277          if self.verbose:
278            KGp=util.tensor_mult(self.__permeability,util.grad(p))
279            def_p=self.__g-(u+KGp)
280            def_v=self.__f-util.div(u, self.__pde_v.getFunctionSpaceForCoefficient("X"))
281            print "DarcyFlux: |g-u-K*grad(p)|_2 = %e (|u|_2 = %e)."%(self.__L2(def_p),self.__L2(u))
282            print "DarcyFlux: |f-div(u)|_2 = %e (|grad(u)|_2 = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))
283          #rescale result
284          p=p*self.scale
285          return u,p
286          
287       def getFlux(self,p, u0=None):
288            """
289            returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
290            on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
291            Notice that ``g`` and ``f`` are used, see `setValue`.
292    
293            :param p: pressure.
294            :type p: scalar value on the domain (e.g. `escript.Data`).
295            :param u0: flux on the locations of the domain marked be ``location_of_fixed_flux``.
296            :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
297            :return: flux
298            :rtype: `escript.Data`
299            """
300            if self.solver  == self.SIMPLE or self.solver  == self.POST  :
301                KGp=util.tensor_mult(self.__permeability,util.grad(p))
302                self.__pde_v.setValue(Y=self.__g-KGp, X=escript.Data())
303                if u0 == None:
304               self.__pde_v.setValue(r=escript.Data())
305            else:
306               self.__pde_v.setValue(r=u0)
307                u= self.__pde_v.getSolution()
308        elif self.solver  == self.POST:
309                self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g)-util.grad(p),
310                                      X=self.lamb * self.__f * util.kronecker(self.domain.getDim()))
311                if u0 == None:
312               self.__pde_v.setValue(r=escript.Data())
313            else:
314               self.__pde_v.setValue(r=u0)
315                u= self.__pde_v.getSolution()
316        elif self.solver  == self.STAB:
317             gp=util.grad(p)
318             self.__pde_v.setValue(Y=0.5*(util.tensor_mult(self.__permeability_inv,self.__g)+gp),
319                                   X= p * util.kronecker(self.domain.getDim()),
320                                   y= - p * self.domain.getNormal())                          
321             if u0 == None:
322               self.__pde_v.setValue(r=escript.Data())
323             else:
324               self.__pde_v.setValue(r=u0)
325             u= self.__pde_v.getSolution()
326        elif  self.solver  == self.SYMSTAB:
327             gp=util.grad(p)
328             self.__pde_v.setValue(Y=0.5*(util.tensor_mult(self.__permeability_inv,self.__g)-gp),
329                                   X= escript.Data() ,
330                                   y= escript.Data() )                          
331             if u0 == None:
332               self.__pde_v.setValue(r=escript.Data())
333             else:
334               self.__pde_v.setValue(r=u0)
335             u= self.__pde_v.getSolution()
336        return u
337          
338        
339       def __solve_STAB(self, u0, p0, max_iter, iter_restart):
340              # p0 is used as an initial guess
341          u=self.getFlux(p0, u0)  
342              self.__pde_p.setValue( Y=self.__f-util.div(u),
343                                     X=0.5*(self.__g - u - util.tensor_mult(self.__permeability,util.grad(p0)) ),
344                                     y= escript.Data(),
345                                     r=escript.Data())
346    
347          dp=self.__pde_p.getSolution()
348          p=GMRES(dp,
349                  self.__STAB_Aprod,
350              p0,
351              self.__inner,
352              atol=self.__norm(p0+dp)*self.getTolerance() ,
353              rtol=0.,
354              iter_max=max_iter,
355              iter_restart=iter_restart,
356              verbose=self.verbose,P_R=None)
357                
358              u=self.getFlux(p, u0)
359              return u,p
360    
361       def __solve_SYMSTAB(self, u0, p0, max_iter, iter_restart):
362              # p0 is used as an initial guess
363          u=self.getFlux(p0, u0)  
364              self.__pde_p.setValue( Y= self.__f,
365                                     X=  0.5*(self.__g + u - util.tensor_mult(self.__permeability,util.grad(p0)) ),
366                                     y=  -  util.inner(self.domain.getNormal(), u),
367                                     r=escript.Data())
368          dp=self.__pde_p.getSolution()
369          
370          p=GMRES(dp,
371                  self.__SYMSTAB_Aprod,
372              p0,
373              self.__inner,
374              atol=self.__norm(p0+dp)*self.getTolerance() ,
375              rtol=0.,
376              iter_max=max_iter,
377              iter_restart=iter_restart,
378              verbose=self.verbose,P_R=None)
379                
380              u=self.getFlux(p, u0)
381              return u,p
382    
383       def __L2(self,v):
384             return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))      
385      
386       def __norm(self,r):
387             return util.sqrt(self.__inner(r,r))
388            
389       def __inner(self,r,s):
390             return util.integrate(util.inner(r,s), escript.Function(self.domain))
391            
392       def __STAB_Aprod(self,p):
393          gp=util.grad(p)
394          self.__pde_v.setValue(Y=-0.5*gp,
395                                X=-p*util.kronecker(self.__pde_v.getDomain()),
396                                y= p * self.domain.getNormal(),  
397                                r=escript.Data())
398          u = -self.__pde_v.getSolution()
399          self.__pde_p.setValue(Y=util.div(u),
400                                X=0.5*(u+util.tensor_mult(self.__permeability,gp)),
401                                y=escript.Data(),
402                                r=escript.Data())
403        
404          return  self.__pde_p.getSolution()
405      
406       def __SYMSTAB_Aprod(self,p):
407          gp=util.grad(p)
408          self.__pde_v.setValue(Y=0.5*gp ,
409                                X=escript.Data(),
410                                y=escript.Data(),  
411                                r=escript.Data())
412          u = -self.__pde_v.getSolution()
413          self.__pde_p.setValue(Y=escript.Data(),
414                                X=0.5*(-u+util.tensor_mult(self.__permeability,gp)),
415                                y=escript.Data(),
416                                r=escript.Data())
417        
418          return  self.__pde_p.getSolution()
419          
420    
421  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
422       """       """
# Line 380  class StokesProblemCartesian(Homogeneous Line 437  class StokesProblemCartesian(Homogeneous
437              sp.initialize(...)              sp.initialize(...)
438              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
439       """       """
440       def __init__(self,domain,adaptSubTolerance=True, **kwargs):       def __init__(self,domain,**kwargs):
441           """           """
442           initialize the Stokes Problem           initialize the Stokes Problem
443    
444           @param domain: domain of the problem. The approximation order needs to be two.           The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
445           @type domain: L{Domain}           LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
446       @param adaptSubTolerance: If True the tolerance for subproblem is set automatically.           with macro elements for the pressure.
447       @type adaptSubTolerance: C{bool}  
448           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.           :param domain: domain of the problem.
449             :type domain: `Domain`
450           """           """
451           HomogeneousSaddlePointProblem.__init__(self,adaptSubTolerance=adaptSubTolerance,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
452           self.domain=domain           self.domain=domain
453           self.vol=util.integrate(1.,Function(self.domain))           self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
454           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_v.setSymmetryOn()
          self.__pde_u.setSymmetryOn()  
455            
456           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
457           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
# Line 409  class StokesProblemCartesian(Homogeneous Line 466  class StokesProblemCartesian(Homogeneous
466           """           """
467       returns the solver options used  solve the equation for velocity.       returns the solver options used  solve the equation for velocity.
468            
469       @rtype: L{SolverOptions}       :rtype: `SolverOptions`
470       """       """
471       return self.__pde_u.getSolverOptions()       return self.__pde_v.getSolverOptions()
472       def setSolverOptionsVelocity(self, options=None):       def setSolverOptionsVelocity(self, options=None):
473           """           """
474       set the solver options for solving the equation for velocity.       set the solver options for solving the equation for velocity.
475            
476       @param options: new solver  options       :param options: new solver  options
477       @type options: L{SolverOptions}       :type options: `SolverOptions`
478       """       """
479           self.__pde_u.setSolverOptions(options)           self.__pde_v.setSolverOptions(options)
480       def getSolverOptionsPressure(self):       def getSolverOptionsPressure(self):
481           """           """
482       returns the solver options used  solve the equation for pressure.       returns the solver options used  solve the equation for pressure.
483       @rtype: L{SolverOptions}       :rtype: `SolverOptions`
484       """       """
485       return self.__pde_prec.getSolverOptions()       return self.__pde_prec.getSolverOptions()
486       def setSolverOptionsPressure(self, options=None):       def setSolverOptionsPressure(self, options=None):
487           """           """
488       set the solver options for solving the equation for pressure.       set the solver options for solving the equation for pressure.
489       @param options: new solver  options       :param options: new solver  options
490       @type options: L{SolverOptions}       :type options: `SolverOptions`
491       """       """
492       self.__pde_prec.setSolverOptions(options)       self.__pde_prec.setSolverOptions(options)
493    
# Line 439  class StokesProblemCartesian(Homogeneous Line 496  class StokesProblemCartesian(Homogeneous
496       set the solver options for solving the equation to project the divergence of       set the solver options for solving the equation to project the divergence of
497       the velocity onto the function space of presure.       the velocity onto the function space of presure.
498            
499       @param options: new solver options       :param options: new solver options
500       @type options: L{SolverOptions}       :type options: `SolverOptions`
501       """       """
502       self.__pde_prec.setSolverOptions(options)       self.__pde_proj.setSolverOptions(options)
503       def getSolverOptionsDiv(self):       def getSolverOptionsDiv(self):
504           """           """
505       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
506       the velocity onto the function space of presure.       the velocity onto the function space of presure.
507            
508       @rtype: L{SolverOptions}       :rtype: `SolverOptions`
509       """       """
510       return self.__pde_prec.getSolverOptions()       return self.__pde_proj.getSolverOptions()
511       def setSubProblemTolerance(self):  
512         def updateStokesEquation(self, v, p):
513           """           """
514       Updates the tolerance for subproblems           updates the Stokes equation to consider dependencies from ``v`` and ``p``
515             :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
516           """           """
517       if self.adaptSubTolerance():           pass
518               sub_tol=self.getSubProblemTolerance()       def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
519           self.getSolverOptionsDiv().setTolerance(sub_tol)          """
520           self.getSolverOptionsDiv().setAbsoluteTolerance(0.)          assigns new values to the model parameters.
521           self.getSolverOptionsPressure().setTolerance(sub_tol)  
522           self.getSolverOptionsPressure().setAbsoluteTolerance(0.)          :param f: external force
523           self.getSolverOptionsVelocity().setTolerance(sub_tol)          :type f: `Vector` object in `FunctionSpace` `Function` or similar
524           self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)          :param fixed_u_mask: mask of locations with fixed velocity.
525                    :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
526            :param eta: viscosity
527            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
528            :param surface_stress: normal surface stress
529            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
530            :param stress: initial stress
531        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
532            """
533            if eta !=None:
534                k=util.kronecker(self.domain.getDim())
535                kk=util.outer(k,k)
536                self.eta=util.interpolate(eta, escript.Function(self.domain))
537            self.__pde_prec.setValue(D=1/self.eta)
538                self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
539            if restoration_factor!=None:
540                n=self.domain.getNormal()
541                self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
542            if fixed_u_mask!=None:
543                self.__pde_v.setValue(q=fixed_u_mask)
544            if f!=None: self.__f=f
545            if surface_stress!=None: self.__surface_stress=surface_stress
546            if stress!=None: self.__stress=stress
547    
548       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):       def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
549          """          """
550          assigns values to the model parameters          assigns values to the model parameters
551    
552          @param f: external force          :param f: external force
553          @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar          :type f: `Vector` object in `FunctionSpace` `Function` or similar
554          @param fixed_u_mask: mask of locations with fixed velocity.          :param fixed_u_mask: mask of locations with fixed velocity.
555          @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar          :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
556          @param eta: viscosity          :param eta: viscosity
557          @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar          :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
558          @param surface_stress: normal surface stress          :param surface_stress: normal surface stress
559          @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
560          @param stress: initial stress          :param stress: initial stress
561      @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar      :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
         @note: All values needs to be set.  
562          """          """
563          self.eta=eta          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
         A =self.__pde_u.createCoefficient("A")  
     self.__pde_u.setValue(A=Data())  
         for i in range(self.domain.getDim()):  
         for j in range(self.domain.getDim()):  
             A[i,j,j,i] += 1.  
             A[i,j,i,j] += 1.  
     self.__pde_prec.setValue(D=1/self.eta)  
         self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask)  
         self.__f=f  
         self.__surface_stress=surface_stress  
         self.__stress=stress  
564    
565       def Bv(self,v):       def Bv(self,v,tol):
566           """           """
567           returns inner product of element p and div(v)           returns inner product of element p and div(v)
568    
569           @param p: a pressure increment           :param v: a residual
570           @param v: a residual           :return: inner product of element p and div(v)
571           @return: inner product of element p and div(v)           :rtype: ``float``
572           @rtype: C{float}           """
573           """           self.__pde_proj.setValue(Y=-util.div(v))
574           self.__pde_proj.setValue(Y=-util.div(v))       self.getSolverOptionsDiv().setTolerance(tol)
575           return self.__pde_proj.getSolution()       self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
576             out=self.__pde_proj.getSolution()
577             return out
578    
579       def inner_pBv(self,p,Bv):       def inner_pBv(self,p,Bv):
580           """           """
581           returns inner product of element p and Bv=-div(v)           returns inner product of element p and Bv=-div(v)
582    
583           @param p: a pressure increment           :param p: a pressure increment
584           @param v: a residual           :param Bv: a residual
585           @return: inner product of element p and Bv=-div(v)           :return: inner product of element p and Bv=-div(v)
586           @rtype: C{float}           :rtype: ``float``
587           """           """
588           return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))           return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
589    
590       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
591           """           """
592           Returns inner product of p0 and p1           Returns inner product of p0 and p1
593    
594           @param p0: a pressure           :param p0: a pressure
595           @param p1: a pressure           :param p1: a pressure
596           @return: inner product of p0 and p1           :return: inner product of p0 and p1
597           @rtype: C{float}           :rtype: ``float``
598           """           """
599           s0=util.interpolate(p0/self.eta,Function(self.domain))           s0=util.interpolate(p0, escript.Function(self.domain))
600           s1=util.interpolate(p1/self.eta,Function(self.domain))           s1=util.interpolate(p1, escript.Function(self.domain))
601           return util.integrate(s0*s1)           return util.integrate(s0*s1)
602    
603       def norm_v(self,v):       def norm_v(self,v):
604           """           """
605           returns the norm of v           returns the norm of v
606    
607           @param v: a velovity           :param v: a velovity
608           @return: norm of v           :return: norm of v
609           @rtype: non-negative C{float}           :rtype: non-negative ``float``
610           """           """
611           return util.sqrt(util.integrate(util.length(util.grad(v))))           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
612    
613       def getV(self, p, v0):  
614         def getDV(self, p, v, tol):
615           """           """
616           return the value for v for a given p (overwrite)           return the value for v for a given p (overwrite)
617    
618           @param p: a pressure           :param p: a pressure
619           @param v0: a initial guess for the value v to return.           :param v: a initial guess for the value v to return.
620           @return: v given as M{v= A^{-1} (f-B^*p)}           :return: dv given as *Adv=(f-Av-B^*p)*
621           """           """
622           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)           self.updateStokesEquation(v,p)
623             self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
624         self.getSolverOptionsVelocity().setTolerance(tol)
625         self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
626           if self.__stress.isEmpty():           if self.__stress.isEmpty():
627              self.__pde_u.setValue(X=p*util.kronecker(self.domain))              self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
628           else:           else:
629              self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))              self.__pde_v.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
630           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
631           return  out           return  out
632    
633       def norm_Bv(self,Bv):       def norm_Bv(self,Bv):
634          """          """
635          Returns Bv (overwrite).          Returns Bv (overwrite).
636    
637          @rtype: equal to the type of p          :rtype: equal to the type of p
638          @note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
639          """          """
640          return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))          return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
641    
642       def solve_AinvBt(self,p):       def solve_AinvBt(self,p, tol):
643           """           """
644           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}           Solves *Av=B^*p* with accuracy `tol`
645    
646           @param p: a pressure increment           :param p: a pressure increment
647           @return: the solution of M{Av=B^*p}           :return: the solution of *Av=B^*p*
648           @note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
649           """           """
650           self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))           self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
651           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
652           return  out           return  out
653    
654       def solve_prec(self,Bv):       def solve_prec(self,Bv, tol):
655           """           """
656           applies preconditioner for for M{BA^{-1}B^*} to M{Bv}           applies preconditioner for for *BA^{-1}B^** to *Bv*
657           with accuracy L{self.getSubProblemTolerance()}           with accuracy `self.getSubProblemTolerance()`
658    
659           @param v: velocity increment           :param Bv: velocity increment
660           @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}           :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
661           @note: boundary conditions on p are zero.           :note: boundary conditions on p are zero.
662           """           """
663           self.__pde_prec.setValue(Y=Bv)           self.__pde_prec.setValue(Y=Bv)
664           return self.__pde_prec.getSolution()       self.getSolverOptionsPressure().setTolerance(tol)
665         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
666             out=self.__pde_prec.getSolution()
667             return out

Legend:
Removed from v.2548  
changed lines
  Added in v.3515

  ViewVC Help
Powered by ViewVC 1.1.26