/[escript]/trunk/escript/py_src/flows.py
ViewVC logotype

Diff of /trunk/escript/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 3885 by gross, Wed Apr 4 22:12:26 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):                self.__pde_v.setValue(D=self.__permeability_inv, A_reduced=self.omega*util.outer(k,k))
175           mg=util.grad(u)           elif self.solver  == self.SMOOTH:
176           return 2.*self.eta*util.symmetric(mg)              self.__pde_v.setValue(D=self.__permeability_inv)
177        def getEtaEffective(self):  
178           return self.eta        if g != None:
179            g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
180        def solve_A(self,u,p):          if g.isEmpty():
181           """               g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
182           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)          else:
183           """               if not g.getShape()==(self.domain.getDim(),): raise ValueError("illegal shape of g")
184           self.__pde_u.setTolerance(self.getSubProblemTolerance())          self.__g=g
185           self.__pde_u.setValue(X=-self.getStress(u),X_reduced=-p*util.kronecker(self.domain))        if f !=None:
186           return  self.__pde_u.getSolution(verbose=self.show_details)           f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
187             if f.isEmpty():      
188                 f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
189        def solve_prec(self,p):           else:
190          a=self.__inv_eta*p               if f.getRank()>0: raise ValueError("illegal rank of f.")
191          return a-util.integrate(a)/self.vol           self.__f=f
   
       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  
192    
193       def getSolverOptionsFlux(self):
194          """
195          Returns the solver options used to solve the flux problems
196          :return: `SolverOptions`
197          """
198          if self.__pde_v == None:
199              return None
200          else:
201              return self.__pde_v.getSolverOptions()
202          
203       def setSolverOptionsFlux(self, options=None):
204          """
205          Sets the solver options used to solve the flux problems
206          If ``options`` is not present, the options are reset to default
207          :param options: `SolverOptions`
208          """
209          if not self.__pde_v == None:
210              self.__pde_v.setSolverOptions(options)
211        
212       def getSolverOptionsPressure(self):
213          """
214          Returns the solver options used to solve the pressure problems
215          :return: `SolverOptions`
216          """
217          return self.__pde_p.getSolverOptions()
218          
219       def setSolverOptionsPressure(self, options=None):
220          """
221          Sets the solver options used to solve the pressure problems
222          If ``options`` is not present, the options are reset to default
223          
224          :param options: `SolverOptions`
225          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
226          """
227          return self.__pde_p.setSolverOptions(options)
228          
229       def solve(self, u0, p0):
230          """
231          solves the problem.
232          
233          :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.
234          :type u0: vector value on the domain (e.g. `escript.Data`).
235          :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.
236          :type p0: scalar value on the domain (e.g. `escript.Data`).
237          :return: flux and pressure
238          :rtype: ``tuple`` of `escript.Data`.
239    
 class StokesProblemCartesian(HomogeneousSaddlePointProblem):  
240        """        """
241        solves        self.__pde_p.setValue(X=self.__g * 1./self.perm_scale,
242                                Y=self.__f * 1./self.perm_scale,
243                                y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux * 1./self.perm_scale ),
244                                r=p0)
245          p=self.__pde_p.getSolution()
246          u = self.getFlux(p, u0)
247          return u,p
248          
249       def getFlux(self,p, u0=None):
250            """
251            returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
252            on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
253            Notice that ``g`` is used, see `setValue`.
254    
255            :param p: pressure.
256            :type p: scalar value on the domain (e.g. `escript.Data`).
257            :param u0: flux on the locations of the domain marked be ``location_of_fixed_flux``.
258            :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
259            :return: flux
260            :rtype: `escript.Data`
261            """
262            if self.solver  == self.EVAL:
263               u = self.__g-self.perm_scale * util.tensor_mult(self.__permeability,util.grad(p))
264            elif self.solver  == self.POST or self.solver  == self.SMOOTH:
265                self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g * 1./self.perm_scale)-util.grad(p))
266                if u0 == None:
267                   self.__pde_v.setValue(r=escript.Data())
268                else:
269                   if not isinstance(u0, escript.Data) : u0 = escript.Vector(u0, escript.Solution(self.domain))
270                   self.__pde_v.setValue(r=1./self.perm_scale * u0)
271                   u= self.__pde_v.getSolution() * self.perm_scale
272            return u
273          
274    class StokesProblemCartesian(HomogeneousSaddlePointProblem):
275         """
276         solves
277    
278            -(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}
279                  u_{i,i}=0                  u_{i,i}=0
280    
281            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
282            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
283    
284        if surface_stress is not give 0 is assumed.       if surface_stress is not given 0 is assumed.
285    
286        typical usage:       typical usage:
287    
288              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
289              sp.setTolerance()              sp.setTolerance()
290              sp.initialize(...)              sp.initialize(...)
291              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
292        """              sp.setStokesEquation(...) # new values for some parameters
293        def __init__(self,domain,**kwargs):              v1,p1=sp.solve(v,p)
294         """
295         def __init__(self,domain,**kwargs):
296             """
297             initialize the Stokes Problem
298    
299             The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
300             LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
301             with macro elements for the pressure.
302    
303             :param domain: domain of the problem.
304             :type domain: `Domain`
305             """
306           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
307           self.domain=domain           self.domain=domain
308           self.vol=util.integrate(1.,Function(self.domain))           self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
309           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_v.setSymmetryOn()
310           self.__pde_u.setSymmetryOn()      
          # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)  
               
