/[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 4326 by caltinay, Wed Mar 20 02:50:46 2013 UTC
# Line 1  Line 1 
1  # $Id:$  # -*- coding: utf-8 -*-
2    ##############################################################################
3  #  #
4  #######################################################  # Copyright (c) 2003-2013 by University of Queensland
5    # http://www.uq.edu.au
6  #  #
7  #       Copyright 2008 by University of Queensland  # Primary Business: Queensland, Australia
8    # Licensed under the Open Software License version 3.0
9    # http://www.opensource.org/licenses/osl-3.0.php
10  #  #
11  #                http://esscc.uq.edu.au  # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12  #        Primary Business: Queensland, Australia  # Development since 2012 by School of Earth Sciences
 #  Licensed under the Open Software License version 3.0  
 #     http://www.opensource.org/licenses/osl-3.0.php  
 #  
 #######################################################  
13  #  #
14    ##############################################################################
15    
16    __copyright__="""Copyright (c) 2003-2013 by University of Queensland
17    http://www.uq.edu.au
18    Primary Business: Queensland, Australia"""
19    __license__="""Licensed under the Open Software License version 3.0
20    http://www.opensource.org/licenses/osl-3.0.php"""
21    __url__="https://launchpad.net/escript-finley"
22    
23  """  """
24  Some models for flow  Some models for flow
25    
26  @var __author__: name of author  :var __author__: name of author
27  @var __copyright__: copyrights  :var __copyright__: copyrights
28  @var __license__: licence agreement  :var __license__: licence agreement
29  @var __url__: url entry point on documentation  :var __url__: url entry point on documentation
30  @var __version__: version  :var __version__: version
31  @var __date__: date of the version  :var __date__: date of the version
32  """  """
33    
34  __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:$"  
   
 from escript import *  
 import util  
 from linearPDEs import LinearPDE  
 from pdetools import HomogeneousSaddlePointProblem,Projector  
35    
36  class StokesProblemCartesian_DC(HomogeneousSaddlePointProblem):  from . import escriptcpp
37    escore=escriptcpp
38    #from . import escript
39    from . import util
40    from .linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
41    from .pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
42    
43    class DarcyFlow(object):
44       """
45       solves the problem
46      
47       *u_i+k_{ij}*p_{,j} = g_i*
48       *u_{i,i} = f*
49      
50       where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
51      
52       :cvar EVAL: direct pressure gradient evaluation for flux
53       :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*
54                   where *l* is the length scale, *K* is the inverse of the permeability tensor, and *w* is a positive weighting factor.
55       :cvar SMOOTH: global smoothing by solving the PDE *K_{ij} u_j= - p_{,j} + K_{ij} g_j*
56       """
57       EVAL="EVAL"
58       SIMPLE="EVAL"
59       POST="POST"
60       SMOOTH="SMOOTH"
61       def __init__(self, domain, useReduced=False, solver="POST", verbose=False, w=1.):
62        """        """
63        solves        initializes the Darcy flux problem.
64    
65            -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i        :param domain: domain of the problem
66                  u_{i,i}=0        :type domain: `Domain`
67          :param useReduced: uses reduced oreder on flux and pressure
68            u=0 where  fixed_u_mask>0        :type useReduced: ``bool``
69            eta*(u_{i,j}+u_{j,i})*n_j=surface_stress        :param solver: solver method
70          :type solver: in [`DarcyFlow.EVAL`, `DarcyFlow.POST`, `DarcyFlow.SMOOTH` ]
71          :param verbose: if ``True`` some information on the iteration progress are printed.
72          :type verbose: ``bool``
73          :param w: weighting factor for `DarcyFlow.POST` solver
74          :type w: ``float``
75          
76          """
77          if not solver in [DarcyFlow.EVAL, DarcyFlow.POST,  DarcyFlow.SMOOTH ] :
78              raise ValueError("unknown solver %d."%solver)
79    
80        if surface_stress is not give 0 is assumed.        self.domain=domain
81          self.solver=solver
82          self.useReduced=useReduced
83          self.verbose=verbose
84          self.l=None
85          self.w=None
86        
87          self.__pde_p=LinearSinglePDE(domain)
88          self.__pde_p.setSymmetryOn()
89          if self.useReduced: self.__pde_p.setReducedOrderOn()
90    
91          if self.solver  == self.EVAL:
92             self.__pde_v=None
93             if self.verbose: print("DarcyFlow: simple solver is used.")
94    
95          elif self.solver  == self.POST:
96             if util.inf(w)<0.:
97                raise ValueError("Weighting factor must be non-negative.")
98             if self.verbose: print("DarcyFlow: global postprocessing of flux is used.")
99             self.__pde_v=LinearPDESystem(domain)
100             self.__pde_v.setSymmetryOn()
101             if self.useReduced: self.__pde_v.setReducedOrderOn()
102             self.w=w
103             x=self.domain.getX()
104             self.l=min( [util.sup(x[i])-util.inf(x[i]) for i in range(self.domain.getDim()) ] )
105             #self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale
106    
107          elif self.solver  == self.SMOOTH:
108             self.__pde_v=LinearPDESystem(domain)
109             self.__pde_v.setSymmetryOn()
110             if self.useReduced: self.__pde_v.setReducedOrderOn()
111             if self.verbose: print("DarcyFlow: flux smoothing is used.")
112             self.w=0
113    
114          self.__f=escore.Data(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
115          self.__g=escore.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
116          self.__permeability_invXg=escore.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
117          self.__permeability_invXg_ref=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)
118          self.ref_point_id=None
119          self.ref_point=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)
120          self.location_of_fixed_pressure = escore.Data(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
121          self.location_of_fixed_flux = escore.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
122          self.perm_scale=1.
123        
124            
125       def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
126          """
127          assigns values to model parameters
128    
129        typical usage:        :param f: volumetic sources/sinks
130          :type f: scalar value on the domain (e.g. `escript.Data`)
131          :param g: flux sources/sinks
132          :type g: vector values on the domain (e.g. `escript.Data`)
133          :param location_of_fixed_pressure: mask for locations where pressure is fixed
134          :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
135          :param location_of_fixed_flux:  mask for locations where flux is fixed.
136          :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
137          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
138          :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
139    
140          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
141          :note: at any point on the boundary of the domain the pressure
142                 (``location_of_fixed_pressure`` >0) or the normal component of the
143                 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
144                 is along the *x_i* axis.
145    
             sp=StokesProblemCartesian(domain)  
             sp.setTolerance()  
             sp.initialize(...)  
             v,p=sp.solve(v0,p0)  
146        """        """
147        def __init__(self,domain,**kwargs):        if location_of_fixed_pressure!=None:
148           HomogeneousSaddlePointProblem.__init__(self,**kwargs)             self.location_of_fixed_pressure=util.wherePositive(util.interpolate(location_of_fixed_pressure, self.__pde_p.getFunctionSpaceForCoefficient("q")))
149           self.domain=domain             self.ref_point_id=self.location_of_fixed_pressure.maxGlobalDataPoint()
150           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.")
151           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)
152           self.__pde_u.setSymmetryOn()             if self.verbose: print(("DarcyFlow: reference point at %s."%(self.ref_point,)))
153           # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)             self.__pde_p.setValue(q=self.location_of_fixed_pressure)
154                      if location_of_fixed_flux!=None:
155           # self.__pde_proj=LinearPDE(domain,numEquations=1,numSolutions=1)            self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
156           # self.__pde_proj.setReducedOrderOn()            if not self.__pde_v == None:
157           # self.__pde_proj.setSymmetryOn()                self.__pde_v.setValue(q=self.location_of_fixed_flux)
158           # self.__pde_proj.setSolverMethod(LinearPDE.LUMPING)              
159          if permeability!=None:
160        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):      
161          self.eta=eta           perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
162          A =self.__pde_u.createCoefficientOfGeneralPDE("A")           self.perm_scale=util.Lsup(util.length(perm))
163      self.__pde_u.setValue(A=Data())           if self.verbose: print(("DarcyFlow: permeability scaling factor = %e."%self.perm_scale))
164          for i in range(self.domain.getDim()):           perm=perm*(1./self.perm_scale)
165          for j in range(self.domain.getDim()):          
166              A[i,j,j,i] += 1.           if perm.getRank()==0:
167              A[i,j,i,j] += 1.  
168          # self.__inv_eta=util.interpolate(self.eta,ReducedFunction(self.domain))              perm_inv=(1./perm)
169          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())
170                perm=perm*util.kronecker(self.domain.getDim())
171          # self.__pde_proj.setValue(D=1/eta)          
172          # self.__pde_proj.setValue(Y=1.)          
173          # self.__inv_eta=util.interpolate(self.__pde_proj.getSolution(),ReducedFunction(self.domain))           elif perm.getRank()==2:
174          self.__inv_eta=util.interpolate(self.eta,ReducedFunction(self.domain))              perm_inv=util.inverse(perm)
175             else:
176        def B(self,arg):              raise ValueError("illegal rank of permeability.")
177           a=util.div(arg, ReducedFunction(self.domain))          
178           return a-util.integrate(a)/self.vol           self.__permeability=perm
179             self.__permeability_inv=perm_inv
180        def inner(self,p0,p1):      
181           return util.integrate(p0*p1)           #====================
182           # return util.integrate(1/self.__inv_eta**2*p0*p1)           self.__pde_p.setValue(A=self.__permeability)
183             if self.solver  == self.EVAL:
184        def getStress(self,u):                pass # no extra work required
185           mg=util.grad(u)           elif self.solver  == self.POST:
186           return 2.*self.eta*util.symmetric(mg)                k=util.kronecker(self.domain.getDim())
187        def getEtaEffective(self):                self.omega = self.w*util.length(perm_inv)*self.l*self.domain.getSize()
188           return self.eta                #self.__pde_v.setValue(D=self.__permeability_inv, A=self.omega*util.outer(k,k))
189                  self.__pde_v.setValue(D=self.__permeability_inv, A_reduced=self.omega*util.outer(k,k))
190        def solve_A(self,u,p):           elif self.solver  == self.SMOOTH:
191           """              self.__pde_v.setValue(D=self.__permeability_inv)
192           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)  
193           """        if g != None:
194           self.__pde_u.setTolerance(self.getSubProblemTolerance())          g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
195           self.__pde_u.setValue(X=-self.getStress(u),X_reduced=-p*util.kronecker(self.domain))          if g.isEmpty():
196           return  self.__pde_u.getSolution(verbose=self.show_details)               g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
197            else:
198                 if not g.getShape()==(self.domain.getDim(),): raise ValueError("illegal shape of g")
199        def solve_prec(self,p):          self.__g=g
200          a=self.__inv_eta*p          self.__permeability_invXg=util.tensor_mult(self.__permeability_inv,self.__g * (1./self.perm_scale ))
201          return a-util.integrate(a)/self.vol          self.__permeability_invXg_ref=util.integrate(self.__permeability_invXg)/util.vol(self.domain)
202          if f !=None:
203        def stoppingcriterium(self,Bv,v,p):           f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
204            n_r=util.sqrt(self.inner(Bv,Bv))           if f.isEmpty():      
205            n_v=util.sqrt(util.integrate(util.length(util.grad(v))**2))               f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
206            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:
207            if self.iter == 0: self.__n_v=n_v;               if f.getRank()>0: raise ValueError("illegal rank of f.")
208            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  
209    
210       def getSolverOptionsFlux(self):
211          """
212          Returns the solver options used to solve the flux problems
213          :return: `SolverOptions`
214          """
215          if self.__pde_v == None:
216              return None
217          else:
218              return self.__pde_v.getSolverOptions()
219          
220       def setSolverOptionsFlux(self, options=None):
221          """
222          Sets the solver options used to solve the flux problems
223          If ``options`` is not present, the options are reset to default
224          :param options: `SolverOptions`
225          """
226          if not self.__pde_v == None:
227              self.__pde_v.setSolverOptions(options)
228        
229       def getSolverOptionsPressure(self):
230          """
231          Returns the solver options used to solve the pressure problems
232          :return: `SolverOptions`
233          """
234          return self.__pde_p.getSolverOptions()
235          
236       def setSolverOptionsPressure(self, options=None):
237          """
238          Sets the solver options used to solve the pressure problems
239          If ``options`` is not present, the options are reset to default
240          
241          :param options: `SolverOptions`
242          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
243          """
244          return self.__pde_p.setSolverOptions(options)
245          
246       def solve(self, u0, p0):
247          """
248          solves the problem.
249          
250          :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.
251          :type u0: vector value on the domain (e.g. `escript.Data`).
252          :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.
253          :type p0: scalar value on the domain (e.g. `escript.Data`).
254          :return: flux and pressure
255          :rtype: ``tuple`` of `escript.Data`.
256    
 class StokesProblemCartesian(HomogeneousSaddlePointProblem):  
257        """        """
258        solves        p0=util.interpolate(p0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
259          if self.ref_point_id == None:
260              p_ref=0
261          else:
262              p_ref=p0.getTupleForGlobalDataPoint(*self.ref_point_id)[0]
263          p0_hydrostatic=p_ref+util.inner(self.__permeability_invXg_ref, self.__pde_p.getFunctionSpaceForCoefficient("q").getX() - self.ref_point)
264          g_2=self.__g - util.tensor_mult(self.__permeability, self.__permeability_invXg_ref * self.perm_scale)
265          self.__pde_p.setValue(X=g_2 * 1./self.perm_scale,
266                                Y=self.__f * 1./self.perm_scale,
267                                y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux * 1./self.perm_scale ),
268                                r=p0 - p0_hydrostatic)
269          pp=self.__pde_p.getSolution()
270          u = self._getFlux(pp, u0)
271          return u,pp + p0_hydrostatic
272          
273       def getFlux(self,p, u0=None):
274            """
275            returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
276            on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
277            Notice that ``g`` is used, see `setValue`.
278    
279            :param p: pressure.
280            :type p: scalar value on the domain (e.g. `escript.Data`).
281            :param u0: flux on the locations of the domain marked be ``location_of_fixed_flux``.
282            :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
283            :return: flux
284            :rtype: `escript.Data`
285            """
286            p=util.interpolate(p, self.__pde_p.getFunctionSpaceForCoefficient("q"))
287            if self.ref_point_id == None:
288                p_ref=0
289            else:
290                p_ref=p.getTupleForGlobalDataPoint(*self.ref_point_id)[0]
291            p_hydrostatic=p_ref+util.inner(self.__permeability_invXg_ref, self.__pde_p.getFunctionSpaceForCoefficient("q").getX() - self.ref_point)
292            return self._getFlux(p-p_hydrostatic, u0)
293    
294       def _getFlux(self, pp, u0=None):
295            """
296            returns the flux for a given pressure ``pp`` where the flux is equal to
297            ``u0`` on locations where ``location_of_fixed_flux`` is positive (see
298            `setValue`). Notice that ``g`` is used, see `setValue`.
299    
300            :param pp: pressure.
301            :type pp: scalar value on the domain (i.e. `escript.Data`).
302            :param u0: flux on the locations of the domain marked in ``location_of_fixed_flux``.
303            :type u0: vector values on the domain (i.e. `escript.Data`) or ``None``
304            :return: flux
305            :rtype: `escript.Data`
306            """
307            if self.solver  == self.EVAL:
308               u = self.__g - util.tensor_mult(self.__permeability, self.perm_scale * (util.grad(pp) + self.__permeability_invXg_ref))
309            elif self.solver  == self.POST or self.solver  == self.SMOOTH:
310                self.__pde_v.setValue(Y= self.__permeability_invXg - (util.grad(pp) + self.__permeability_invXg_ref))
311    
312                if u0 == None:
313                   self.__pde_v.setValue(r=escore.Data())
314                else:
315                   if not isinstance(u0, escore.Data) : u0 = escore.Vector(u0, escore.Solution(self.domain))
316                   self.__pde_v.setValue(r=1./self.perm_scale * u0)
317                u= self.__pde_v.getSolution() * self.perm_scale
318            return u
319          
320    class StokesProblemCartesian(HomogeneousSaddlePointProblem):
321         """
322         solves
323    
324            -(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}
325                  u_{i,i}=0                  u_{i,i}=0
326    
327            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
328            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
329    
330        if surface_stress is not give 0 is assumed.       if surface_stress is not given 0 is assumed.
331    
332        typical usage:       typical usage:
333    
334              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
335              sp.setTolerance()              sp.setTolerance()
336              sp.initialize(...)              sp.initialize(...)
337              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
338        """              sp.setStokesEquation(...) # new values for some parameters
339        def __init__(self,domain,**kwargs):              v1,p1=sp.solve(v,p)
340         """
341         def __init__(self,domain,**kwargs):
342             """
343             initialize the Stokes Problem
344    
345             The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
346             LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
347             with macro elements for the pressure.
348    
349             :param domain: domain of the problem.
350             :type domain: `Domain`
351             """
352           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
353           self.domain=domain           self.domain=domain
354           self.vol=util.integrate(1.,Function(self.domain))           self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
355           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_v.setSymmetryOn()
356           self.__pde_u.setSymmetryOn()      
          # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)  
               
