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

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

  ViewVC Help
Powered by ViewVC 1.1.26