/[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 1481 by artak, Wed Apr 9 00:45:47 2008 UTC revision 3510 by gross, Fri May 13 06:09:49 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          assigns values to model parameters
111    
112          :param f: volumetic sources/sinks
113          :type f: scalar value on the domain (e.g. `escript.Data`)
114          :param g: flux sources/sinks
115          :type g: vector values on the domain (e.g. `escript.Data`)
116          :param location_of_fixed_pressure: mask for locations where pressure is fixed
117          :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
118          :param location_of_fixed_flux:  mask for locations where flux is fixed.
119          :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
120          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
121          :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
122    
123          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
124          :note: at any point on the boundary of the domain the pressure
125                 (``location_of_fixed_pressure`` >0) or the normal component of the
126                 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
127                 is along the *x_i* axis.
128    
129          """
130          if location_of_fixed_pressure!=None:
131               self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)
132               self.__pde_p.setValue(q=self.location_of_fixed_pressure)
133          if location_of_fixed_flux!=None:
134              self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
135              self.__pde_v.setValue(q=self.location_of_fixed_flux)
136          
137                
138          # pressure  is rescaled by the factor 1/self.scale
139          if permeability!=None:
140        
141         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        solves        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    
347          dp=self.__pde_p.getSolution()
348          p=GMRES(dp,
349                  self.__STAB_Aprod,
350              p0,
351              self.__inner,
352              atol=self.__norm(p0+dp)*self.getTolerance() ,
353              rtol=0.,
354              iter_max=max_iter,
355              iter_restart=iter_restart,
356              verbose=self.verbose,P_R=None)
357                
358              u=self.getFlux(p, u0)
359              return u,p
360    
361       def __solve_SYMSTAB(self, u0, p0, max_iter, iter_restart):
362              # p0 is used as an initial guess
363          u=self.getFlux(p0, u0)  
364              self.__pde_p.setValue( Y= self.__f,
365                                     X=  0.5*(self.__g + u - util.tensor_mult(self.__permeability,util.grad(p0)) ),
366                                     y=  -  util.inner(self.domain.getNormal(), u),
367                                     r=escript.Data())
368          dp=self.__pde_p.getSolution()
369          
370          print dp
371              print p0+dp
372              
373          p=GMRES(dp,
374                  self.__SYMSTAB_Aprod,
375              p0,
376              self.__inner,
377              atol=self.__norm(p0+dp)*self.getTolerance() ,
378              rtol=0.,
379              iter_max=max_iter,
380              iter_restart=iter_restart,
381              verbose=self.verbose,P_R=None)
382                
383              u=self.getFlux(p, u0)
384              return u,p
385    
386            -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i     def __L2(self,v):
387             return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))      
388      
389       def __norm(self,r):
390             return util.sqrt(self.__inner(r,r))
391            
392       def __inner(self,r,s):
393             return util.integrate(util.inner(r,s), escript.Function(self.domain))
394            
395       def __STAB_Aprod(self,p):
396          gp=util.grad(p)
397          self.__pde_v.setValue(Y=-0.5*gp,
398                                X=-p*util.kronecker(self.__pde_v.getDomain()),
399                                y= p * self.domain.getNormal(),  
400                                r=escript.Data())
401          u = -self.__pde_v.getSolution()
402          self.__pde_p.setValue(Y=util.div(u),
403                                X=0.5*(u+util.tensor_mult(self.__permeability,gp)),
404                                y=escript.Data(),
405                                r=escript.Data())
406        
407          return  self.__pde_p.getSolution()
408      
409       def __SYMSTAB_Aprod(self,p):
410          gp=util.grad(p)
411          self.__pde_v.setValue(Y=0.5*gp ,
412                                X=escript.Data(),
413                                y=escript.Data(),  
414                                r=escript.Data())
415          u = -self.__pde_v.getSolution()
416          self.__pde_p.setValue(Y=escript.Data(),
417                                X=0.5*(-u+util.tensor_mult(self.__permeability,gp)),
418                                y=escript.Data(),
419                                r=escript.Data())
420        
421          return  self.__pde_p.getSolution()
422          
423    
424    class StokesProblemCartesian(HomogeneousSaddlePointProblem):
425         """
426         solves
427    
428              -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
429                  u_{i,i}=0                  u_{i,i}=0
430    
431            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
432            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
433    
434        if surface_stress is not give 0 is assumed.       if surface_stress is not given 0 is assumed.
435    
436        typical usage:       typical usage:
437    
438              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
439              sp.setTolerance()              sp.setTolerance()
440              sp.initialize(...)              sp.initialize(...)
441              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
442        """       """
443        def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
444             """
445             initialize the Stokes Problem
446    
447             The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
448             LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
449             with macro elements for the pressure.
450    
451             :param domain: domain of the problem.
452             :type domain: `Domain`
453             """
454           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
455           self.domain=domain           self.domain=domain
456           self.vol=util.integrate(1.,Function(self.domain))           self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
457           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_v.setSymmetryOn()
458           self.__pde_u.setSymmetryOn()      
          self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)  
               
