/[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 1519 by artak, Tue Apr 22 03:45:36 2008 UTC revision 3502 by gross, Thu Apr 28 05:06:24 2011 UTC
# Line 1  Line 1 
1  # $Id:$  # -*- coding: utf-8 -*-
2    ########################################################
3  #  #
4  #######################################################  # Copyright (c) 2003-2010 by University of Queensland
5    # Earth Systems Science Computational Center (ESSCC)
6    # http://www.uq.edu.au/esscc
7  #  #
8  #       Copyright 2008 by University of Queensland  # Primary Business: Queensland, Australia
9  #  # Licensed under the Open Software License version 3.0
10  #                http://esscc.uq.edu.au  # http://www.opensource.org/licenses/osl-3.0.php
 #        Primary Business: Queensland, Australia  
 #  Licensed under the Open Software License version 3.0  
 #     http://www.opensource.org/licenses/osl-3.0.php  
 #  
 #######################################################  
11  #  #
12    ########################################################
13    
14    __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15    Earth Systems Science Computational Center (ESSCC)
16    http://www.uq.edu.au/esscc
17    Primary Business: Queensland, Australia"""
18    __license__="""Licensed under the Open Software License version 3.0
19    http://www.opensource.org/licenses/osl-3.0.php"""
20    __url__="https://launchpad.net/escript-finley"
21    
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"
 __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:$"  
34    
35  from escript import *  import escript
36  import util  import util
37  from linearPDEs import LinearPDE  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
38  from pdetools import HomogeneousSaddlePointProblem  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39    
40  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class DarcyFlow(object):
41       """
42       solves the problem
43      
44       *u_i+k_{ij}*p_{,j} = g_i*
45       *u_{i,i} = f*
46      
47       where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
48      
49       :cvar SIMPLE: simple solver
50       :cvar POST: solver using global postprocessing of flux
51       :cvar STAB: solver uses (non-symmetric) stabilization
52       :cvar SYMSTAB: solver uses symmetric stabilization
53       """
54       SIMPLE="SIMPLE"
55       POST="POST"
56       STAB="STAB"
57       SYMSTAB="SYMSTAB"
58       def __init__(self, domain, useReduced=False, solver="SYMSTAB", verbose=False, w=1.):
59          """
60          initializes the Darcy flux problem
61          :param domain: domain of the problem
62          :type domain: `Domain`
63          :param useReduced: uses reduced oreder on flux and pressure
64          :type useReduced: ``bool``
65          :param solver: solver method
66          :type solver: in [`DarcyFlow.SIMPLE`, `DarcyFlow.POST', `DarcyFlow.STAB`, `DarcyFlow.SYMSTAB` ]
67          :param verbose: if ``True`` some information on the iteration progress are printed.
68          :type verbose: ``bool``
69          :param w: weighting factor for `DarcyFlow.POST` solver
70          :type w: ``float``
71          
72          """
73          self.domain=domain
74          self.solver=solver
75          self.useReduced=useReduced
76          self.verbose=verbose
77          self.scale=1.
78          
79          
80          self.__pde_v=LinearPDESystem(domain)
81          self.__pde_v.setSymmetryOn()
82          if self.useReduced: self.__pde_v.setReducedOrderOn()
83          self.__pde_p=LinearSinglePDE(domain)
84          self.__pde_p.setSymmetryOn()
85          if self.useReduced: self.__pde_p.setReducedOrderOn()
86          
87          if self.solver  == self.SIMPLE:
88         if self.verbose: print "DarcyFlow: simple solver is used."
89             self.__pde_v.setValue(D=util.kronecker(self.domain.getDim()))
90          elif self.solver  == self.POST:
91         self.w=w
92         if util.inf(w)<0.:
93            raise ValueError,"Weighting factor must be non-negative."
94         if self.verbose: print "DarcyFlow: global postprocessing of flux is used."
95          elif self.solver  == self.STAB:
96          if self.verbose: print "DarcyFlow: (non-symmetric) stabilization is used."
97          elif  self.solver  == self.SYMSTAB:
98          if self.verbose: print "DarcyFlow: symmetric stabilization is used."
99          else:
100        raise ValueError,"unknown solver %s"%self.solver
101          self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
102          self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
103          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          self.setTolerance()
106        
107            
108       def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
109        """        """
110        solves        assigns values to model parameters
111    
112            -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i        :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         perm=util.interpolate(permeability,self.__pde_v.getFunctionSpaceForCoefficient("A"))
142             V=util.vol(self.domain)
143             l=V**(1./self.domain.getDim())
144            
145         if perm.getRank()==0:
146            perm_inv=(1./perm)
147                self.scale=util.integrate(perm_inv)/V*l
148            perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
149            perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
150            
151            
152         elif perm.getRank()==2:
153            perm_inv=util.inverse(perm)
154                self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)/V)*l
155            perm_inv*=(1./self.scale)
156            perm=perm*self.scale
157         else:
158            raise ValueError,"illegal rank of permeability."
159            
160         self.__permeability=perm
161         self.__permeability_inv=perm_inv
162         if self.verbose: print "DarcyFlow: scaling factor for pressure is %e."%self.scale
163        
164         if self.solver  == self.SIMPLE:
165            self.__pde_p.setValue(A=self.__permeability)
166         elif self.solver  == self.POST:
167            self.__pde_p.setValue(A=self.__permeability)
168            k=util.kronecker(self.domain.getDim())
169            self.lamb = self.w*util.length(perm_inv)*l
170            self.__pde_v.setValue(D=self.__permeability_inv, A=self.lamb*self.domain.getSize()*util.outer(k,k))
171         elif self.solver  == self.STAB:
172            self.__pde_p.setValue(A=0.5*self.__permeability)
173            self.__pde_v.setValue(D=0.5*self.__permeability_inv)
174         elif  self.solver  == self.SYMSTAB:
175            self.__pde_p.setValue(A=0.5*self.__permeability)
176            self.__pde_v.setValue(D=0.5*self.__permeability_inv)
177    
178          if g != None:
179        g=util.interpolate(g, self.__pde_v.getFunctionSpaceForCoefficient("Y"))
180        if g.isEmpty():
181              g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
182        else:
183            if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g"
184        self.__g=g
185          if f !=None:
186         f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
187         if f.isEmpty():      
188              f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
189         else:
190             if f.getRank()>0: raise ValueError,"illegal rank of f."
191         self.__f=f
192       def getSolverOptionsFlux(self):
193          """
194          Returns the solver options used to solve the flux problems
195          :return: `SolverOptions`
196          """
197          return self.__pde_v.getSolverOptions()
198          
199       def setSolverOptionsFlux(self, options=None):
200          """
201          Sets the solver options used to solve the flux problems
202          If ``options`` is not present, the options are reset to default
203          :param options: `SolverOptions`
204          """
205          return self.__pde_v.setSolverOptions(options)
206        
207       def getSolverOptionsPressure(self):
208          """
209          Returns the solver options used to solve the pressure problems
210          :return: `SolverOptions`
211          """
212          return self.__pde_p.getSolverOptions()
213          
214       def setSolverOptionsPressure(self, options=None):
215          """
216          Sets the solver options used to solve the pressure problems
217          If ``options`` is not present, the options are reset to default
218          
219          :param options: `SolverOptions`
220          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
221          """
222          return self.__pde_p.setSolverOptions(options)
223          
224       def setTolerance(self,rtol=1e-4):
225          """
226          sets the relative tolerance ``rtol`` for the pressure for the stabelized solvers.
227          
228          :param rtol: relative tolerance for the pressure
229          :type rtol: non-negative ``float``
230          """
231          if rtol<0:
232         raise ValueError,"Relative tolerance needs to be non-negative."
233          self.__rtol=rtol
234          
235       def getTolerance(self):
236          """
237          returns the relative tolerance
238          :return: current relative tolerance
239          :rtype: ``float``
240          """
241          return self.__rtol
242          
243       def solve(self,u0,p0, max_iter=100, iter_restart=20):
244          """
245          solves the problem.
246          
247          The iteration is terminated if the residual norm is less then self.getTolerance().
248    
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          dp=self.__pde_p.getSolution()
347          p=GMRES(dp,
348                  self.__STAB_Aprod,
349              p0,
350              self.__inner,
351              atol=self.__norm(p0+dp)*self.getTolerance() ,
352              rtol=0.,
353              iter_max=max_iter,
354              iter_restart=iter_restart,
355              verbose=self.verbose,P_R=None)
356                
357              u=self.getFlux(p, u0)
358              return u,p
359    
360       def __solve_SYMSTAB(self, u0, p0, max_iter, iter_restart):
361              # p0 is used as an initial guess
362          u=self.getFlux(p0, u0)  
363              self.__pde_p.setValue( Y= self.__f,
364                                     X=  0.5*(self.__g + u - util.tensor_mult(self.__permeability,util.grad(p0)) ),
365                                     y= -   util.inner(self.domain.getNormal(), u),
366                                     r=escript.Data())
367          dp=self.__pde_p.getSolution()
368          p=GMRES(dp,
369                  self.__SYMSTAB_Aprod,
370              p0,
371              self.__inner,
372              atol=self.__norm(p0+dp)*self.getTolerance() ,
373              rtol=0.,
374              iter_max=max_iter,
375              iter_restart=iter_restart,
376              verbose=self.verbose,P_R=None)
377                
378              u=self.getFlux(p, u0)
379              return u,p
380    
381       def __L2(self,v):
382             return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))      
383      
384       def __norm(self,r):
385             return util.sqrt(self.__inner(r,r))
386            
387       def __inner(self,r,s):
388             return util.integrate(util.inner(r,s), escript.Function(self.domain))
389            
390       def __STAB_Aprod(self,p):
391          gp=util.grad(p)
392          self.__pde_v.setValue(Y=-0.5*gp,
393                                X=-p*util.kronecker(self.__pde_v.getDomain()),
394                                y= p * self.domain.getNormal(),  
395                                r=escript.Data())
396          u = -self.__pde_v.getSolution()
397          self.__pde_p.setValue(Y=util.div(u),
398                                X=0.5*(u+util.tensor_mult(self.__permeability,gp)),
399                                y=escript.Data(),
400                                r=escript.Data())
401        
402          return  self.__pde_p.getSolution()
403      
404       def __SYMSTAB_Aprod(self,p):
405          gp=util.grad(p)
406          self.__pde_v.setValue(Y=0.5*gp ,
407                                X=escript.Data(),
408                                y=escript.Data(),  
409                                r=escript.Data())
410          u = -self.__pde_v.getSolution()
411          self.__pde_p.setValue(Y=escript.Data(),
412                                X=0.5*(-u+util.tensor_mult(self.__permeability,gp)),
413                                y=   util.inner(self.domain.getNormal(), u),
414                                r=escript.Data())
415        
416          return  self.__pde_p.getSolution()
417          
418    
419    class StokesProblemCartesian(HomogeneousSaddlePointProblem):
420         """
421         solves
422    
423              -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
424                  u_{i,i}=0                  u_{i,i}=0
425    
426            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
427            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
428    
429        if surface_stress is not give 0 is assumed.       if surface_stress is not given 0 is assumed.
430    
431        typical usage:       typical usage:
432    
433              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
434              sp.setTolerance()              sp.setTolerance()
435              sp.initialize(...)              sp.initialize(...)
436              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
437        """       """
438        def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
439             """
440             initialize the Stokes Problem
441    
442             The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
443             LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
444             with macro elements for the pressure.
445    
446             :param domain: domain of the problem.
447             :type domain: `Domain`
448             """
449           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
450           self.domain=domain           self.domain=domain
451           self.vol=util.integrate(1.,Function(self.domain))           self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
452           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_v.setSymmetryOn()
453           self.__pde_u.setSymmetryOn()      
          self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)  
               
