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

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

  ViewVC Help
Powered by ViewVC 1.1.26