/[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 3905 by gross, Tue Jun 5 08:33:41 2012 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:$"  
   
 from escript import *  
 import util  
 from linearPDEs import LinearPDE  
 from pdetools import HomogeneousSaddlePointProblem  
34    
35  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  from . import escript
36    from . import util
37    from .linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
38    from .pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
39    
40    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 EVAL: direct pressure gradient evaluation for flux
50       :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*
51                   where *l* is the length scale, *K* is the inverse of the permeability tensor, and *w* is a positive weighting factor.
52       :cvar SMOOTH: global smoothing by solving the PDE *K_{ij} u_j= - p_{,j} + K_{ij} g_j*
53       """
54       EVAL="EVAL"
55       SIMPLE="EVAL"
56       POST="POST"
57       SMOOTH="SMOOTH"
58       def __init__(self, domain, useReduced=False, solver="POST", 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.EVAL`, `DarcyFlow.POST',  `DarcyFlow.SMOOTH' ]
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          if not solver in [DarcyFlow.EVAL, DarcyFlow.POST,  DarcyFlow.SMOOTH ] :
74              raise ValueError("unknown solver %d."%solver)
75    
76          self.domain=domain
77          self.solver=solver
78          self.useReduced=useReduced
79          self.verbose=verbose
80          self.l=None
81          self.w=None
82        
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.EVAL:
88             self.__pde_v=None
89             if self.verbose: print("DarcyFlow: simple solver is used.")
90    
91          elif self.solver  == self.POST:
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             self.__pde_v=LinearPDESystem(domain)
96             self.__pde_v.setSymmetryOn()
97             if self.useReduced: self.__pde_v.setReducedOrderOn()
98             self.w=w
99             x=self.domain.getX()
100             self.l=min( [util.sup(x[i])-util.inf(x[i]) for i in xrange(self.domain.getDim()) ] )
101             #self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale
102    
103          elif self.solver  == self.SMOOTH:
104             self.__pde_v=LinearPDESystem(domain)
105             self.__pde_v.setSymmetryOn()
106             if self.useReduced: self.__pde_v.setReducedOrderOn()
107             if self.verbose: print("DarcyFlow: flux smoothing is used.")
108             self.w=0
109    
110          self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
111          self.__g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
112          self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
113          self.location_of_fixed_flux = escript.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
114          self.perm_scale=1.
115        
116            
117       def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
118        """        """
119        solves        assigns values to model parameters
120    
121            -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i        :param f: volumetic sources/sinks
122          :type f: scalar value on the domain (e.g. `escript.Data`)
123          :param g: flux sources/sinks
124          :type g: vector values on the domain (e.g. `escript.Data`)
125          :param location_of_fixed_pressure: mask for locations where pressure is fixed
126          :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
127          :param location_of_fixed_flux:  mask for locations where flux is fixed.
128          :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
129          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
130          :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
131    
132          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
133          :note: at any point on the boundary of the domain the pressure
134                 (``location_of_fixed_pressure`` >0) or the normal component of the
135                 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
136                 is along the *x_i* axis.
137    
138          """
139          if location_of_fixed_pressure!=None:
140               self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)
141               self.__pde_p.setValue(q=self.location_of_fixed_pressure)
142          if location_of_fixed_flux!=None:
143              self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
144              if not self.__pde_v == None: self.__pde_v.setValue(q=self.location_of_fixed_flux)
145                
146          if permeability!=None:
147        
148             perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
149             self.perm_scale=util.Lsup(util.length(perm))
150             if self.verbose: print(("DarcyFlow: permeability scaling factor = %e."%self.perm_scale))
151             perm=perm*(1./self.perm_scale)
152            
153             if perm.getRank()==0:
154    
155                perm_inv=(1./perm)
156                perm_inv=perm_inv*util.kronecker(self.domain.getDim())
157                perm=perm*util.kronecker(self.domain.getDim())
158            
159            
160             elif perm.getRank()==2:
161                perm_inv=util.inverse(perm)
162             else:
163                raise ValueError("illegal rank of permeability.")
164            
165             self.__permeability=perm
166             self.__permeability_inv=perm_inv
167        
168             #====================
169             self.__pde_p.setValue(A=self.__permeability)
170             if self.solver  == self.EVAL:
171                  pass # no extra work required
172             elif self.solver  == self.POST:
173                  k=util.kronecker(self.domain.getDim())
174                  self.omega = self.w*util.length(perm_inv)*self.l*self.domain.getSize()
175                  #self.__pde_v.setValue(D=self.__permeability_inv, A=self.omega*util.outer(k,k))
176                  self.__pde_v.setValue(D=self.__permeability_inv, A_reduced=self.omega*util.outer(k,k))
177             elif self.solver  == self.SMOOTH:
178                self.__pde_v.setValue(D=self.__permeability_inv)
179    
180          if g != None:
181            g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
182            if g.isEmpty():
183                 g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
184            else:
185                 if not g.getShape()==(self.domain.getDim(),): raise ValueError("illegal shape of g")
186            self.__g=g
187          if f !=None:
188             f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
189             if f.isEmpty():      
190                 f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
191             else:
192                 if f.getRank()>0: raise ValueError("illegal rank of f.")
193             self.__f=f
194    
195       def getSolverOptionsFlux(self):
196          """
197          Returns the solver options used to solve the flux problems
198          :return: `SolverOptions`
199          """
200          if self.__pde_v == None:
201              return None
202          else:
203              return self.__pde_v.getSolverOptions()
204          
205       def setSolverOptionsFlux(self, options=None):
206          """
207          Sets the solver options used to solve the flux problems
208          If ``options`` is not present, the options are reset to default
209          :param options: `SolverOptions`
210          """
211          if not self.__pde_v == None:
212              self.__pde_v.setSolverOptions(options)
213        
214       def getSolverOptionsPressure(self):
215          """
216          Returns the solver options used to solve the pressure problems
217          :return: `SolverOptions`
218          """
219          return self.__pde_p.getSolverOptions()
220          
221       def setSolverOptionsPressure(self, options=None):
222          """
223          Sets the solver options used to solve the pressure problems
224          If ``options`` is not present, the options are reset to default
225          
226          :param options: `SolverOptions`
227          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
228          """
229          return self.__pde_p.setSolverOptions(options)
230          
231       def solve(self, u0, p0):
232          """
233          solves the problem.
234          
235          :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.
236          :type u0: vector value on the domain (e.g. `escript.Data`).
237          :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.
238          :type p0: scalar value on the domain (e.g. `escript.Data`).
239          :return: flux and pressure
240          :rtype: ``tuple`` of `escript.Data`.
241    
242          """
243          self.__pde_p.setValue(X=self.__g * 1./self.perm_scale,
244                                Y=self.__f * 1./self.perm_scale,
245                                y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux * 1./self.perm_scale ),
246                                r=p0)
247          p=self.__pde_p.getSolution()
248          u = self.getFlux(p, u0)
249          return u,p
250          
251       def getFlux(self,p, u0=None):
252            """
253            returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
254            on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
255            Notice that ``g`` is used, see `setValue`.
256    
257            :param p: pressure.
258            :type p: scalar value on the domain (e.g. `escript.Data`).
259            :param u0: flux on the locations of the domain marked be ``location_of_fixed_flux``.
260            :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
261            :return: flux
262            :rtype: `escript.Data`
263            """
264            if self.solver  == self.EVAL:
265               u = self.__g-self.perm_scale * util.tensor_mult(self.__permeability,util.grad(p))
266            elif self.solver  == self.POST or self.solver  == self.SMOOTH:
267                self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g * 1./self.perm_scale)-util.grad(p))
268                if u0 == None:
269                   self.__pde_v.setValue(r=escript.Data())
270                else:
271                   if not isinstance(u0, escript.Data) : u0 = escript.Vector(u0, escript.Solution(self.domain))
272                   self.__pde_v.setValue(r=1./self.perm_scale * u0)
273                   u= self.__pde_v.getSolution() * self.perm_scale
274            return u
275          
276    class StokesProblemCartesian(HomogeneousSaddlePointProblem):
277         """
278         solves
279    
280              -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
281                  u_{i,i}=0                  u_{i,i}=0
282    
283            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
284            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
285    
286        if surface_stress is not give 0 is assumed.       if surface_stress is not given 0 is assumed.
287    
288        typical usage:       typical usage:
289    
290              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
291              sp.setTolerance()              sp.setTolerance()
292              sp.initialize(...)              sp.initialize(...)
293              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
294        """              sp.setStokesEquation(...) # new values for some parameters
295        def __init__(self,domain,**kwargs):              v1,p1=sp.solve(v,p)
296         """
297         def __init__(self,domain,**kwargs):
298             """
299             initialize the Stokes Problem
300    
301             The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
302             LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
303             with macro elements for the pressure.
304    
305             :param domain: domain of the problem.
306             :type domain: `Domain`
307             """
308           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
309           self.domain=domain           self.domain=domain
310           self.vol=util.integrate(1.,Function(self.domain))           self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
311           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_v.setSymmetryOn()
312           self.__pde_u.setSymmetryOn()      
          self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)  
               
313           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
314           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
315           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
316    
317           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
318           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
319             self.__pde_proj.setValue(D=1)
320           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
          self.__pde_proj.setValue(D=1.)  
321    
322        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):       def getSolverOptionsVelocity(self):
323          self.eta=eta           """
324          A =self.__pde_u.createCoefficientOfGeneralPDE("A")       returns the solver options used  solve the equation for velocity.
325      self.__pde_u.setValue(A=Data())      
326          for i in range(self.domain.getDim()):       :rtype: `SolverOptions`
327          for j in range(self.domain.getDim()):       """
328              A[i,j,j,i] += 1.           return self.__pde_v.getSolverOptions()
329              A[i,j,i,j] += 1.       def setSolverOptionsVelocity(self, options=None):
330      self.__pde_prec.setValue(D=1./self.eta)           """
331          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.
332        
333        def B(self,arg):       :param options: new solver  options
334           d=util.div(arg)       :type options: `SolverOptions`
335           self.__pde_proj.setValue(Y=d)       """
336           self.__pde_proj.setTolerance(self.getSubProblemTolerance())           self.__pde_v.setSolverOptions(options)
337           return self.__pde_proj.getSolution(verbose=self.show_details)       def getSolverOptionsPressure(self):
338             """
339        def inner(self,p0,p1):       returns the solver options used  solve the equation for pressure.
340           s0=util.interpolate(p0,Function(self.domain))       :rtype: `SolverOptions`
341           s1=util.interpolate(p1,Function(self.domain))       """
342             return self.__pde_prec.getSolverOptions()
343         def setSolverOptionsPressure(self, options=None):
344             """
345         set the solver options for solving the equation for pressure.
346         :param options: new solver  options
347         :type options: `SolverOptions`
348         """
349             self.__pde_prec.setSolverOptions(options)
350    
351         def setSolverOptionsDiv(self, options=None):
352             """
353         set the solver options for solving the equation to project the divergence of
354         the velocity onto the function space of presure.
355        
356         :param options: new solver options
357         :type options: `SolverOptions`
358         """
359             self.__pde_proj.setSolverOptions(options)
360         def getSolverOptionsDiv(self):
361             """
362         returns the solver options for solving the equation to project the divergence of
363         the velocity onto the function space of presure.
364        
365         :rtype: `SolverOptions`
366         """
367             return self.__pde_proj.getSolverOptions()
368    
369         def updateStokesEquation(self, v, p):
370             """
371             updates the Stokes equation to consider dependencies from ``v`` and ``p``
372             :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values to model parameters.
373             """
374             pass
375         def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
376            """
377            assigns new values to the model parameters.
378    
379            :param f: external force
380            :type f: `Vector` object in `FunctionSpace` `Function` or similar
381            :param fixed_u_mask: mask of locations with fixed velocity.
382            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
383            :param eta: viscosity
384            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
385            :param surface_stress: normal surface stress
386            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
387            :param stress: initial stress
388        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
389            """
390            if eta !=None:
391                k=util.kronecker(self.domain.getDim())
392                kk=util.outer(k,k)
393                self.eta=util.interpolate(eta, escript.Function(self.domain))
394                self.__pde_prec.setValue(D=1/self.eta)
395                self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
396            if restoration_factor!=None:
397                n=self.domain.getNormal()
398                self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
399            if fixed_u_mask!=None:
400                self.__pde_v.setValue(q=fixed_u_mask)
401            if f!=None: self.__f=f
402            if surface_stress!=None: self.__surface_stress=surface_stress
403            if stress!=None: self.__stress=stress
404    
405         def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
406            """
407            assigns values to the model parameters
408    
409            :param f: external force
410            :type f: `Vector` object in `FunctionSpace` `Function` or similar
411            :param fixed_u_mask: mask of locations with fixed velocity.
412            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
413            :param eta: viscosity
414            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
415            :param surface_stress: normal surface stress
416            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
417            :param stress: initial stress
418            :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
419            """
420            self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
421    
422         def Bv(self,v,tol):
423             """
424             returns inner product of element p and div(v)
425    
426             :param v: a residual
427             :return: inner product of element p and div(v)
428             :rtype: ``float``
429             """
430             self.__pde_proj.setValue(Y=-util.div(v))
431             self.getSolverOptionsDiv().setTolerance(tol)
432             self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
433             out=self.__pde_proj.getSolution()
434             return out
435    
436         def inner_pBv(self,p,Bv):
437             """
438             returns inner product of element p and Bv=-div(v)
439    
440             :param p: a pressure increment
441             :param Bv: a residual
442             :return: inner product of element p and Bv=-div(v)
443             :rtype: ``float``
444             """
445             return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
446    
447         def inner_p(self,p0,p1):
448             """
449             Returns inner product of p0 and p1
450    
451             :param p0: a pressure
452             :param p1: a pressure
453             :return: inner product of p0 and p1
454             :rtype: ``float``
455             """
456             s0=util.interpolate(p0, escript.Function(self.domain))
457             s1=util.interpolate(p1, escript.Function(self.domain))
458           return util.integrate(s0*s1)           return util.integrate(s0*s1)
459    
460        def getStress(self,u):       def norm_v(self,v):
461           mg=util.grad(u)           """
462           return 2.*self.eta*util.symmetric(mg)           returns the norm of v
463    
464        def solve_A(self,u,p):           :param v: a velovity
465           """           :return: norm of v
466           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           :rtype: non-negative ``float``
467           """           """
468           self.__pde_u.setTolerance(self.getSubProblemTolerance())           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
469           self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain))  
470           return  self.__pde_u.getSolution(verbose=self.show_details)  
471         def getDV(self, p, v, tol):
472        def solve_prec(self,p):           """
473           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           return the value for v for a given p
474           self.__pde_prec.setValue(Y=p)  
475           q=self.__pde_prec.getSolution(verbose=self.show_details)           :param p: a pressure
476           return q           :param v: a initial guess for the value v to return.
477        def stoppingcriterium(self,Bv,v,p):           :return: dv given as *Adv=(f-Av-B^*p)*
478            n_r=util.sqrt(self.inner(Bv,Bv))           """
479            n_v=util.Lsup(v)           self.updateStokesEquation(v,p)
480            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)
481            self.iter+=1           self.getSolverOptionsVelocity().setTolerance(tol)
482            if n_r <= self.vol**(1./2.-1./self.domain.getDim())*n_v*self.getTolerance():           self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
483                if self.verbose: print "PCG terminated after %s steps."%self.iter           if self.__stress.isEmpty():
484                return True              self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
485            else:           else:
486                return False              self.__pde_v.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
487        def stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):           out=self.__pde_v.getSolution()
488        if TOL==None:           return  out
489               TOL=self.getTolerance()  
490            if self.verbose: print "%s step %s: L2(r) = %s, L2(b)*TOL=%s"%(solver,self.iter,norm_r,norm_b*TOL)       def norm_Bv(self,Bv):
491            self.iter+=1          """
492                      Returns Bv (overwrite).
493            if norm_r <= norm_b*TOL:  
494                if self.verbose: print "%s terminated after %s steps."%(solver,self.iter)          :rtype: equal to the type of p
495                return True          :note: boundary conditions on p should be zero!
496            else:          """
497                return False          return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
498    
499         def solve_AinvBt(self,p, tol):
500             """
501             Solves *Av=B^*p* with accuracy `tol`
502    
503             :param p: a pressure increment
504             :return: the solution of *Av=B^*p*
505             :note: boundary conditions on v should be zero!
506             """
507             self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
508             out=self.__pde_v.getSolution()
509             return  out
510    
511         def solve_prec(self,Bv, tol):
512             """
513             applies preconditioner for for *BA^{-1}B^** to *Bv*
514             with accuracy `self.getSubProblemTolerance()`
515    
516             :param Bv: velocity increment
517             :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
518             :note: boundary conditions on p are zero.
519             """
520             self.__pde_prec.setValue(Y=Bv)
521             self.getSolverOptionsPressure().setTolerance(tol)
522             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
523             out=self.__pde_prec.getSolution()
524             return out

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

  ViewVC Help
Powered by ViewVC 1.1.26