/[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 3619 by gross, Wed Oct 5 03:53:34 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 DarcyFlow(object):  class DarcyFlow(object):
44     """     """
# Line 55  class DarcyFlow(object): Line 58  class DarcyFlow(object):
58     SIMPLE="EVAL"     SIMPLE="EVAL"
59     POST="POST"     POST="POST"
60     SMOOTH="SMOOTH"     SMOOTH="SMOOTH"
61     def __init__(self, domain, useReduced=False, solver="EVAL", verbose=False, w=1.):     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 solver: solver method        :param solver: solver method
70        :type solver: in [`DarcyFlow.EVAL`, `DarcyFlow.POST',  `DarcyFlow.SMOOTH' ]        :type solver: in [`DarcyFlow.EVAL`, `DarcyFlow.POST`, `DarcyFlow.SMOOTH` ]
71        :param verbose: if ``True`` some information on the iteration progress are printed.        :param verbose: if ``True`` some information on the iteration progress are printed.
72        :type verbose: ``bool``        :type verbose: ``bool``
73        :param w: weighting factor for `DarcyFlow.POST` solver        :param w: weighting factor for `DarcyFlow.POST` solver
# Line 71  class DarcyFlow(object): Line 75  class DarcyFlow(object):
75                
76        """        """
77        if not solver in [DarcyFlow.EVAL, DarcyFlow.POST,  DarcyFlow.SMOOTH ] :        if not solver in [DarcyFlow.EVAL, DarcyFlow.POST,  DarcyFlow.SMOOTH ] :
78            raise ValueError,"unknown solver %d."%solver            raise ValueError("unknown solver %d."%solver)
79    
80        self.domain=domain        self.domain=domain
81        self.solver=solver        self.solver=solver
# Line 86  class DarcyFlow(object): Line 90  class DarcyFlow(object):
90    
91        if self.solver  == self.EVAL:        if self.solver  == self.EVAL:
92           self.__pde_v=None           self.__pde_v=None
93       if self.verbose: print "DarcyFlow: simple solver is used."           if self.verbose: print("DarcyFlow: simple solver is used.")
94    
95        elif self.solver  == self.POST:        elif self.solver  == self.POST:
96       if util.inf(w)<0.:           if util.inf(w)<0.:
97          raise ValueError,"Weighting factor must be non-negative."              raise ValueError("Weighting factor must be non-negative.")
98       if self.verbose: print "DarcyFlow: global postprocessing of flux is used."           if self.verbose: print("DarcyFlow: global postprocessing of flux is used.")
99           self.__pde_v=LinearPDESystem(domain)           self.__pde_v=LinearPDESystem(domain)
100           self.__pde_v.setSymmetryOn()           self.__pde_v.setSymmetryOn()
101           if self.useReduced: self.__pde_v.setReducedOrderOn()           if self.useReduced: self.__pde_v.setReducedOrderOn()
102       self.w=w           self.w=w
103           self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale           x=self.domain.getX()
104             self.l=min( [util.sup(x[i])-util.inf(x[i]) for i in range(self.domain.getDim()) ] )
105             #self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale
106    
107        elif self.solver  == self.SMOOTH:        elif self.solver  == self.SMOOTH:
108           self.__pde_v=LinearPDESystem(domain)           self.__pde_v=LinearPDESystem(domain)
109           self.__pde_v.setSymmetryOn()           self.__pde_v.setSymmetryOn()
110           if self.useReduced: self.__pde_v.setReducedOrderOn()           if self.useReduced: self.__pde_v.setReducedOrderOn()
111       if self.verbose: print "DarcyFlow: flux smoothing is used."           if self.verbose: print("DarcyFlow: flux smoothing is used.")
112       self.w=0           self.w=0
113    
114        self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))        self.__f=escore.Data(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
115        self.__g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))        self.__g=escore.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
116        self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))        self.__permeability_invXg=escore.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
117        self.location_of_fixed_flux = escript.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))        self.__permeability_invXg_ref=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)
118          self.ref_point_id=None
119          self.ref_point=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)
120          self.location_of_fixed_pressure = escore.Data(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
121          self.location_of_fixed_flux = escore.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
122        self.perm_scale=1.        self.perm_scale=1.
123            
124                    
# Line 135  class DarcyFlow(object): Line 145  class DarcyFlow(object):
145    
146        """        """
147        if location_of_fixed_pressure!=None:        if location_of_fixed_pressure!=None:
148             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")))
149               self.ref_point_id=self.location_of_fixed_pressure.maxGlobalDataPoint()
150               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               self.ref_point=self.__pde_p.getFunctionSpaceForCoefficient("q").getX().getTupleForGlobalDataPoint(*self.ref_point_id)
152               if self.verbose: print(("DarcyFlow: reference point at %s."%(self.ref_point,)))
153             self.__pde_p.setValue(q=self.location_of_fixed_pressure)             self.__pde_p.setValue(q=self.location_of_fixed_pressure)
154        if location_of_fixed_flux!=None:        if location_of_fixed_flux!=None:
155            self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)            self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
156            if not self.__pde_v == None: self.__pde_v.setValue(q=self.location_of_fixed_flux)            if not self.__pde_v == None:
157                  self.__pde_v.setValue(q=self.location_of_fixed_flux)
158                            
159        if permeability!=None:        if permeability!=None:
160            
161       perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))           perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
162           self.perm_scale=util.Lsup(util.length(perm))           self.perm_scale=util.Lsup(util.length(perm))
163             if self.verbose: print(("DarcyFlow: permeability scaling factor = %e."%self.perm_scale))
164           perm=perm*(1./self.perm_scale)           perm=perm*(1./self.perm_scale)
165                    
166       if perm.getRank()==0:           if perm.getRank()==0:
167    
168          perm_inv=(1./perm)              perm_inv=(1./perm)
169          perm_inv=perm_inv*util.kronecker(self.domain.getDim())              perm_inv=perm_inv*util.kronecker(self.domain.getDim())
170          perm=perm*util.kronecker(self.domain.getDim())              perm=perm*util.kronecker(self.domain.getDim())
171                    
172                    
173       elif perm.getRank()==2:           elif perm.getRank()==2:
174          perm_inv=util.inverse(perm)              perm_inv=util.inverse(perm)
175       else:           else:
176          raise ValueError,"illegal rank of permeability."              raise ValueError("illegal rank of permeability.")
177                    
178       self.__permeability=perm           self.__permeability=perm
179       self.__permeability_inv=perm_inv           self.__permeability_inv=perm_inv
180            
181           #====================           #====================
182       self.__pde_p.setValue(A=self.__permeability)           self.__pde_p.setValue(A=self.__permeability)
183           if self.solver  == self.EVAL:           if self.solver  == self.EVAL:
184                pass # no extra work required                pass # no extra work required
185           elif self.solver  == self.POST:           elif self.solver  == self.POST:
186          k=util.kronecker(self.domain.getDim())                k=util.kronecker(self.domain.getDim())
187          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()
188          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))
189                  self.__pde_v.setValue(D=self.__permeability_inv, A_reduced=self.omega*util.outer(k,k))
190           elif self.solver  == self.SMOOTH:           elif self.solver  == self.SMOOTH:
191          self.__pde_v.setValue(D=self.__permeability_inv)              self.__pde_v.setValue(D=self.__permeability_inv)
192    
193        if g != None:        if g != None:
194      g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))          g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
195      if g.isEmpty():          if g.isEmpty():
196            g=Vector(0,self.__pde_p.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_p.getFunctionSpaceForCoefficient("Y"))           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"))               f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
206       else:           else:
207           if f.getRank()>0: raise ValueError,"illegal rank of f."               if f.getRank()>0: raise ValueError("illegal rank of f.")
208       self.__f=f           self.__f=f
209    
210     def getSolverOptionsFlux(self):     def getSolverOptionsFlux(self):
211        """        """
# Line 236  class DarcyFlow(object): Line 255  class DarcyFlow(object):
255        :rtype: ``tuple`` of `escript.Data`.        :rtype: ``tuple`` of `escript.Data`.
256    
257        """        """
258        self.__pde_p.setValue(X=self.__g * 1./self.perm_scale,        p0=util.interpolate(p0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
259                              Y=self.__f,        if self.ref_point_id == None:
260                              y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux),            p_ref=0
261                              r=p0)        else:
262        p=self.__pde_p.getSolution()            p_ref=p0.getTupleForGlobalDataPoint(*self.ref_point_id)[0]
263        u = self.getFlux(p, u0)        p0_hydrostatic=p_ref+util.inner(self.__permeability_invXg_ref, self.__pde_p.getFunctionSpaceForCoefficient("q").getX() - self.ref_point)
264        return u,p        g_2=self.__g - util.tensor_mult(self.__permeability, self.__permeability_invXg_ref * self.perm_scale)
265          self.__pde_p.setValue(X=g_2 * 1./self.perm_scale,
266                                Y=self.__f * 1./self.perm_scale,
267                                y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux * 1./self.perm_scale ),
268                                r=p0 - p0_hydrostatic)
269          pp=self.__pde_p.getSolution()
270          u = self._getFlux(pp, u0)
271          return u,pp + p0_hydrostatic
272                
273     def getFlux(self,p, u0=None):     def getFlux(self,p, u0=None):
274          """          """
# Line 257  class DarcyFlow(object): Line 283  class DarcyFlow(object):
283          :return: flux          :return: flux
284          :rtype: `escript.Data`          :rtype: `escript.Data`
285          """          """
286            p=util.interpolate(p, self.__pde_p.getFunctionSpaceForCoefficient("q"))
287            if self.ref_point_id == None:
288                p_ref=0
289            else:
290                p_ref=p.getTupleForGlobalDataPoint(*self.ref_point_id)[0]
291            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:          if self.solver  == self.EVAL:
308             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))
309          elif self.solver  == self.POST or self.solver  == self.SMOOTH:          elif self.solver  == self.POST or self.solver  == self.SMOOTH:
310              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))
311                print
312              if u0 == None:              if u0 == None:
313             self.__pde_v.setValue(r=escript.Data())                 self.__pde_v.setValue(r=escore.Data())
314          else:              else:
315             self.__pde_v.setValue(r=u0)                 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              u= self.__pde_v.getSolution() * self.perm_scale
318      return u          return u
319                
320  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
321       """       """
# Line 311  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 320  class StokesProblemCartesian(Homogeneous Line 369  class StokesProblemCartesian(Homogeneous
369            
370       :rtype: `SolverOptions`       :rtype: `SolverOptions`
371       """       """
372       return self.__pde_v.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 334  class StokesProblemCartesian(Homogeneous Line 383  class StokesProblemCartesian(Homogeneous
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 351  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 359  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           """           """
# Line 385  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_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)))
440          if restoration_factor!=None:          if restoration_factor!=None:
441              n=self.domain.getNormal()              n=self.domain.getNormal()
# Line 397  class StokesProblemCartesian(Homogeneous Line 446  class StokesProblemCartesian(Homogeneous
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 410  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 423  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 437  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 448  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 473  class StokesProblemCartesian(Homogeneous Line 522  class StokesProblemCartesian(Homogeneous
522           """           """
523           self.updateStokesEquation(v,p)           self.updateStokesEquation(v,p)
524           self.__pde_v.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_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)))
529           else:           else:
# Line 489  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 499  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_v.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_v.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.3619  
changed lines
  Added in v.4154

  ViewVC Help
Powered by ViewVC 1.1.26