357           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
358           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
359           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
360    
361           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
362           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
363             self.__pde_proj.setValue(D=1)
364           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
          self.__pde_proj.setValue(D=1.)  
365    
366        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data()):       def getSolverOptionsVelocity(self):
367          self.eta=eta           """
368          A =self.__pde_u.createCoefficientOfGeneralPDE("A")       returns the solver options used  solve the equation for velocity.
369      self.__pde_u.setValue(A=Data())      
370          for i in range(self.domain.getDim()):       :rtype: `SolverOptions`
371          for j in range(self.domain.getDim()):       """
372              A[i,j,j,i] += 1.           return self.__pde_v.getSolverOptions()
373              A[i,j,i,j] += 1.       def setSolverOptionsVelocity(self, options=None):
374      self.__pde_prec.setValue(D=1/self.eta)           """
375          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.
376        
377        def B(self,arg):       :param options: new solver  options
378           d=util.div(arg)       :type options: `SolverOptions`
379           self.__pde_proj.setValue(Y=d)       """
380           self.__pde_proj.setTolerance(self.getSubProblemTolerance())           self.__pde_v.setSolverOptions(options)
381           return self.__pde_proj.getSolution(verbose=self.show_details)       def getSolverOptionsPressure(self):
382             """
383        def inner(self,p0,p1):       returns the solver options used  solve the equation for pressure.
384           s0=util.interpolate(p0,Function(self.domain))       :rtype: `SolverOptions`
385           s1=util.interpolate(p1,Function(self.domain))       """
386             return self.__pde_prec.getSolverOptions()
387         def setSolverOptionsPressure(self, options=None):
388             """
389         set the solver options for solving the equation for pressure.
390         :param options: new solver  options
391         :type options: `SolverOptions`
392         """
393             self.__pde_prec.setSolverOptions(options)
394    
395         def setSolverOptionsDiv(self, options=None):
396             """
397         set the solver options for solving the equation to project the divergence of
398         the velocity onto the function space of presure.
399        
400         :param options: new solver options
401         :type options: `SolverOptions`
402         """
403             self.__pde_proj.setSolverOptions(options)
404         def getSolverOptionsDiv(self):
405             """
406         returns the solver options for solving the equation to project the divergence of
407         the velocity onto the function space of presure.
408        
409         :rtype: `SolverOptions`
410         """
411             return self.__pde_proj.getSolverOptions()
412    
413         def updateStokesEquation(self, v, p):
414             """
415             updates the Stokes equation to consider dependencies from ``v`` and ``p``
416             :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values to model parameters.
417             """
418             pass
419         def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
420            """
421            assigns new values to the model parameters.
422    
423            :param f: external force
424            :type f: `Vector` object in `FunctionSpace` `Function` or similar
425            :param fixed_u_mask: mask of locations with fixed velocity.
426            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
427            :param eta: viscosity
428            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
429            :param surface_stress: normal surface stress
430            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
431            :param stress: initial stress
432        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
433            """
434            if eta !=None:
435                k=util.kronecker(self.domain.getDim())
436                kk=util.outer(k,k)
437                self.eta=util.interpolate(eta, escore.Function(self.domain))
438                self.__pde_prec.setValue(D=1/self.eta)
439                self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
440            if restoration_factor!=None:
441                n=self.domain.getNormal()
442                self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
443            if fixed_u_mask!=None:
444                self.__pde_v.setValue(q=fixed_u_mask)
445            if f!=None: self.__f=f
446            if surface_stress!=None: self.__surface_stress=surface_stress
447            if stress!=None: self.__stress=stress
448    
449         def initialize(self,f=escore.Data(),fixed_u_mask=escore.Data(),eta=1,surface_stress=escore.Data(),stress=escore.Data(), restoration_factor=0):
450            """
451            assigns values to the model parameters
452    
453            :param f: external force
454            :type f: `Vector` object in `FunctionSpace` `Function` or similar
455            :param fixed_u_mask: mask of locations with fixed velocity.
456            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
457            :param eta: viscosity
458            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
459            :param surface_stress: normal surface stress
460            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
461            :param stress: initial stress
462            :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
463            """
464            self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
465    
466         def Bv(self,v,tol):
467             """
468             returns inner product of element p and div(v)
469    
470             :param v: a residual
471             :return: inner product of element p and div(v)
472             :rtype: ``float``
473             """
474             self.__pde_proj.setValue(Y=-util.div(v))
475             self.getSolverOptionsDiv().setTolerance(tol)
476             self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
477             out=self.__pde_proj.getSolution()
478             return out
479    
480         def inner_pBv(self,p,Bv):
481             """
482             returns inner product of element p and Bv=-div(v)
483    
484             :param p: a pressure increment
485             :param Bv: a residual
486             :return: inner product of element p and Bv=-div(v)
487             :rtype: ``float``
488             """
489             return util.integrate(util.interpolate(p,escore.Function(self.domain))*util.interpolate(Bv, escore.Function(self.domain)))
490    
491         def inner_p(self,p0,p1):
492             """
493             Returns inner product of p0 and p1
494    
495             :param p0: a pressure
496             :param p1: a pressure
497             :return: inner product of p0 and p1
498             :rtype: ``float``
499             """
500             s0=util.interpolate(p0, escore.Function(self.domain))
501             s1=util.interpolate(p1, escore.Function(self.domain))
502           return util.integrate(s0*s1)           return util.integrate(s0*s1)
503    
504        def inner_a(self,a0,a1):       def norm_v(self,v):
505           p0=util.interpolate(a0[1],Function(self.domain))           """
506           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  
507    
508             :param v: a velovity
509             :return: norm of v
510             :rtype: non-negative ``float``
511             """
512             return util.sqrt(util.integrate(util.length(util.grad(v))**2))
513    
514    
515         def getDV(self, p, v, tol):
516             """
517             return the value for v for a given p
518    
519             :param p: a pressure
520             :param v: a initial guess for the value v to return.
521             :return: dv given as *Adv=(f-Av-B^*p)*
522             """
523             self.updateStokesEquation(v,p)
524             self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
525             self.getSolverOptionsVelocity().setTolerance(tol)
526             self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
527             if self.__stress.isEmpty():
528                self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
529             else:
530                self.__pde_v.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
531             out=self.__pde_v.getSolution()
532             return  out
533    
534         def norm_Bv(self,Bv):
535            """
536            Returns Bv (overwrite).
537    
538            :rtype: equal to the type of p
539            :note: boundary conditions on p should be zero!
540            """
541            return util.sqrt(util.integrate(util.interpolate(Bv, escore.Function(self.domain))**2))
542    
543         def solve_AinvBt(self,p, tol):
544             """
545             Solves *Av=B^*p* with accuracy `tol`
546    
547             :param p: a pressure increment
548             :return: the solution of *Av=B^*p*
549             :note: boundary conditions on v should be zero!
550             """
551             self.__pde_v.setValue(Y=escore.Data(), y=escore.Data(), X=-p*util.kronecker(self.domain))
552             out=self.__pde_v.getSolution()
553             return  out
554    
555         def solve_prec(self,Bv, tol):
556             """
557             applies preconditioner for for *BA^{-1}B^** to *Bv*
558             with accuracy ``self.getSubProblemTolerance()``
559    
560             :param Bv: velocity increment
561             :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
562             :note: boundary conditions on p are zero.
563             """
564             self.__pde_prec.setValue(Y=Bv)
565             self.getSolverOptionsPressure().setTolerance(tol)
566             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
567             out=self.__pde_prec.getSolution()
568             return out

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

  ViewVC Help
Powered by ViewVC 1.1.26