311           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
312           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
313           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
314    
315           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
316           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
317             self.__pde_proj.setValue(D=1)
318           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
          self.__pde_proj.setValue(D=1.)  
319    
320        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):       def getSolverOptionsVelocity(self):
321          self.eta=eta           """
322          A =self.__pde_u.createCoefficientOfGeneralPDE("A")       returns the solver options used  solve the equation for velocity.
323      self.__pde_u.setValue(A=Data())      
324          for i in range(self.domain.getDim()):       :rtype: `SolverOptions`
325          for j in range(self.domain.getDim()):       """
326              A[i,j,j,i] += 1.           return self.__pde_v.getSolverOptions()
327              A[i,j,i,j] += 1.       def setSolverOptionsVelocity(self, options=None):
328      self.__pde_prec.setValue(D=1/self.eta)           """
329          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.
330        
331        def B(self,arg):       :param options: new solver  options
332           d=util.div(arg)       :type options: `SolverOptions`
333           self.__pde_proj.setValue(Y=d)       """
334           self.__pde_proj.setTolerance(self.getSubProblemTolerance())           self.__pde_v.setSolverOptions(options)
335           return self.__pde_proj.getSolution(verbose=self.show_details)       def getSolverOptionsPressure(self):
336             """
337        def inner(self,p0,p1):       returns the solver options used  solve the equation for pressure.
338           s0=util.interpolate(p0,Function(self.domain))       :rtype: `SolverOptions`
339           s1=util.interpolate(p1,Function(self.domain))       """
340             return self.__pde_prec.getSolverOptions()
341         def setSolverOptionsPressure(self, options=None):
342             """
343         set the solver options for solving the equation for pressure.
344         :param options: new solver  options
345         :type options: `SolverOptions`
346         """
347             self.__pde_prec.setSolverOptions(options)
348    
349         def setSolverOptionsDiv(self, options=None):
350             """
351         set the solver options for solving the equation to project the divergence of
352         the velocity onto the function space of presure.
353        
354         :param options: new solver options
355         :type options: `SolverOptions`
356         """
357             self.__pde_proj.setSolverOptions(options)
358         def getSolverOptionsDiv(self):
359             """
360         returns the solver options for solving the equation to project the divergence of
361         the velocity onto the function space of presure.
362        
363         :rtype: `SolverOptions`
364         """
365             return self.__pde_proj.getSolverOptions()
366    
367         def updateStokesEquation(self, v, p):
368             """
369             updates the Stokes equation to consider dependencies from ``v`` and ``p``
370             :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values to model parameters.
371             """
372             pass
373         def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
374            """
375            assigns new values to the model parameters.
376    
377            :param f: external force
378            :type f: `Vector` object in `FunctionSpace` `Function` or similar
379            :param fixed_u_mask: mask of locations with fixed velocity.
380            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
381            :param eta: viscosity
382            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
383            :param surface_stress: normal surface stress
384            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
385            :param stress: initial stress
386        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
387            """
388            if eta !=None:
389                k=util.kronecker(self.domain.getDim())
390                kk=util.outer(k,k)
391                self.eta=util.interpolate(eta, escript.Function(self.domain))
392                self.__pde_prec.setValue(D=1/self.eta)
393                self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
394            if restoration_factor!=None:
395                n=self.domain.getNormal()
396                self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
397            if fixed_u_mask!=None:
398                self.__pde_v.setValue(q=fixed_u_mask)
399            if f!=None: self.__f=f
400            if surface_stress!=None: self.__surface_stress=surface_stress
401            if stress!=None: self.__stress=stress
402    
403         def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
404            """
405            assigns values to the model parameters
406    
407            :param f: external force
408            :type f: `Vector` object in `FunctionSpace` `Function` or similar
409            :param fixed_u_mask: mask of locations with fixed velocity.
410            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
411            :param eta: viscosity
412            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
413            :param surface_stress: normal surface stress
414            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
415            :param stress: initial stress
416            :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
417            """
418            self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
419    
420         def Bv(self,v,tol):
421             """
422             returns inner product of element p and div(v)
423    
424             :param v: a residual
425             :return: inner product of element p and div(v)
426             :rtype: ``float``
427             """
428             self.__pde_proj.setValue(Y=-util.div(v))
429             self.getSolverOptionsDiv().setTolerance(tol)
430             self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
431             out=self.__pde_proj.getSolution()
432             return out
433    
434         def inner_pBv(self,p,Bv):
435             """
436             returns inner product of element p and Bv=-div(v)
437    
438             :param p: a pressure increment
439             :param Bv: a residual
440             :return: inner product of element p and Bv=-div(v)
441             :rtype: ``float``
442             """
443             return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
444    
445         def inner_p(self,p0,p1):
446             """
447             Returns inner product of p0 and p1
448    
449             :param p0: a pressure
450             :param p1: a pressure
451             :return: inner product of p0 and p1
452             :rtype: ``float``
453             """
454             s0=util.interpolate(p0, escript.Function(self.domain))
455             s1=util.interpolate(p1, escript.Function(self.domain))
456           return util.integrate(s0*s1)           return util.integrate(s0*s1)
457    
458        def inner_a(self,a0,a1):       def norm_v(self,v):
459           p0=util.interpolate(a0[1],Function(self.domain))           """
460           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  
461    
462             :param v: a velovity
463             :return: norm of v
464             :rtype: non-negative ``float``
465             """
466             return util.sqrt(util.integrate(util.length(util.grad(v))**2))
467    
468    
469         def getDV(self, p, v, tol):
470             """
471             return the value for v for a given p
472    
473             :param p: a pressure
474             :param v: a initial guess for the value v to return.
475             :return: dv given as *Adv=(f-Av-B^*p)*
476             """
477             self.updateStokesEquation(v,p)
478             self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
479             self.getSolverOptionsVelocity().setTolerance(tol)
480             self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
481             if self.__stress.isEmpty():
482                self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
483             else:
484                self.__pde_v.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
485             out=self.__pde_v.getSolution()
486             return  out
487    
488         def norm_Bv(self,Bv):
489            """
490            Returns Bv (overwrite).
491    
492            :rtype: equal to the type of p
493            :note: boundary conditions on p should be zero!
494            """
495            return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
496    
497         def solve_AinvBt(self,p, tol):
498             """
499             Solves *Av=B^*p* with accuracy `tol`
500    
501             :param p: a pressure increment
502             :return: the solution of *Av=B^*p*
503             :note: boundary conditions on v should be zero!
504             """
505             self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
506             out=self.__pde_v.getSolution()
507             return  out
508    
509         def solve_prec(self,Bv, tol):
510             """
511             applies preconditioner for for *BA^{-1}B^** to *Bv*
512             with accuracy `self.getSubProblemTolerance()`
513    
514             :param Bv: velocity increment
515             :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
516             :note: boundary conditions on p are zero.
517             """
518             self.__pde_prec.setValue(Y=Bv)
519             self.getSolverOptionsPressure().setTolerance(tol)
520             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
521             out=self.__pde_prec.getSolution()
522             return out

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

  ViewVC Help
Powered by ViewVC 1.1.26