454           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
455           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
456           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
457    
458           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
459           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
460         self.__pde_proj.setValue(D=1)
461           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
          self.__pde_proj.setValue(D=1.)  
462    
463        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):       def getSolverOptionsVelocity(self):
464          self.eta=eta           """
465          A =self.__pde_u.createCoefficientOfGeneralPDE("A")       returns the solver options used  solve the equation for velocity.
466      self.__pde_u.setValue(A=Data())      
467          for i in range(self.domain.getDim()):       :rtype: `SolverOptions`
468          for j in range(self.domain.getDim()):       """
469              A[i,j,j,i] += 1.       return self.__pde_v.getSolverOptions()
470              A[i,j,i,j] += 1.       def setSolverOptionsVelocity(self, options=None):
471      self.__pde_prec.setValue(D=1./self.eta)           """
472          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.
473        
474        def B(self,arg):       :param options: new solver  options
475           d=util.div(arg)       :type options: `SolverOptions`
476           self.__pde_proj.setValue(Y=d)       """
477           self.__pde_proj.setTolerance(self.getSubProblemTolerance())           self.__pde_v.setSolverOptions(options)
478           return self.__pde_proj.getSolution(verbose=self.show_details)       def getSolverOptionsPressure(self):
479             """
480        def inner(self,p0,p1):       returns the solver options used  solve the equation for pressure.
481           s0=util.interpolate(p0,Function(self.domain))       :rtype: `SolverOptions`
482           s1=util.interpolate(p1,Function(self.domain))       """
483         return self.__pde_prec.getSolverOptions()
484         def setSolverOptionsPressure(self, options=None):
485             """
486         set the solver options for solving the equation for pressure.
487         :param options: new solver  options
488         :type options: `SolverOptions`
489         """
490         self.__pde_prec.setSolverOptions(options)
491    
492         def setSolverOptionsDiv(self, options=None):
493             """
494         set the solver options for solving the equation to project the divergence of
495         the velocity onto the function space of presure.
496        
497         :param options: new solver options
498         :type options: `SolverOptions`
499         """
500         self.__pde_proj.setSolverOptions(options)
501         def getSolverOptionsDiv(self):
502             """
503         returns the solver options for solving the equation to project the divergence of
504         the velocity onto the function space of presure.
505        
506         :rtype: `SolverOptions`
507         """
508         return self.__pde_proj.getSolverOptions()
509    
510         def updateStokesEquation(self, v, p):
511             """
512             updates the Stokes equation to consider dependencies from ``v`` and ``p``
513             :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
514             """
515             pass
516         def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
517            """
518            assigns new values to the model parameters.
519    
520            :param f: external force
521            :type f: `Vector` object in `FunctionSpace` `Function` or similar
522            :param fixed_u_mask: mask of locations with fixed velocity.
523            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
524            :param eta: viscosity
525            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
526            :param surface_stress: normal surface stress
527            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
528            :param stress: initial stress
529        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
530            """
531            if eta !=None:
532                k=util.kronecker(self.domain.getDim())
533                kk=util.outer(k,k)
534                self.eta=util.interpolate(eta, escript.Function(self.domain))
535            self.__pde_prec.setValue(D=1/self.eta)
536                self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
537            if restoration_factor!=None:
538                n=self.domain.getNormal()
539                self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
540            if fixed_u_mask!=None:
541                self.__pde_v.setValue(q=fixed_u_mask)
542            if f!=None: self.__f=f
543            if surface_stress!=None: self.__surface_stress=surface_stress
544            if stress!=None: self.__stress=stress
545    
546         def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
547            """
548            assigns values to the model parameters
549    
550            :param f: external force
551            :type f: `Vector` object in `FunctionSpace` `Function` or similar
552            :param fixed_u_mask: mask of locations with fixed velocity.
553            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
554            :param eta: viscosity
555            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
556            :param surface_stress: normal surface stress
557            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
558            :param stress: initial stress
559        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
560            """
561            self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
562    
563         def Bv(self,v,tol):
564             """
565             returns inner product of element p and div(v)
566    
567             :param v: a residual
568             :return: inner product of element p and div(v)
569             :rtype: ``float``
570             """
571             self.__pde_proj.setValue(Y=-util.div(v))
572         self.getSolverOptionsDiv().setTolerance(tol)
573         self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
574             out=self.__pde_proj.getSolution()
575             return out
576    
577         def inner_pBv(self,p,Bv):
578             """
579             returns inner product of element p and Bv=-div(v)
580    
581             :param p: a pressure increment
582             :param Bv: a residual
583             :return: inner product of element p and Bv=-div(v)
584             :rtype: ``float``
585             """
586             return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
587    
588         def inner_p(self,p0,p1):
589             """
590             Returns inner product of p0 and p1
591    
592             :param p0: a pressure
593             :param p1: a pressure
594             :return: inner product of p0 and p1
595             :rtype: ``float``
596             """
597             s0=util.interpolate(p0, escript.Function(self.domain))
598             s1=util.interpolate(p1, escript.Function(self.domain))
599           return util.integrate(s0*s1)           return util.integrate(s0*s1)
600    
601        def getStress(self,u):       def norm_v(self,v):
602           mg=util.grad(u)           """
603           return 2.*self.eta*util.symmetric(mg)           returns the norm of v
604    
605        def solve_A(self,u,p):           :param v: a velovity
606           """           :return: norm of v
607           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           :rtype: non-negative ``float``
608           """           """
609           self.__pde_u.setTolerance(self.getSubProblemTolerance())           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
610           self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain))  
611           return  self.__pde_u.getSolution(verbose=self.show_details)  
612         def getDV(self, p, v, tol):
613        def solve_prec(self,p):           """
614           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           return the value for v for a given p (overwrite)
          self.__pde_prec.setValue(Y=p)  
          q=self.__pde_prec.getSolution(verbose=self.show_details)  
          return q  
       def stoppingcriterium(self,Bv,v,p):  
           n_r=util.sqrt(self.inner(Bv,Bv))  
           n_v=util.Lsup(v)  
           if self.verbose: print "PCG step %s: L2(div(v)) = %s, Lsup(v)=%s"%(self.iter,n_r,n_v)  
           self.iter+=1  
           if n_r <= self.vol**(1./2.-1./self.domain.getDim())*n_v*self.getTolerance():  
               if self.verbose: print "PCG terminated after %s steps."%self.iter  
               return True  
           else:  
               return False  
       def stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):  
       if TOL==None:  
              TOL=self.getTolerance()  
           if self.verbose: print "%s step %s: L2(r) = %s, L2(b)*TOL=%s"%(solver,self.iter,norm_r,norm_b*TOL)  
           self.iter+=1  
             
           if norm_r <= norm_b*TOL:  
               if self.verbose: print "%s terminated after %s steps."%(solver,self.iter)  
               return True  
           else:  
               return False  
615    
616             :param p: a pressure
617             :param v: a initial guess for the value v to return.
618             :return: dv given as *Adv=(f-Av-B^*p)*
619             """
620             self.updateStokesEquation(v,p)
621             self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
622         self.getSolverOptionsVelocity().setTolerance(tol)
623         self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
624             if self.__stress.isEmpty():
625                self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
626             else:
627                self.__pde_v.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
628             out=self.__pde_v.getSolution()
629             return  out
630    
631         def norm_Bv(self,Bv):
632            """
633            Returns Bv (overwrite).
634    
635            :rtype: equal to the type of p
636            :note: boundary conditions on p should be zero!
637            """
638            return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
639    
640         def solve_AinvBt(self,p, tol):
641             """
642             Solves *Av=B^*p* with accuracy `tol`
643    
644             :param p: a pressure increment
645             :return: the solution of *Av=B^*p*
646             :note: boundary conditions on v should be zero!
647             """
648             self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
649             out=self.__pde_v.getSolution()
650             return  out
651    
652         def solve_prec(self,Bv, tol):
653             """
654             applies preconditioner for for *BA^{-1}B^** to *Bv*
655             with accuracy `self.getSubProblemTolerance()`
656    
657             :param Bv: velocity increment
658             :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
659             :note: boundary conditions on p are zero.
660             """
661             self.__pde_prec.setValue(Y=Bv)
662         self.getSolverOptionsPressure().setTolerance(tol)
663         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
664             out=self.__pde_prec.getSolution()
665             return out

Legend:
Removed from v.1519  
changed lines
  Added in v.3502

  ViewVC Help
Powered by ViewVC 1.1.26