/[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 1517 by artak, Fri Apr 18 02:36:37 2008 UTC revision 3990 by caltinay, Tue Sep 25 05:03:20 2012 UTC
# Line 1  Line 1 
1  # $Id:$  # -*- coding: utf-8 -*-
2    ##############################################################################
3  #  #
4  #######################################################  # Copyright (c) 2003-2012 by University of Queensland
5    # http://www.uq.edu.au
6  #  #
7  #       Copyright 2008 by University of Queensland  # Primary Business: Queensland, Australia
8    # Licensed under the Open Software License version 3.0
9    # http://www.opensource.org/licenses/osl-3.0.php
10  #  #
11  #                http://esscc.uq.edu.au  # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12  #        Primary Business: Queensland, Australia  # Development since 2012 by School of Earth Sciences
 #  Licensed under the Open Software License version 3.0  
 #     http://www.opensource.org/licenses/osl-3.0.php  
 #  
 #######################################################  
13  #  #
14    ##############################################################################
15    
16    __copyright__="""Copyright (c) 2003-2012 by University of Queensland
17    http://www.uq.edu.au
18    Primary Business: Queensland, Australia"""
19    __license__="""Licensed under the Open Software License version 3.0
20    http://www.opensource.org/licenses/osl-3.0.php"""
21    __url__="https://launchpad.net/escript-finley"
22    
23  """  """
24  Some models for flow  Some models for flow
25    
26  @var __author__: name of author  :var __author__: name of author
27  @var __copyright__: copyrights  :var __copyright__: copyrights
28  @var __license__: licence agreement  :var __license__: licence agreement
29  @var __url__: url entry point on documentation  :var __url__: url entry point on documentation
30  @var __version__: version  :var __version__: version
31  @var __date__: date of the version  :var __date__: date of the version
32  """  """
33    
34  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
 __copyright__="""  Copyright (c) 2008 by ACcESS MNRF  
                     http://www.access.edu.au  
                 Primary Business: Queensland, Australia"""  
 __license__="""Licensed under the Open Software License version 3.0  
              http://www.opensource.org/licenses/osl-3.0.php"""  
 __url__="http://www.iservo.edu.au/esys"  
 __version__="$Revision:$"  
 __date__="$Date:$"  
   
 from escript import *  
 import util  
 from linearPDEs import LinearPDE  
 from pdetools import HomogeneousSaddlePointProblem  
35    
36  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  from . import escript
37    from . import util
38    from .linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
39    from .pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
40    
41    class DarcyFlow(object):
42       """
43       solves the problem
44      
45       *u_i+k_{ij}*p_{,j} = g_i*
46       *u_{i,i} = f*
47      
48       where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
49      
50       :cvar EVAL: direct pressure gradient evaluation for flux
51       :cvar POST: global postprocessing of flux by solving the PDE *K_{ij} u_j + (w * K * l u_{k,k})_{,i}= - p_{,j} + K_{ij} g_j*
52                   where *l* is the length scale, *K* is the inverse of the permeability tensor, and *w* is a positive weighting factor.
53       :cvar SMOOTH: global smoothing by solving the PDE *K_{ij} u_j= - p_{,j} + K_{ij} g_j*
54       """
55       EVAL="EVAL"
56       SIMPLE="EVAL"
57       POST="POST"
58       SMOOTH="SMOOTH"
59       def __init__(self, domain, useReduced=False, solver="POST", verbose=False, w=1.):
60          """
61          initializes the Darcy flux problem.
62    
63          :param domain: domain of the problem
64          :type domain: `Domain`
65          :param useReduced: uses reduced oreder on flux and pressure
66          :type useReduced: ``bool``
67          :param solver: solver method
68          :type solver: in [`DarcyFlow.EVAL`, `DarcyFlow.POST`, `DarcyFlow.SMOOTH` ]
69          :param verbose: if ``True`` some information on the iteration progress are printed.
70          :type verbose: ``bool``
71          :param w: weighting factor for `DarcyFlow.POST` solver
72          :type w: ``float``
73          
74          """
75          if not solver in [DarcyFlow.EVAL, DarcyFlow.POST,  DarcyFlow.SMOOTH ] :
76              raise ValueError("unknown solver %d."%solver)
77    
78          self.domain=domain
79          self.solver=solver
80          self.useReduced=useReduced
81          self.verbose=verbose
82          self.l=None
83          self.w=None
84        
85          self.__pde_p=LinearSinglePDE(domain)
86          self.__pde_p.setSymmetryOn()
87          if self.useReduced: self.__pde_p.setReducedOrderOn()
88    
89          if self.solver  == self.EVAL:
90             self.__pde_v=None
91             if self.verbose: print("DarcyFlow: simple solver is used.")
92    
93          elif self.solver  == self.POST:
94             if util.inf(w)<0.:
95                raise ValueError("Weighting factor must be non-negative.")
96             if self.verbose: print("DarcyFlow: global postprocessing of flux is used.")
97             self.__pde_v=LinearPDESystem(domain)
98             self.__pde_v.setSymmetryOn()
99             if self.useReduced: self.__pde_v.setReducedOrderOn()
100             self.w=w
101             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:
106             self.__pde_v=LinearPDESystem(domain)
107             self.__pde_v.setSymmetryOn()
108             if self.useReduced: self.__pde_v.setReducedOrderOn()
109             if self.verbose: print("DarcyFlow: flux smoothing is used.")
110             self.w=0
111    
112          self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
113          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"))
119          self.location_of_fixed_flux = escript.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
120          self.perm_scale=1.
121        
122            
123       def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
124        """        """
125        solves        assigns values to model parameters
126    
127            -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i        :param f: volumetic sources/sinks
128          :type f: scalar value on the domain (e.g. `escript.Data`)
129          :param g: flux sources/sinks
130          :type g: vector values on the domain (e.g. `escript.Data`)
131          :param location_of_fixed_pressure: mask for locations where pressure is fixed
132          :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
133          :param location_of_fixed_flux:  mask for locations where flux is fixed.
134          :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
135          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
136          :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
137    
138          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
139          :note: at any point on the boundary of the domain the pressure
140                 (``location_of_fixed_pressure`` >0) or the normal component of the
141                 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
142                 is along the *x_i* axis.
143    
144          """
145          if location_of_fixed_pressure!=None:
146               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)
152          if location_of_fixed_flux!=None:
153              self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
154              if not self.__pde_v == None:
155                  self.__pde_v.setValue(q=self.location_of_fixed_flux)
156                
157          if permeability!=None:
158        
159             perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
160             self.perm_scale=util.Lsup(util.length(perm))
161             if self.verbose: print(("DarcyFlow: permeability scaling factor = %e."%self.perm_scale))
162             perm=perm*(1./self.perm_scale)
163            
164             if perm.getRank()==0:
165    
166                perm_inv=(1./perm)
167                perm_inv=perm_inv*util.kronecker(self.domain.getDim())
168                perm=perm*util.kronecker(self.domain.getDim())
169            
170            
171             elif perm.getRank()==2:
172                perm_inv=util.inverse(perm)
173             else:
174                raise ValueError("illegal rank of permeability.")
175            
176             self.__permeability=perm
177             self.__permeability_inv=perm_inv
178        
179             #====================
180             self.__pde_p.setValue(A=self.__permeability)
181             if self.solver  == self.EVAL:
182                  pass # no extra work required
183             elif self.solver  == self.POST:
184                  k=util.kronecker(self.domain.getDim())
185                  self.omega = self.w*util.length(perm_inv)*self.l*self.domain.getSize()
186                  #self.__pde_v.setValue(D=self.__permeability_inv, A=self.omega*util.outer(k,k))
187                  self.__pde_v.setValue(D=self.__permeability_inv, A_reduced=self.omega*util.outer(k,k))
188             elif self.solver  == self.SMOOTH:
189                self.__pde_v.setValue(D=self.__permeability_inv)
190    
191          if g != None:
192            g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
193            if g.isEmpty():
194                 g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
195            else:
196                 if not g.getShape()==(self.domain.getDim(),): raise ValueError("illegal shape of g")
197            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:
201             f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
202             if f.isEmpty():      
203                 f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
204             else:
205                 if f.getRank()>0: raise ValueError("illegal rank of f.")
206             self.__f=f
207    
208       def getSolverOptionsFlux(self):
209          """
210          Returns the solver options used to solve the flux problems
211          :return: `SolverOptions`
212          """
213          if self.__pde_v == None:
214              return None
215          else:
216              return self.__pde_v.getSolverOptions()
217          
218       def setSolverOptionsFlux(self, options=None):
219          """
220          Sets the solver options used to solve the flux problems
221          If ``options`` is not present, the options are reset to default
222          :param options: `SolverOptions`
223          """
224          if not self.__pde_v == None:
225              self.__pde_v.setSolverOptions(options)
226        
227       def getSolverOptionsPressure(self):
228          """
229          Returns the solver options used to solve the pressure problems
230          :return: `SolverOptions`
231          """
232          return self.__pde_p.getSolverOptions()
233          
234       def setSolverOptionsPressure(self, options=None):
235          """
236          Sets the solver options used to solve the pressure problems
237          If ``options`` is not present, the options are reset to default
238          
239          :param options: `SolverOptions`
240          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
241          """
242          return self.__pde_p.setSolverOptions(options)
243          
244       def solve(self, u0, p0):
245          """
246          solves the problem.
247          
248          :param u0: initial guess for the flux. At locations in the domain marked by ``location_of_fixed_flux`` the value of ``u0`` is kept unchanged.
249          :type u0: vector value on the domain (e.g. `escript.Data`).
250          :param p0: initial guess for the pressure. At locations in the domain marked by ``location_of_fixed_pressure`` the value of ``p0`` is kept unchanged.
251          :type p0: scalar value on the domain (e.g. `escript.Data`).
252          :return: flux and pressure
253          :rtype: ``tuple`` of `escript.Data`.
254    
255          """
256          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,
265                                y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux * 1./self.perm_scale ),
266                                r=p0 - p0_hydrostatic)
267          pp=self.__pde_p.getSolution()
268          u = self._getFlux(pp, u0)
269          return u,pp + p0_hydrostatic
270          
271       def getFlux(self,p, u0=None):
272            """
273            returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
274            on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
275            Notice that ``g`` is used, see `setValue`.
276    
277            :param p: pressure.
278            :type p: scalar value on the domain (e.g. `escript.Data`).
279            :param u0: flux on the locations of the domain marked be ``location_of_fixed_flux``.
280            :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
281            :return: flux
282            :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:
306               u = self.__g - util.tensor_mult(self.__permeability, self.perm_scale * (util.grad(pp) + self.__permeability_invXg_ref))
307            elif self.solver  == self.POST or self.solver  == self.SMOOTH:
308                self.__pde_v.setValue(Y= self.__permeability_invXg - (util.grad(pp) + self.__permeability_invXg_ref))
309                print
310                if u0 == None:
311                   self.__pde_v.setValue(r=escript.Data())
312                else:
313                   if not isinstance(u0, escript.Data) : u0 = escript.Vector(u0, escript.Solution(self.domain))
314                   self.__pde_v.setValue(r=1./self.perm_scale * u0)
315                u= self.__pde_v.getSolution() * self.perm_scale
316            return u
317          
318    class StokesProblemCartesian(HomogeneousSaddlePointProblem):
319         """
320         solves
321    
322              -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
323                  u_{i,i}=0                  u_{i,i}=0
324    
325            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
326            eta*(u_{i,j}+u_{j,i})*n_j=surface_stress            eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
327    
328        if surface_stress is not give 0 is assumed.       if surface_stress is not given 0 is assumed.
329    
330        typical usage:       typical usage:
331    
332              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
333              sp.setTolerance()              sp.setTolerance()
334              sp.initialize(...)              sp.initialize(...)
335              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
336        """              sp.setStokesEquation(...) # new values for some parameters
337        def __init__(self,domain,**kwargs):              v1,p1=sp.solve(v,p)
338         """
339         def __init__(self,domain,**kwargs):
340             """
341             initialize the Stokes Problem
342    
343             The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
344             LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
345             with macro elements for the pressure.
346    
347             :param domain: domain of the problem.
348             :type domain: `Domain`
349             """
350           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
351           self.domain=domain           self.domain=domain
352           self.vol=util.integrate(1.,Function(self.domain))           self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
353           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_v.setSymmetryOn()
354           self.__pde_u.setSymmetryOn()      
 #         self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)  
               
355           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
356           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
357           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
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)
362           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
          self.__pde_proj.setValue(D=1.)  
363    
364        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):       def getSolverOptionsVelocity(self):
365          self.eta=eta           """
366          A =self.__pde_u.createCoefficientOfGeneralPDE("A")       returns the solver options used  solve the equation for velocity.
367      self.__pde_u.setValue(A=Data())      
368          for i in range(self.domain.getDim()):       :rtype: `SolverOptions`
369          for j in range(self.domain.getDim()):       """
370              A[i,j,j,i] += 1.           return self.__pde_v.getSolverOptions()
371              A[i,j,i,j] += 1.       def setSolverOptionsVelocity(self, options=None):
372      self.__pde_prec.setValue(D=1./self.eta)           """
373          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)       set the solver options for solving the equation for velocity.
374        
375        def B(self,arg):       :param options: new solver  options
376           d=util.div(arg)       :type options: `SolverOptions`
377           self.__pde_proj.setValue(Y=d)       """
378           self.__pde_proj.setTolerance(self.getSubProblemTolerance())           self.__pde_v.setSolverOptions(options)
379           return self.__pde_proj.getSolution(verbose=self.show_details)       def getSolverOptionsPressure(self):
380             """
381        def inner(self,p0,p1):       returns the solver options used  solve the equation for pressure.
382           s0=util.interpolate(p0,Function(self.domain))       :rtype: `SolverOptions`
383           s1=util.interpolate(p1,Function(self.domain))       """
384             return self.__pde_prec.getSolverOptions()
385         def setSolverOptionsPressure(self, options=None):
386             """
387         set the solver options for solving the equation for pressure.
388         :param options: new solver  options
389         :type options: `SolverOptions`
390         """
391             self.__pde_prec.setSolverOptions(options)
392    
393         def setSolverOptionsDiv(self, options=None):
394             """
395         set the solver options for solving the equation to project the divergence of
396         the velocity onto the function space of presure.
397        
398         :param options: new solver options
399         :type options: `SolverOptions`
400         """
401             self.__pde_proj.setSolverOptions(options)
402         def getSolverOptionsDiv(self):
403             """
404         returns the solver options for solving the equation to project the divergence of
405         the velocity onto the function space of presure.
406        
407         :rtype: `SolverOptions`
408         """
409             return self.__pde_proj.getSolverOptions()
410    
411         def updateStokesEquation(self, v, p):
412             """
413             updates the Stokes equation to consider dependencies from ``v`` and ``p``
414             :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values to model parameters.
415             """
416             pass
417         def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
418            """
419            assigns new values to the model parameters.
420    
421            :param f: external force
422            :type f: `Vector` object in `FunctionSpace` `Function` or similar
423            :param fixed_u_mask: mask of locations with fixed velocity.
424            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
425            :param eta: viscosity
426            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
427            :param surface_stress: normal surface stress
428            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
429            :param stress: initial stress
430        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
431            """
432            if eta !=None:
433                k=util.kronecker(self.domain.getDim())
434                kk=util.outer(k,k)
435                self.eta=util.interpolate(eta, escript.Function(self.domain))
436                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)))
438            if restoration_factor!=None:
439                n=self.domain.getNormal()
440                self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
441            if fixed_u_mask!=None:
442                self.__pde_v.setValue(q=fixed_u_mask)
443            if f!=None: self.__f=f
444            if surface_stress!=None: self.__surface_stress=surface_stress
445            if stress!=None: self.__stress=stress
446    
447         def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
448            """
449            assigns values to the model parameters
450    
451            :param f: external force
452            :type f: `Vector` object in `FunctionSpace` `Function` or similar
453            :param fixed_u_mask: mask of locations with fixed velocity.
454            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
455            :param eta: viscosity
456            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
457            :param surface_stress: normal surface stress
458            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
459            :param stress: initial stress
460            :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
461            """
462            self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
463    
464         def Bv(self,v,tol):
465             """
466             returns inner product of element p and div(v)
467    
468             :param v: a residual
469             :return: inner product of element p and div(v)
470             :rtype: ``float``
471             """
472             self.__pde_proj.setValue(Y=-util.div(v))
473             self.getSolverOptionsDiv().setTolerance(tol)
474             self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
475             out=self.__pde_proj.getSolution()
476             return out
477    
478         def inner_pBv(self,p,Bv):
479             """
480             returns inner product of element p and Bv=-div(v)
481    
482             :param p: a pressure increment
483             :param Bv: a residual
484             :return: inner product of element p and Bv=-div(v)
485             :rtype: ``float``
486             """
487             return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
488    
489         def inner_p(self,p0,p1):
490             """
491             Returns inner product of p0 and p1
492    
493             :param p0: a pressure
494             :param p1: a pressure
495             :return: inner product of p0 and p1
496             :rtype: ``float``
497             """
498             s0=util.interpolate(p0, escript.Function(self.domain))
499             s1=util.interpolate(p1, escript.Function(self.domain))
500           return util.integrate(s0*s1)           return util.integrate(s0*s1)
501    
502        def getStress(self,u):       def norm_v(self,v):
503           mg=util.grad(u)           """
504           return 2.*self.eta*util.symmetric(mg)           returns the norm of v
505    
506        def solve_A(self,u,p):           :param v: a velovity
507           """           :return: norm of v
508           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           :rtype: non-negative ``float``
509           """           """
510           self.__pde_u.setTolerance(self.getSubProblemTolerance())           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
511           self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain))  
512           return  self.__pde_u.getSolution(verbose=self.show_details)  
513         def getDV(self, p, v, tol):
514        def solve_prec(self,p):           """
515           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           return the value for v for a given p
516           self.__pde_prec.setValue(Y=p)  
517           q=self.__pde_prec.getSolution(verbose=self.show_details)           :param p: a pressure
518           return q           :param v: a initial guess for the value v to return.
519        def stoppingcriterium(self,Bv,v,p):           :return: dv given as *Adv=(f-Av-B^*p)*
520            n_r=util.sqrt(self.inner(Bv,Bv))           """
521            n_v=util.Lsup(v)           self.updateStokesEquation(v,p)
522            if self.verbose: print "PCG step %s: L2(div(v)) = %s, Lsup(v)=%s"%(self.iter,n_r,n_v)           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
523            self.iter+=1           self.getSolverOptionsVelocity().setTolerance(tol)
524            if n_r <= self.vol**(1./2.-1./self.domain.getDim())*n_v*self.getTolerance():           self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
525                if self.verbose: print "PCG terminated after %s steps."%self.iter           if self.__stress.isEmpty():
526                return True              self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
527            else:           else:
528                return False              self.__pde_v.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
529        def stoppingcriterium2(self,norm_r,norm_b,solver='GMRES'):           out=self.__pde_v.getSolution()
530        if self.verbose: print "%s step %s: L2(r) = %s, L2(b)*TOL=%s"%(solver,self.iter,norm_r,norm_b*self.getTolerance())           return  out
531            self.iter+=1  
532            if norm_r <= norm_b*self.getTolerance():       def norm_Bv(self,Bv):
533                if self.verbose: print "%s terminated after %s steps."%(solver,self.iter)          """
534                return True          Returns Bv (overwrite).
535            else:  
536                return False          :rtype: equal to the type of p
537            :note: boundary conditions on p should be zero!
538            """
539            return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
540    
541         def solve_AinvBt(self,p, tol):
542             """
543             Solves *Av=B^*p* with accuracy `tol`
544    
545             :param p: a pressure increment
546             :return: the solution of *Av=B^*p*
547             :note: boundary conditions on v should be zero!
548             """
549             self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
550             out=self.__pde_v.getSolution()
551             return  out
552    
553         def solve_prec(self,Bv, tol):
554             """
555             applies preconditioner for for *BA^{-1}B^** to *Bv*
556             with accuracy ``self.getSubProblemTolerance()``
557    
558             :param Bv: velocity increment
559             :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
560             :note: boundary conditions on p are zero.
561             """
562             self.__pde_prec.setValue(Y=Bv)
563             self.getSolverOptionsPressure().setTolerance(tol)
564             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
565             out=self.__pde_prec.getSolution()
566             return out

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

  ViewVC Help
Powered by ViewVC 1.1.26