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

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

  ViewVC Help
Powered by ViewVC 1.1.26