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

Legend:
Removed from v.1673  
changed lines
  Added in v.3852

  ViewVC Help
Powered by ViewVC 1.1.26