/[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 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:$"  
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             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          assigns values to model parameters
120    
121        typical usage:        :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    
             sp=StokesProblemCartesian(domain)  
             sp.setTolerance()  
             sp.initialize(...)  
             v,p=sp.solve(v0,p0)  
138        """        """
139        def __init__(self,domain,**kwargs):        if location_of_fixed_pressure!=None:
140           HomogeneousSaddlePointProblem.__init__(self,**kwargs)             self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)
141           self.domain=domain             self.__pde_p.setValue(q=self.location_of_fixed_pressure)
142           self.vol=util.integrate(1.,Function(self.domain))        if location_of_fixed_flux!=None:
143           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())            self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
144           self.__pde_u.setSymmetryOn()            if not self.__pde_v == None: self.__pde_v.setValue(q=self.location_of_fixed_flux)
145           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)              
146                      if permeability!=None:
147           # self.__pde_proj=LinearPDE(domain,numEquations=1,numSolutions=1)      
148           # self.__pde_proj.setReducedOrderOn()           perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
149           # self.__pde_proj.setSymmetryOn()           self.perm_scale=util.Lsup(util.length(perm))
150           # self.__pde_proj.setSolverMethod(LinearPDE.LUMPING)           if self.verbose: print(("DarcyFlow: permeability scaling factor = %e."%self.perm_scale))
151             perm=perm*(1./self.perm_scale)
152        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):          
153          self.eta=eta           if perm.getRank()==0:
154          A =self.__pde_u.createCoefficientOfGeneralPDE("A")  
155      self.__pde_u.setValue(A=Data())              perm_inv=(1./perm)
156          for i in range(self.domain.getDim()):              perm_inv=perm_inv*util.kronecker(self.domain.getDim())
157          for j in range(self.domain.getDim()):              perm=perm*util.kronecker(self.domain.getDim())
158              A[i,j,j,i] += 1.          
159              A[i,j,i,j] += 1.          
160          # self.__inv_eta=util.interpolate(self.eta,ReducedFunction(self.domain))           elif perm.getRank()==2:
161          self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask,Y=f,y=surface_stress)              perm_inv=util.inverse(perm)
162             else:
163          # self.__pde_proj.setValue(D=1/eta)              raise ValueError("illegal rank of permeability.")
164          # self.__pde_proj.setValue(Y=1.)          
165          # self.__inv_eta=util.interpolate(self.__pde_proj.getSolution(),ReducedFunction(self.domain))           self.__permeability=perm
166          self.__inv_eta=util.interpolate(self.eta,ReducedFunction(self.domain))           self.__permeability_inv=perm_inv
167        
168        def B(self,arg):           #====================
169           a=util.div(arg, ReducedFunction(self.domain))           self.__pde_p.setValue(A=self.__permeability)
170           return a-util.integrate(a)/self.vol           if self.solver  == self.EVAL:
171                  pass # no extra work required
172        def inner(self,p0,p1):           elif self.solver  == self.POST:
173           return util.integrate(p0*p1)                k=util.kronecker(self.domain.getDim())
174           # return util.integrate(1/self.__inv_eta**2*p0*p1)                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        def getStress(self,u):                self.__pde_v.setValue(D=self.__permeability_inv, A_reduced=self.omega*util.outer(k,k))
177           mg=util.grad(u)           elif self.solver  == self.SMOOTH:
178           return 2.*self.eta*util.symmetric(mg)              self.__pde_v.setValue(D=self.__permeability_inv)
179        def getEtaEffective(self):  
180           return self.eta        if g != None:
181            g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
182        def solve_A(self,u,p):          if g.isEmpty():
183           """               g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
184           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)          else:
185           """               if not g.getShape()==(self.domain.getDim(),): raise ValueError("illegal shape of g")
186           self.__pde_u.setTolerance(self.getSubProblemTolerance())          self.__g=g
187           self.__pde_u.setValue(X=-self.getStress(u),X_reduced=-p*util.kronecker(self.domain))        if f !=None:
188           return  self.__pde_u.getSolution(verbose=self.show_details)           f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
189             if f.isEmpty():      
190                 f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
191        def solve_prec(self,p):           else:
192          a=self.__inv_eta*p               if f.getRank()>0: raise ValueError("illegal rank of f.")
193          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  
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    
 class StokesProblemCartesian(HomogeneousSaddlePointProblem):  
242        """        """
243        solves        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            -(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 inner_a(self,a0,a1):       def norm_v(self,v):
461           p0=util.interpolate(a0[1],Function(self.domain))           """
462           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  
463    
464             :param v: a velovity
465             :return: norm of v
466             :rtype: non-negative ``float``
467             """
468             return util.sqrt(util.integrate(util.length(util.grad(v))**2))
469    
470    
471         def getDV(self, p, v, tol):
472             """
473             return the value for v for a given p
474    
475             :param p: a pressure
476             :param v: a initial guess for the value v to return.
477             :return: dv given as *Adv=(f-Av-B^*p)*
478             """
479             self.updateStokesEquation(v,p)
480             self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
481             self.getSolverOptionsVelocity().setTolerance(tol)
482             self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
483             if self.__stress.isEmpty():
484                self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
485             else:
486                self.__pde_v.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
487             out=self.__pde_v.getSolution()
488             return  out
489    
490         def norm_Bv(self,Bv):
491            """
492            Returns Bv (overwrite).
493    
494            :rtype: equal to the type of p
495            :note: boundary conditions on p should be zero!
496            """
497            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.1673  
changed lines
  Added in v.3905

  ViewVC Help
Powered by ViewVC 1.1.26