459           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
460           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
461           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
462    
463           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
464           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
465         self.__pde_proj.setValue(D=1)
466           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
          self.__pde_proj.setValue(D=1.)  
467    
468        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):       def getSolverOptionsVelocity(self):
469          self.eta=eta           """
470          A =self.__pde_u.createCoefficientOfGeneralPDE("A")       returns the solver options used  solve the equation for velocity.
471      self.__pde_u.setValue(A=Data())      
472          for i in range(self.domain.getDim()):       :rtype: `SolverOptions`
473          for j in range(self.domain.getDim()):       """
474              A[i,j,j,i] += 1.       return self.__pde_v.getSolverOptions()
475              A[i,j,i,j] += 1.       def setSolverOptionsVelocity(self, options=None):
476      self.__pde_prec.setValue(D=1./self.eta)           """
477          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.
478        
479        def B(self,arg):       :param options: new solver  options
480           d=util.div(arg)       :type options: `SolverOptions`
481           self.__pde_proj.setValue(Y=d)       """
482           self.__pde_proj.setTolerance(self.getSubProblemTolerance())           self.__pde_v.setSolverOptions(options)
483           return self.__pde_proj.getSolution(verbose=self.show_details)       def getSolverOptionsPressure(self):
484             """
485        def inner(self,p0,p1):       returns the solver options used  solve the equation for pressure.
486           s0=util.interpolate(p0,Function(self.domain))       :rtype: `SolverOptions`
487           s1=util.interpolate(p1,Function(self.domain))       """
488         return self.__pde_prec.getSolverOptions()
489         def setSolverOptionsPressure(self, options=None):
490             """
491         set the solver options for solving the equation for pressure.
492         :param options: new solver  options
493         :type options: `SolverOptions`
494         """
495         self.__pde_prec.setSolverOptions(options)
496    
497         def setSolverOptionsDiv(self, options=None):
498             """
499         set the solver options for solving the equation to project the divergence of
500         the velocity onto the function space of presure.
501        
502         :param options: new solver options
503         :type options: `SolverOptions`
504         """
505         self.__pde_proj.setSolverOptions(options)
506         def getSolverOptionsDiv(self):
507             """
508         returns the solver options for solving the equation to project the divergence of
509         the velocity onto the function space of presure.
510        
511         :rtype: `SolverOptions`
512         """
513         return self.__pde_proj.getSolverOptions()
514    
515         def updateStokesEquation(self, v, p):
516             """
517             updates the Stokes equation to consider dependencies from ``v`` and ``p``
518             :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
519             """
520             pass
521         def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
522            """
523            assigns new values to the model parameters.
524    
525            :param f: external force
526            :type f: `Vector` object in `FunctionSpace` `Function` or similar
527            :param fixed_u_mask: mask of locations with fixed velocity.
528            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
529            :param eta: viscosity
530            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
531            :param surface_stress: normal surface stress
532            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
533            :param stress: initial stress
534        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
535            """
536            if eta !=None:
537                k=util.kronecker(self.domain.getDim())
538                kk=util.outer(k,k)
539                self.eta=util.interpolate(eta, escript.Function(self.domain))
540            self.__pde_prec.setValue(D=1/self.eta)
541                self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
542            if restoration_factor!=None:
543                n=self.domain.getNormal()
544                self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
545            if fixed_u_mask!=None:
546                self.__pde_v.setValue(q=fixed_u_mask)
547            if f!=None: self.__f=f
548            if surface_stress!=None: self.__surface_stress=surface_stress
549            if stress!=None: self.__stress=stress
550    
551         def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
552            """
553            assigns values to the model parameters
554    
555            :param f: external force
556            :type f: `Vector` object in `FunctionSpace` `Function` or similar
557            :param fixed_u_mask: mask of locations with fixed velocity.
558            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
559            :param eta: viscosity
560            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
561            :param surface_stress: normal surface stress
562            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
563            :param stress: initial stress
564        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
565            """
566            self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
567    
568         def Bv(self,v,tol):
569             """
570             returns inner product of element p and div(v)
571    
572             :param v: a residual
573             :return: inner product of element p and div(v)
574             :rtype: ``float``
575             """
576             self.__pde_proj.setValue(Y=-util.div(v))
577         self.getSolverOptionsDiv().setTolerance(tol)
578         self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
579             out=self.__pde_proj.getSolution()
580             return out
581    
582         def inner_pBv(self,p,Bv):
583             """
584             returns inner product of element p and Bv=-div(v)
585    
586             :param p: a pressure increment
587             :param Bv: a residual
588             :return: inner product of element p and Bv=-div(v)
589             :rtype: ``float``
590             """
591             return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
592    
593         def inner_p(self,p0,p1):
594             """
595             Returns inner product of p0 and p1
596    
597             :param p0: a pressure
598             :param p1: a pressure
599             :return: inner product of p0 and p1
600             :rtype: ``float``
601             """
602             s0=util.interpolate(p0, escript.Function(self.domain))
603             s1=util.interpolate(p1, escript.Function(self.domain))
604           return util.integrate(s0*s1)           return util.integrate(s0*s1)
605    
606        def getStress(self,u):       def norm_v(self,v):
607           mg=util.grad(u)           """
608           return 2.*self.eta*util.symmetric(mg)           returns the norm of v
609    
610        def solve_A(self,u,p):           :param v: a velovity
611           """           :return: norm of v
612           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           :rtype: non-negative ``float``
613           """           """
614           self.__pde_u.setTolerance(self.getSubProblemTolerance())           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
615           self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain))  
616           return  self.__pde_u.getSolution(verbose=self.show_details)  
617         def getDV(self, p, v, tol):
618        def solve_prec(self,p):           """
619           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           return the value for v for a given p (overwrite)
620           self.__pde_prec.setValue(Y=p)  
621           q=self.__pde_prec.getSolution(verbose=self.show_details)           :param p: a pressure
622           return q           :param v: a initial guess for the value v to return.
623        def stoppingcriterium(self,Bv,v,p):           :return: dv given as *Adv=(f-Av-B^*p)*
624            n_r=util.sqrt(self.inner(Bv,Bv))           """
625            n_v=util.Lsup(v)           self.updateStokesEquation(v,p)
626            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)
627            self.iter+=1       self.getSolverOptionsVelocity().setTolerance(tol)
628            if n_r <= self.vol**(1./2.-1./self.domain.getDim())*n_v*self.getTolerance():       self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
629                if self.verbose: print "PCG terminated after %s steps."%self.iter           if self.__stress.isEmpty():
630                return True              self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
631            else:           else:
632                return False              self.__pde_v.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
633        def stoppingcriterium_GMRES(self,norm_r,norm_b):           out=self.__pde_v.getSolution()
634            if self.verbose: print "GMRES step %s: L2(r) = %s, L2(b)*TOL=%s"%(self.iter,norm_r,norm_b*self.getTolerance())           return  out
635            self.iter+=1  
636            if norm_r <= norm_b*self.getTolerance():       def norm_Bv(self,Bv):
637                if self.verbose: print "GMRES terminated after %s steps."%self.iter          """
638                return True          Returns Bv (overwrite).
639            else:  
640                return False          :rtype: equal to the type of p
641            :note: boundary conditions on p should be zero!
642        def stoppingcriterium_MINRES(self,norm_r,norm_Ax):          """
643            if self.verbose: print "MINRES step %s: L2(r) = %s, L2(b)*TOL=%s"%(self.iter,norm_r,norm_Ax*self.getTolerance())          return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
644            self.iter+=1  
645            if norm_r <= norm_Ax*self.getTolerance():       def solve_AinvBt(self,p, tol):
646                if self.verbose: print "MINRES terminated after %s steps."%self.iter           """
647                return True           Solves *Av=B^*p* with accuracy `tol`
648            else:  
649                return False           :param p: a pressure increment
650             :return: the solution of *Av=B^*p*
651             :note: boundary conditions on v should be zero!
652             """
653             self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
654             out=self.__pde_v.getSolution()
655             return  out
656    
657         def solve_prec(self,Bv, tol):
658             """
659             applies preconditioner for for *BA^{-1}B^** to *Bv*
660             with accuracy `self.getSubProblemTolerance()`
661    
662             :param Bv: velocity increment
663             :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
664             :note: boundary conditions on p are zero.
665             """
666             self.__pde_prec.setValue(Y=Bv)
667         self.getSolverOptionsPressure().setTolerance(tol)
668         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
669             out=self.__pde_prec.getSolution()
670             return out

Legend:
Removed from v.1481  
changed lines
  Added in v.3510

  ViewVC Help
Powered by ViewVC 1.1.26