/[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 3630 by gross, Wed Oct 19 03:04:53 2011 UTC revision 3990 by caltinay, Tue Sep 25 05:03:20 2012 UTC
# Line 1  Line 1 
1  # -*- coding: utf-8 -*-  # -*- coding: utf-8 -*-
2  ########################################################  ##############################################################################
3  #  #
4  # Copyright (c) 2003-2010 by University of Queensland  # Copyright (c) 2003-2012 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # http://www.uq.edu.au
 # http://www.uq.edu.au/esscc  
6  #  #
7  # Primary Business: Queensland, Australia  # Primary Business: Queensland, Australia
8  # Licensed under the Open Software License version 3.0  # Licensed under the Open Software License version 3.0
9  # http://www.opensource.org/licenses/osl-3.0.php  # http://www.opensource.org/licenses/osl-3.0.php
10  #  #
11  ########################################################  # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    # Development since 2012 by School of Earth Sciences
13    #
14    ##############################################################################
15    
16  __copyright__="""Copyright (c) 2003-2010 by University of Queensland  __copyright__="""Copyright (c) 2003-2012 by University of Queensland
17  Earth Systems Science Computational Center (ESSCC)  http://www.uq.edu.au
 http://www.uq.edu.au/esscc  
18  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
19  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
20  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
# Line 32  Some models for flow Line 33  Some models for flow
33    
34  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
35    
36  import escript  from . import escript
37  import util  from . import util
38  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions  from .linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
39  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES  from .pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
40    
41  class DarcyFlow(object):  class DarcyFlow(object):
42     """     """
# Line 57  class DarcyFlow(object): Line 58  class DarcyFlow(object):
58     SMOOTH="SMOOTH"     SMOOTH="SMOOTH"
59     def __init__(self, domain, useReduced=False, solver="POST", verbose=False, w=1.):     def __init__(self, domain, useReduced=False, solver="POST", verbose=False, w=1.):
60        """        """
61        initializes the Darcy flux problem        initializes the Darcy flux problem.
62    
63        :param domain: domain of the problem        :param domain: domain of the problem
64        :type domain: `Domain`        :type domain: `Domain`
65        :param useReduced: uses reduced oreder on flux and pressure        :param useReduced: uses reduced oreder on flux and pressure
66        :type useReduced: ``bool``        :type useReduced: ``bool``
67        :param solver: solver method        :param solver: solver method
68        :type solver: in [`DarcyFlow.EVAL`, `DarcyFlow.POST',  `DarcyFlow.SMOOTH' ]        :type solver: in [`DarcyFlow.EVAL`, `DarcyFlow.POST`, `DarcyFlow.SMOOTH` ]
69        :param verbose: if ``True`` some information on the iteration progress are printed.        :param verbose: if ``True`` some information on the iteration progress are printed.
70        :type verbose: ``bool``        :type verbose: ``bool``
71        :param w: weighting factor for `DarcyFlow.POST` solver        :param w: weighting factor for `DarcyFlow.POST` solver
# Line 71  class DarcyFlow(object): Line 73  class DarcyFlow(object):
73                
74        """        """
75        if not solver in [DarcyFlow.EVAL, DarcyFlow.POST,  DarcyFlow.SMOOTH ] :        if not solver in [DarcyFlow.EVAL, DarcyFlow.POST,  DarcyFlow.SMOOTH ] :
76            raise ValueError,"unknown solver %d."%solver            raise ValueError("unknown solver %d."%solver)
77    
78        self.domain=domain        self.domain=domain
79        self.solver=solver        self.solver=solver
# Line 86  class DarcyFlow(object): Line 88  class DarcyFlow(object):
88    
89        if self.solver  == self.EVAL:        if self.solver  == self.EVAL:
90           self.__pde_v=None           self.__pde_v=None
91       if self.verbose: print "DarcyFlow: simple solver is used."           if self.verbose: print("DarcyFlow: simple solver is used.")
92    
93        elif self.solver  == self.POST:        elif self.solver  == self.POST:
94       if util.inf(w)<0.:           if util.inf(w)<0.:
95          raise ValueError,"Weighting factor must be non-negative."              raise ValueError("Weighting factor must be non-negative.")
96       if self.verbose: print "DarcyFlow: global postprocessing of flux is used."           if self.verbose: print("DarcyFlow: global postprocessing of flux is used.")
97           self.__pde_v=LinearPDESystem(domain)           self.__pde_v=LinearPDESystem(domain)
98           self.__pde_v.setSymmetryOn()           self.__pde_v.setSymmetryOn()
99           if self.useReduced: self.__pde_v.setReducedOrderOn()           if self.useReduced: self.__pde_v.setReducedOrderOn()
100       self.w=w           self.w=w
101           self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale           x=self.domain.getX()
102             self.l=min( [util.sup(x[i])-util.inf(x[i]) for i in xrange(self.domain.getDim()) ] )
103             #self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale
104    
105        elif self.solver  == self.SMOOTH:        elif self.solver  == self.SMOOTH:
106           self.__pde_v=LinearPDESystem(domain)           self.__pde_v=LinearPDESystem(domain)
107           self.__pde_v.setSymmetryOn()           self.__pde_v.setSymmetryOn()
108           if self.useReduced: self.__pde_v.setReducedOrderOn()           if self.useReduced: self.__pde_v.setReducedOrderOn()
109       if self.verbose: print "DarcyFlow: flux smoothing is used."           if self.verbose: print("DarcyFlow: flux smoothing is used.")
110       self.w=0           self.w=0
111    
112        self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))        self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
113        self.__g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))        self.__g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
114          self.__permeability_invXg=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
115          self.__permeability_invXg_ref=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)
116          self.ref_point_id=None
117          self.ref_point=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)
118        self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))        self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
119        self.location_of_fixed_flux = escript.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))        self.location_of_fixed_flux = escript.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
120        self.perm_scale=1.        self.perm_scale=1.
# Line 135  class DarcyFlow(object): Line 143  class DarcyFlow(object):
143    
144        """        """
145        if location_of_fixed_pressure!=None:        if location_of_fixed_pressure!=None:
146             self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)             self.location_of_fixed_pressure=util.wherePositive(util.interpolate(location_of_fixed_pressure, self.__pde_p.getFunctionSpaceForCoefficient("q")))
147               self.ref_point_id=self.location_of_fixed_pressure.maxGlobalDataPoint()
148               if not self.location_of_fixed_pressure.getTupleForGlobalDataPoint(*self.ref_point_id)[0] > 0: raise ValueError("pressure needs to be fixed at least one point.")
149               self.ref_point=self.__pde_p.getFunctionSpaceForCoefficient("q").getX().getTupleForGlobalDataPoint(*self.ref_point_id)
150               if self.verbose: print(("DarcyFlow: reference point at %s."%(self.ref_point,)))
151             self.__pde_p.setValue(q=self.location_of_fixed_pressure)             self.__pde_p.setValue(q=self.location_of_fixed_pressure)
152        if location_of_fixed_flux!=None:        if location_of_fixed_flux!=None:
153            self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)            self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
154            if not self.__pde_v == None: self.__pde_v.setValue(q=self.location_of_fixed_flux)            if not self.__pde_v == None:
155                  self.__pde_v.setValue(q=self.location_of_fixed_flux)
156                            
157        if permeability!=None:        if permeability!=None:
158            
159       perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))           perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
160           self.perm_scale=util.Lsup(util.length(perm))           self.perm_scale=util.Lsup(util.length(perm))
161       if self.verbose: print "DarcyFlow: permeability scaling factor = %e."%self.perm_scale           if self.verbose: print(("DarcyFlow: permeability scaling factor = %e."%self.perm_scale))
162           perm=perm*(1./self.perm_scale)           perm=perm*(1./self.perm_scale)
163                    
164       if perm.getRank()==0:           if perm.getRank()==0:
165    
166          perm_inv=(1./perm)              perm_inv=(1./perm)
167          perm_inv=perm_inv*util.kronecker(self.domain.getDim())              perm_inv=perm_inv*util.kronecker(self.domain.getDim())
168          perm=perm*util.kronecker(self.domain.getDim())              perm=perm*util.kronecker(self.domain.getDim())
169                    
170                    
171       elif perm.getRank()==2:           elif perm.getRank()==2:
172          perm_inv=util.inverse(perm)              perm_inv=util.inverse(perm)
173       else:           else:
174          raise ValueError,"illegal rank of permeability."              raise ValueError("illegal rank of permeability.")
175                    
176       self.__permeability=perm           self.__permeability=perm
177       self.__permeability_inv=perm_inv           self.__permeability_inv=perm_inv
178            
179           #====================           #====================
180       self.__pde_p.setValue(A=self.__permeability)           self.__pde_p.setValue(A=self.__permeability)
181           if self.solver  == self.EVAL:           if self.solver  == self.EVAL:
182                pass # no extra work required                pass # no extra work required
183           elif self.solver  == self.POST:           elif self.solver  == self.POST:
184          k=util.kronecker(self.domain.getDim())                k=util.kronecker(self.domain.getDim())
185          self.omega = self.w*util.length(perm_inv)*self.l*self.domain.getSize()                self.omega = self.w*util.length(perm_inv)*self.l*self.domain.getSize()
186          self.__pde_v.setValue(D=self.__permeability_inv, A=self.omega*util.outer(k,k))                #self.__pde_v.setValue(D=self.__permeability_inv, A=self.omega*util.outer(k,k))
187                  self.__pde_v.setValue(D=self.__permeability_inv, A_reduced=self.omega*util.outer(k,k))
188           elif self.solver  == self.SMOOTH:           elif self.solver  == self.SMOOTH:
189          self.__pde_v.setValue(D=self.__permeability_inv)              self.__pde_v.setValue(D=self.__permeability_inv)
190    
191        if g != None:        if g != None:
192      g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))          g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
193      if g.isEmpty():          if g.isEmpty():
194            g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))               g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
195      else:          else:
196          if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g"               if not g.getShape()==(self.domain.getDim(),): raise ValueError("illegal shape of g")
197      self.__g=g          self.__g=g
198            self.__permeability_invXg=util.tensor_mult(self.__permeability_inv,self.__g * (1./self.perm_scale ))
199            self.__permeability_invXg_ref=util.integrate(self.__permeability_invXg)/util.vol(self.domain)
200        if f !=None:        if f !=None:
201       f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))           f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
202       if f.isEmpty():                 if f.isEmpty():      
203            f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))               f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
204       else:           else:
205           if f.getRank()>0: raise ValueError,"illegal rank of f."               if f.getRank()>0: raise ValueError("illegal rank of f.")
206       self.__f=f           self.__f=f
207    
208     def getSolverOptionsFlux(self):     def getSolverOptionsFlux(self):
209        """        """
# Line 237  class DarcyFlow(object): Line 253  class DarcyFlow(object):
253        :rtype: ``tuple`` of `escript.Data`.        :rtype: ``tuple`` of `escript.Data`.
254    
255        """        """
256        self.__pde_p.setValue(X=self.__g * 1./self.perm_scale,        p0=util.interpolate(p0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
257          if self.ref_point_id == None:
258              p_ref=0
259          else:
260              p_ref=p0.getTupleForGlobalDataPoint(*self.ref_point_id)[0]
261          p0_hydrostatic=p_ref+util.inner(self.__permeability_invXg_ref, self.__pde_p.getFunctionSpaceForCoefficient("q").getX() - self.ref_point)
262          g_2=self.__g - util.tensor_mult(self.__permeability, self.__permeability_invXg_ref * self.perm_scale)
263          self.__pde_p.setValue(X=g_2 * 1./self.perm_scale,
264                              Y=self.__f * 1./self.perm_scale,                              Y=self.__f * 1./self.perm_scale,
265                              y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux * 1./self.perm_scale ),                              y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux * 1./self.perm_scale ),
266                              r=p0)                              r=p0 - p0_hydrostatic)
267        p=self.__pde_p.getSolution()        pp=self.__pde_p.getSolution()
268        u = self.getFlux(p, u0)        u = self._getFlux(pp, u0)
269        return u,p        return u,pp + p0_hydrostatic
270                
271     def getFlux(self,p, u0=None):     def getFlux(self,p, u0=None):
272          """          """
# Line 258  class DarcyFlow(object): Line 281  class DarcyFlow(object):
281          :return: flux          :return: flux
282          :rtype: `escript.Data`          :rtype: `escript.Data`
283          """          """
284            p=util.interpolate(p, self.__pde_p.getFunctionSpaceForCoefficient("q"))
285            if self.ref_point_id == None:
286                p_ref=0
287            else:
288                p_ref=p.getTupleForGlobalDataPoint(*self.ref_point_id)[0]
289            p_hydrostatic=p_ref+util.inner(self.__permeability_invXg_ref, self.__pde_p.getFunctionSpaceForCoefficient("q").getX() - self.ref_point)
290            return self._getFlux(p-p_hydrostatic, u0)
291    
292       def _getFlux(self, pp, u0=None):
293            """
294            returns the flux for a given pressure ``pp`` where the flux is equal to
295            ``u0`` on locations where ``location_of_fixed_flux`` is positive (see
296            `setValue`). Notice that ``g`` is used, see `setValue`.
297    
298            :param pp: pressure.
299            :type pp: scalar value on the domain (i.e. `escript.Data`).
300            :param u0: flux on the locations of the domain marked in ``location_of_fixed_flux``.
301            :type u0: vector values on the domain (i.e. `escript.Data`) or ``None``
302            :return: flux
303            :rtype: `escript.Data`
304            """
305          if self.solver  == self.EVAL:          if self.solver  == self.EVAL:
306             u = self.__g-self.perm_scale * util.tensor_mult(self.__permeability,util.grad(p))             u = self.__g - util.tensor_mult(self.__permeability, self.perm_scale * (util.grad(pp) + self.__permeability_invXg_ref))
307          elif self.solver  == self.POST or self.solver  == self.SMOOTH:          elif self.solver  == self.POST or self.solver  == self.SMOOTH:
308              self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g * 1./self.perm_scale)-util.grad(p))              self.__pde_v.setValue(Y= self.__permeability_invXg - (util.grad(pp) + self.__permeability_invXg_ref))
309                print
310              if u0 == None:              if u0 == None:
311             self.__pde_v.setValue(r=escript.Data())                 self.__pde_v.setValue(r=escript.Data())
312          else:              else:
313                 if not isinstance(u0, escript.Data) : u0 = escript.Vector(u0, escript.Solution(self.domain))                 if not isinstance(u0, escript.Data) : u0 = escript.Vector(u0, escript.Solution(self.domain))
314             self.__pde_v.setValue(r=1./self.perm_scale * u0)                 self.__pde_v.setValue(r=1./self.perm_scale * u0)
315              u= self.__pde_v.getSolution() * self.perm_scale              u= self.__pde_v.getSolution() * self.perm_scale
316      return u          return u
317                
318  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
319       """       """
# Line 313  class StokesProblemCartesian(Homogeneous Line 358  class StokesProblemCartesian(Homogeneous
358    
359           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
360           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
361       self.__pde_proj.setValue(D=1)           self.__pde_proj.setValue(D=1)
362           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
363    
364       def getSolverOptionsVelocity(self):       def getSolverOptionsVelocity(self):
# Line 322  class StokesProblemCartesian(Homogeneous Line 367  class StokesProblemCartesian(Homogeneous
367            
368       :rtype: `SolverOptions`       :rtype: `SolverOptions`
369       """       """
370       return self.__pde_v.getSolverOptions()           return self.__pde_v.getSolverOptions()
371       def setSolverOptionsVelocity(self, options=None):       def setSolverOptionsVelocity(self, options=None):
372           """           """
373       set the solver options for solving the equation for velocity.       set the solver options for solving the equation for velocity.
# Line 336  class StokesProblemCartesian(Homogeneous Line 381  class StokesProblemCartesian(Homogeneous
381       returns the solver options used  solve the equation for pressure.       returns the solver options used  solve the equation for pressure.
382       :rtype: `SolverOptions`       :rtype: `SolverOptions`
383       """       """
384       return self.__pde_prec.getSolverOptions()           return self.__pde_prec.getSolverOptions()
385       def setSolverOptionsPressure(self, options=None):       def setSolverOptionsPressure(self, options=None):
386           """           """
387       set the solver options for solving the equation for pressure.       set the solver options for solving the equation for pressure.
388       :param options: new solver  options       :param options: new solver  options
389       :type options: `SolverOptions`       :type options: `SolverOptions`
390       """       """
391       self.__pde_prec.setSolverOptions(options)           self.__pde_prec.setSolverOptions(options)
392    
393       def setSolverOptionsDiv(self, options=None):       def setSolverOptionsDiv(self, options=None):
394           """           """
# Line 353  class StokesProblemCartesian(Homogeneous Line 398  class StokesProblemCartesian(Homogeneous
398       :param options: new solver options       :param options: new solver options
399       :type options: `SolverOptions`       :type options: `SolverOptions`
400       """       """
401       self.__pde_proj.setSolverOptions(options)           self.__pde_proj.setSolverOptions(options)
402       def getSolverOptionsDiv(self):       def getSolverOptionsDiv(self):
403           """           """
404       returns the solver options for solving the equation to project the divergence of       returns the solver options for solving the equation to project the divergence of
# Line 361  class StokesProblemCartesian(Homogeneous Line 406  class StokesProblemCartesian(Homogeneous
406            
407       :rtype: `SolverOptions`       :rtype: `SolverOptions`
408       """       """
409       return self.__pde_proj.getSolverOptions()           return self.__pde_proj.getSolverOptions()
410    
411       def updateStokesEquation(self, v, p):       def updateStokesEquation(self, v, p):
412           """           """
# Line 388  class StokesProblemCartesian(Homogeneous Line 433  class StokesProblemCartesian(Homogeneous
433              k=util.kronecker(self.domain.getDim())              k=util.kronecker(self.domain.getDim())
434              kk=util.outer(k,k)              kk=util.outer(k,k)
435              self.eta=util.interpolate(eta, escript.Function(self.domain))              self.eta=util.interpolate(eta, escript.Function(self.domain))
436          self.__pde_prec.setValue(D=1/self.eta)              self.__pde_prec.setValue(D=1/self.eta)
437              self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))              self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
438          if restoration_factor!=None:          if restoration_factor!=None:
439              n=self.domain.getNormal()              n=self.domain.getNormal()
# Line 412  class StokesProblemCartesian(Homogeneous Line 457  class StokesProblemCartesian(Homogeneous
457          :param surface_stress: normal surface stress          :param surface_stress: normal surface stress
458          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
459          :param stress: initial stress          :param stress: initial stress
460      :type stress: `Tensor` object on `FunctionSpace` `Function` or similar          :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
461          """          """
462          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
463    
# Line 425  class StokesProblemCartesian(Homogeneous Line 470  class StokesProblemCartesian(Homogeneous
470           :rtype: ``float``           :rtype: ``float``
471           """           """
472           self.__pde_proj.setValue(Y=-util.div(v))           self.__pde_proj.setValue(Y=-util.div(v))
473       self.getSolverOptionsDiv().setTolerance(tol)           self.getSolverOptionsDiv().setTolerance(tol)
474       self.getSolverOptionsDiv().setAbsoluteTolerance(0.)           self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
475           out=self.__pde_proj.getSolution()           out=self.__pde_proj.getSolution()
476           return out           return out
477    
# Line 475  class StokesProblemCartesian(Homogeneous Line 520  class StokesProblemCartesian(Homogeneous
520           """           """
521           self.updateStokesEquation(v,p)           self.updateStokesEquation(v,p)
522           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
523       self.getSolverOptionsVelocity().setTolerance(tol)           self.getSolverOptionsVelocity().setTolerance(tol)
524       self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)           self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
525           if self.__stress.isEmpty():           if self.__stress.isEmpty():
526              self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))              self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
527           else:           else:
# Line 508  class StokesProblemCartesian(Homogeneous Line 553  class StokesProblemCartesian(Homogeneous
553       def solve_prec(self,Bv, tol):       def solve_prec(self,Bv, tol):
554           """           """
555           applies preconditioner for for *BA^{-1}B^** to *Bv*           applies preconditioner for for *BA^{-1}B^** to *Bv*
556           with accuracy `self.getSubProblemTolerance()`           with accuracy ``self.getSubProblemTolerance()``
557    
558           :param Bv: velocity increment           :param Bv: velocity increment
559           :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*           :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
560           :note: boundary conditions on p are zero.           :note: boundary conditions on p are zero.
561           """           """
562           self.__pde_prec.setValue(Y=Bv)           self.__pde_prec.setValue(Y=Bv)
563       self.getSolverOptionsPressure().setTolerance(tol)           self.getSolverOptionsPressure().setTolerance(tol)
564       self.getSolverOptionsPressure().setAbsoluteTolerance(0.)           self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
565           out=self.__pde_prec.getSolution()           out=self.__pde_prec.getSolution()
566           return out           return out

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

  ViewVC Help
Powered by ViewVC 1.1.26