/[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

trunk/escript/py_src/flows.py revision 3515 by gross, Thu May 19 08:20:57 2011 UTC trunk/escriptcore/py_src/flows.py revision 4507 by jfenwick, Wed Jul 24 02:50:22 2013 UTC
# Line 1  Line 1 
1  # -*- coding: utf-8 -*-  # -*- coding: utf-8 -*-
2  ########################################################  ##############################################################################
3  #  #
4  # Copyright (c) 2003-2010 by University of Queensland  # Copyright (c) 2003-2013 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # http://www.uq.edu.au
 # http://www.uq.edu.au/esscc  
6  #  #
7  # Primary Business: Queensland, Australia  # Primary Business: Queensland, Australia
8  # Licensed under the Open Software License version 3.0  # Licensed under the Open Software License version 3.0
9  # http://www.opensource.org/licenses/osl-3.0.php  # http://www.opensource.org/licenses/osl-3.0.php
10  #  #
11  ########################################################  # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    # Development since 2012 by School of Earth Sciences
13    #
14    ##############################################################################
15    
16  __copyright__="""Copyright (c) 2003-2010 by University of Queensland  __copyright__="""Copyright (c) 2003-2013 by University of Queensland
17  Earth Systems Science Computational Center (ESSCC)  http://www.uq.edu.au
 http://www.uq.edu.au/esscc  
18  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
19  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
20  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
# Line 32  Some models for flow Line 33  Some models for flow
33    
34  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
35    
36  import escript  from . import escriptcpp as escore
37  import util  from . import util
38  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions  from . import linearPDEs as lpe
39  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES  from . import pdetools as pdt
40    
41  class DarcyFlow(object):  class DarcyFlow(object):
42     """     """
# Line 46  class DarcyFlow(object): Line 47  class DarcyFlow(object):
47        
48     where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,     where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
49        
50     :cvar SIMPLE: simple solver     :cvar EVAL: direct pressure gradient evaluation for flux
51     :cvar POST: solver using global postprocessing of flux     :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*
52     :cvar STAB: solver uses (non-symmetric) stabilization                 where *l* is the length scale, *K* is the inverse of the permeability tensor, and *w* is a positive weighting factor.
53     :cvar SYMSTAB: solver uses symmetric stabilization     :cvar SMOOTH: global smoothing by solving the PDE *K_{ij} u_j= - p_{,j} + K_{ij} g_j*
54     """     """
55     SIMPLE="SIMPLE"     EVAL="EVAL"
56       SIMPLE="EVAL"
57     POST="POST"     POST="POST"
58     STAB="STAB"     SMOOTH="SMOOTH"
59     SYMSTAB="SYMSTAB"     def __init__(self, domain, useReduced=False, solver="POST", verbose=False, w=1.):
    def __init__(self, domain, useReduced=False, solver="SYMSTAB", verbose=False, w=1.):  
60        """        """
61        initializes the Darcy flux problem        initializes the Darcy flux problem.
62    
63        :param domain: domain of the problem        :param domain: domain of the problem
64        :type domain: `Domain`        :type domain: `Domain`
65        :param useReduced: uses reduced oreder on flux and pressure        :param useReduced: uses reduced oreder on flux and pressure
66        :type useReduced: ``bool``        :type useReduced: ``bool``
67        :param solver: solver method        :param solver: solver method
68        :type solver: in [`DarcyFlow.SIMPLE`, `DarcyFlow.POST', `DarcyFlow.STAB`, `DarcyFlow.SYMSTAB` ]        :type solver: in [`DarcyFlow.EVAL`, `DarcyFlow.POST`, `DarcyFlow.SMOOTH` ]
69        :param verbose: if ``True`` some information on the iteration progress are printed.        :param verbose: if ``True`` some information on the iteration progress are printed.
70        :type verbose: ``bool``        :type verbose: ``bool``
71        :param w: weighting factor for `DarcyFlow.POST` solver        :param w: weighting factor for `DarcyFlow.POST` solver
72        :type w: ``float``        :type w: ``float``
73                
74        """        """
75          if not solver in [DarcyFlow.EVAL, DarcyFlow.POST,  DarcyFlow.SMOOTH ] :
76              raise ValueError("unknown solver %d."%solver)
77    
78        self.domain=domain        self.domain=domain
79        self.solver=solver        self.solver=solver
80        self.useReduced=useReduced        self.useReduced=useReduced
81        self.verbose=verbose        self.verbose=verbose
82        self.scale=1.        self.l=None
83                self.w=None
84              
85        self.__pde_v=LinearPDESystem(domain)        self.__pde_p=lpe.LinearSinglePDE(domain)
       self.__pde_v.setSymmetryOn()  
       if self.useReduced: self.__pde_v.setReducedOrderOn()  
       self.__pde_p=LinearSinglePDE(domain)  
86        self.__pde_p.setSymmetryOn()        self.__pde_p.setSymmetryOn()
87        if self.useReduced: self.__pde_p.setReducedOrderOn()        if self.useReduced: self.__pde_p.setReducedOrderOn()
88          
89        if self.solver  == self.SIMPLE:        if self.solver  == self.EVAL:
90       if self.verbose: print "DarcyFlow: simple solver is used."           self.__pde_v=None
91           self.__pde_v.setValue(D=util.kronecker(self.domain.getDim()))           if self.verbose: print("DarcyFlow: simple solver is used.")
92    
93        elif self.solver  == self.POST:        elif self.solver  == self.POST:
94       self.w=w           if util.inf(w)<0.:
95       if util.inf(w)<0.:              raise ValueError("Weighting factor must be non-negative.")
96          raise ValueError,"Weighting factor must be non-negative."           if self.verbose: print("DarcyFlow: global postprocessing of flux is used.")
97       if self.verbose: print "DarcyFlow: global postprocessing of flux is used."           self.__pde_v=lpe.LinearPDESystem(domain)
98        elif self.solver  == self.STAB:           self.__pde_v.setSymmetryOn()
99        if self.verbose: print "DarcyFlow: (non-symmetric) stabilization is used."           if self.useReduced: self.__pde_v.setReducedOrderOn()
100        elif  self.solver  == self.SYMSTAB:           self.w=w
101        if self.verbose: print "DarcyFlow: symmetric stabilization is used."           x=self.domain.getX()
102        else:           self.l=min( [util.sup(x[i])-util.inf(x[i]) for i in range(self.domain.getDim()) ] )
103      raise ValueError,"unknown solver %s"%self.solver           #self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale
104        self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))  
105        self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))        elif self.solver  == self.SMOOTH:
106        self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))           self.__pde_v=lpe.LinearPDESystem(domain)
107        self.location_of_fixed_flux = escript.Vector(0, self.__pde_v.getFunctionSpaceForCoefficient("q"))           self.__pde_v.setSymmetryOn()
108        self.setTolerance()           if self.useReduced: self.__pde_v.setReducedOrderOn()
109             if self.verbose: print("DarcyFlow: flux smoothing is used.")
110             self.w=0
111    
112          self.__f=escore.Data(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
113          self.__g=escore.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
114          self.__permeability_invXg=escore.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
115          self.__permeability_invXg_ref=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)
116          self.ref_point_id=None
117          self.ref_point=util.numpy.zeros((self.domain.getDim()),util.numpy.float64)
118          self.location_of_fixed_pressure = escore.Data(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
119          self.location_of_fixed_flux = escore.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
120          self.perm_scale=1.
121            
122                    
123     def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):     def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
# Line 128  class DarcyFlow(object): Line 143  class DarcyFlow(object):
143    
144        """        """
145        if location_of_fixed_pressure!=None:        if location_of_fixed_pressure!=None:
146             self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)             self.location_of_fixed_pressure=util.wherePositive(util.interpolate(location_of_fixed_pressure, self.__pde_p.getFunctionSpaceForCoefficient("q")))
147               self.ref_point_id=self.location_of_fixed_pressure.maxGlobalDataPoint()
148               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.")
149               self.ref_point=self.__pde_p.getFunctionSpaceForCoefficient("q").getX().getTupleForGlobalDataPoint(*self.ref_point_id)
150               if self.verbose: print(("DarcyFlow: reference point at %s."%(self.ref_point,)))
151             self.__pde_p.setValue(q=self.location_of_fixed_pressure)             self.__pde_p.setValue(q=self.location_of_fixed_pressure)
152        if location_of_fixed_flux!=None:        if location_of_fixed_flux!=None:
153            self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)            self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
154            self.__pde_v.setValue(q=self.location_of_fixed_flux)            if not self.__pde_v == None:
155                        self.__pde_v.setValue(q=self.location_of_fixed_flux)
156                            
       # pressure  is rescaled by the factor 1/self.scale  
157        if permeability!=None:        if permeability!=None:
158            
159       perm=util.interpolate(permeability,self.__pde_v.getFunctionSpaceForCoefficient("A"))           perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
160           V=util.vol(self.domain)           self.perm_scale=util.Lsup(util.length(perm))
161           l=V**(1./self.domain.getDim())           if self.verbose: print(("DarcyFlow: permeability scaling factor = %e."%self.perm_scale))
162             perm=perm*(1./self.perm_scale)
163                    
164       if perm.getRank()==0:           if perm.getRank()==0:
165          perm_inv=(1./perm)  
166              self.scale=util.integrate(perm_inv)/V*l              perm_inv=(1./perm)
167          perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))              perm_inv=perm_inv*util.kronecker(self.domain.getDim())
168          perm=perm*(self.scale*util.kronecker(self.domain.getDim()))              perm=perm*util.kronecker(self.domain.getDim())
169                    
170                    
171       elif perm.getRank()==2:           elif perm.getRank()==2:
172          perm_inv=util.inverse(perm)              perm_inv=util.inverse(perm)
173              self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)/V)*l           else:
174          perm_inv*=(1./self.scale)              raise ValueError("illegal rank of permeability.")
         perm=perm*self.scale  
      else:  
         raise ValueError,"illegal rank of permeability."  
175                    
176       self.__permeability=perm           self.__permeability=perm
177       self.__permeability_inv=perm_inv           self.__permeability_inv=perm_inv
      if self.verbose: print "DarcyFlow: scaling factor for pressure is %e."%self.scale  
178            
179       if self.solver  == self.SIMPLE:           #====================
180          self.__pde_p.setValue(A=self.__permeability)           self.__pde_p.setValue(A=self.__permeability)
181       elif self.solver  == self.POST:           if self.solver  == self.EVAL:
182          self.__pde_p.setValue(A=self.__permeability)                pass # no extra work required
183          k=util.kronecker(self.domain.getDim())           elif self.solver  == self.POST:
184          self.lamb = self.w*util.length(perm_inv)*l                k=util.kronecker(self.domain.getDim())
185          self.__pde_v.setValue(D=self.__permeability_inv, A=self.lamb*self.domain.getSize()*util.outer(k,k))                self.omega = self.w*util.length(perm_inv)*self.l*self.domain.getSize()
186       elif self.solver  == self.STAB:                #self.__pde_v.setValue(D=self.__permeability_inv, A=self.omega*util.outer(k,k))
187          self.__pde_p.setValue(A=0.5*self.__permeability)                self.__pde_v.setValue(D=self.__permeability_inv, A_reduced=self.omega*util.outer(k,k))
188          self.__pde_v.setValue(D=0.5*self.__permeability_inv)           elif self.solver  == self.SMOOTH:
189       elif  self.solver  == self.SYMSTAB:              self.__pde_v.setValue(D=self.__permeability_inv)
         self.__pde_p.setValue(A=0.5*self.__permeability)  
         self.__pde_v.setValue(D=0.5*self.__permeability_inv)  
190    
191        if g != None:        if g != None:
192      g=util.interpolate(g, self.__pde_v.getFunctionSpaceForCoefficient("Y"))          g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
193      if g.isEmpty():          if g.isEmpty():
194            g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))               g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
195      else:          else:
196          if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g"               if not g.getShape()==(self.domain.getDim(),): raise ValueError("illegal shape of g")
197      self.__g=g          self.__g=g
198            self.__permeability_invXg=util.tensor_mult(self.__permeability_inv,self.__g * (1./self.perm_scale ))
199            self.__permeability_invXg_ref=util.integrate(self.__permeability_invXg)/util.vol(self.domain)
200        if f !=None:        if f !=None:
201       f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))           f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
202       if f.isEmpty():                 if f.isEmpty():      
203            f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))               f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
204       else:           else:
205           if f.getRank()>0: raise ValueError,"illegal rank of f."               if f.getRank()>0: raise ValueError("illegal rank of f.")
206       self.__f=f           self.__f=f
207    
208     def getSolverOptionsFlux(self):     def getSolverOptionsFlux(self):
209        """        """
210        Returns the solver options used to solve the flux problems        Returns the solver options used to solve the flux problems
211        :return: `SolverOptions`        :return: `SolverOptions`
212        """        """
213        return self.__pde_v.getSolverOptions()        if self.__pde_v == None:
214              return None
215          else:
216              return self.__pde_v.getSolverOptions()
217                
218     def setSolverOptionsFlux(self, options=None):     def setSolverOptionsFlux(self, options=None):
219        """        """
# Line 202  class DarcyFlow(object): Line 221  class DarcyFlow(object):
221        If ``options`` is not present, the options are reset to default        If ``options`` is not present, the options are reset to default
222        :param options: `SolverOptions`        :param options: `SolverOptions`
223        """        """
224        return self.__pde_v.setSolverOptions(options)        if not self.__pde_v == None:
225              self.__pde_v.setSolverOptions(options)
226            
227     def getSolverOptionsPressure(self):     def getSolverOptionsPressure(self):
228        """        """
# Line 221  class DarcyFlow(object): Line 241  class DarcyFlow(object):
241        """        """
242        return self.__pde_p.setSolverOptions(options)        return self.__pde_p.setSolverOptions(options)
243                
244     def setTolerance(self,rtol=1e-4):     def solve(self, u0, p0):
       """  
       sets the relative tolerance ``rtol`` for the pressure for the stabelized solvers.  
         
       :param rtol: relative tolerance for the pressure  
       :type rtol: non-negative ``float``  
       """  
       if rtol<0:  
      raise ValueError,"Relative tolerance needs to be non-negative."  
       self.__rtol=rtol  
         
    def getTolerance(self):  
       """  
       returns the relative tolerance  
       :return: current relative tolerance  
       :rtype: ``float``  
       """  
       return self.__rtol  
         
    def solve(self,u0,p0, max_iter=100, iter_restart=20):  
245        """        """
246        solves the problem.        solves the problem.
247                
       The iteration is terminated if the residual norm is less then self.getTolerance().  
   
248        :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.        :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.
249        :type u0: vector value on the domain (e.g. `escript.Data`).        :type u0: vector value on the domain (e.g. `escript.Data`).
250        :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.        :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.
251        :type p0: scalar value on the domain (e.g. `escript.Data`).        :type p0: scalar value on the domain (e.g. `escript.Data`).
       :param max_iter: maximum number of (outer) iteration steps for the stabilization solvers,  
       :type max_iter: ``int``  
       :param iter_restart: number of steps after which the iteration is restarted. The larger ``iter_restart`` the larger the required memory.  
                            A small value for ``iter_restart`` may require a large number of iteration steps or may even lead to a failure  
                            of the iteration. ``iter_restart`` is relevant for the stabilization solvers only.  
       :type iter_restart: ``int``  
252        :return: flux and pressure        :return: flux and pressure
253        :rtype: ``tuple`` of `escript.Data`.        :rtype: ``tuple`` of `escript.Data`.
254    
255        """        """
256        # rescale initial guess:        p0=util.interpolate(p0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
257        p0=p0/self.scale        if self.ref_point_id == None:
258        if self.solver  == self.SIMPLE or self.solver  == self.POST :            p_ref=0
259          self.__pde_p.setValue(X=self.__g ,        else:
260                                Y=self.__f,            p_ref=p0.getTupleForGlobalDataPoint(*self.ref_point_id)[0]
261                                y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux),        p0_hydrostatic=p_ref+util.inner(self.__permeability_invXg_ref, self.__pde_p.getFunctionSpaceForCoefficient("q").getX() - self.ref_point)
262                                r=p0)        g_2=self.__g - util.tensor_mult(self.__permeability, self.__permeability_invXg_ref * self.perm_scale)
263          p=self.__pde_p.getSolution()        self.__pde_p.setValue(X=g_2 * 1./self.perm_scale,
264          u = self.getFlux(p, u0)                              Y=self.__f * 1./self.perm_scale,
265        elif  self.solver  == self.STAB:                              y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux * 1./self.perm_scale ),
266      u,p = self.__solve_STAB(u0,p0, max_iter, iter_restart)                              r=p0 - p0_hydrostatic)
267        elif  self.solver  == self.SYMSTAB:        pp=self.__pde_p.getSolution()
268      u,p = self.__solve_SYMSTAB(u0,p0, max_iter, iter_restart)        u = self._getFlux(pp, u0)
269              return u,pp + p0_hydrostatic
       if self.verbose:  
         KGp=util.tensor_mult(self.__permeability,util.grad(p))  
         def_p=self.__g-(u+KGp)  
         def_v=self.__f-util.div(u, self.__pde_v.getFunctionSpaceForCoefficient("X"))  
         print "DarcyFlux: |g-u-K*grad(p)|_2 = %e (|u|_2 = %e)."%(self.__L2(def_p),self.__L2(u))  
         print "DarcyFlux: |f-div(u)|_2 = %e (|grad(u)|_2 = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))  
       #rescale result  
       p=p*self.scale  
       return u,p  
270                
271     def getFlux(self,p, u0=None):     def getFlux(self,p, u0=None):
272          """          """
273          returns the flux for a given pressure ``p`` where the flux is equal to ``u0``          returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
274          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
275          Notice that ``g`` and ``f`` are used, see `setValue`.          Notice that ``g`` is used, see `setValue`.
276    
277          :param p: pressure.          :param p: pressure.
278          :type p: scalar value on the domain (e.g. `escript.Data`).          :type p: scalar value on the domain (e.g. `escript.Data`).
# Line 297  class DarcyFlow(object): Line 281  class DarcyFlow(object):
281          :return: flux          :return: flux
282          :rtype: `escript.Data`          :rtype: `escript.Data`
283          """          """
284          if self.solver  == self.SIMPLE or self.solver  == self.POST  :          p=util.interpolate(p, self.__pde_p.getFunctionSpaceForCoefficient("q"))
285              KGp=util.tensor_mult(self.__permeability,util.grad(p))          if self.ref_point_id == None:
286              self.__pde_v.setValue(Y=self.__g-KGp, X=escript.Data())              p_ref=0
287              if u0 == None:          else:
288             self.__pde_v.setValue(r=escript.Data())              p_ref=p.getTupleForGlobalDataPoint(*self.ref_point_id)[0]
289          else:          p_hydrostatic=p_ref+util.inner(self.__permeability_invXg_ref, self.__pde_p.getFunctionSpaceForCoefficient("q").getX() - self.ref_point)
290             self.__pde_v.setValue(r=u0)          return self._getFlux(p-p_hydrostatic, u0)
             u= self.__pde_v.getSolution()  
     elif self.solver  == self.POST:  
             self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g)-util.grad(p),  
                                   X=self.lamb * self.__f * util.kronecker(self.domain.getDim()))  
             if u0 == None:  
            self.__pde_v.setValue(r=escript.Data())  
         else:  
            self.__pde_v.setValue(r=u0)  
             u= self.__pde_v.getSolution()  
     elif self.solver  == self.STAB:  
          gp=util.grad(p)  
          self.__pde_v.setValue(Y=0.5*(util.tensor_mult(self.__permeability_inv,self.__g)+gp),  
                                X= p * util.kronecker(self.domain.getDim()),  
                                y= - p * self.domain.getNormal())                            
          if u0 == None:  
            self.__pde_v.setValue(r=escript.Data())  
          else:  
            self.__pde_v.setValue(r=u0)  
          u= self.__pde_v.getSolution()  
     elif  self.solver  == self.SYMSTAB:  
          gp=util.grad(p)  
          self.__pde_v.setValue(Y=0.5*(util.tensor_mult(self.__permeability_inv,self.__g)-gp),  
                                X= escript.Data() ,  
                                y= escript.Data() )                            
          if u0 == None:  
            self.__pde_v.setValue(r=escript.Data())  
          else:  
            self.__pde_v.setValue(r=u0)  
          u= self.__pde_v.getSolution()  
     return u  
         
       
    def __solve_STAB(self, u0, p0, max_iter, iter_restart):  
           # p0 is used as an initial guess  
       u=self.getFlux(p0, u0)    
           self.__pde_p.setValue( Y=self.__f-util.div(u),  
                                  X=0.5*(self.__g - u - util.tensor_mult(self.__permeability,util.grad(p0)) ),  
                                  y= escript.Data(),  
                                  r=escript.Data())  
   
       dp=self.__pde_p.getSolution()  
       p=GMRES(dp,  
               self.__STAB_Aprod,  
           p0,  
           self.__inner,  
           atol=self.__norm(p0+dp)*self.getTolerance() ,  
           rtol=0.,  
           iter_max=max_iter,  
           iter_restart=iter_restart,  
           verbose=self.verbose,P_R=None)  
               
           u=self.getFlux(p, u0)  
           return u,p  
   
    def __solve_SYMSTAB(self, u0, p0, max_iter, iter_restart):  
           # p0 is used as an initial guess  
       u=self.getFlux(p0, u0)    
           self.__pde_p.setValue( Y= self.__f,  
                                  X=  0.5*(self.__g + u - util.tensor_mult(self.__permeability,util.grad(p0)) ),  
                                  y=  -  util.inner(self.domain.getNormal(), u),  
                                  r=escript.Data())  
       dp=self.__pde_p.getSolution()  
         
       p=GMRES(dp,  
               self.__SYMSTAB_Aprod,  
           p0,  
           self.__inner,  
           atol=self.__norm(p0+dp)*self.getTolerance() ,  
           rtol=0.,  
           iter_max=max_iter,  
           iter_restart=iter_restart,  
           verbose=self.verbose,P_R=None)  
               
           u=self.getFlux(p, u0)  
           return u,p  
291    
292     def __L2(self,v):     def _getFlux(self, pp, u0=None):
293           return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))                """
294              returns the flux for a given pressure ``pp`` where the flux is equal to
295     def __norm(self,r):          ``u0`` on locations where ``location_of_fixed_flux`` is positive (see
296           return util.sqrt(self.__inner(r,r))          `setValue`). Notice that ``g`` is used, see `setValue`.
297            
298     def __inner(self,r,s):          :param pp: pressure.
299           return util.integrate(util.inner(r,s), escript.Function(self.domain))          :type pp: scalar value on the domain (i.e. `escript.Data`).
300                    :param u0: flux on the locations of the domain marked in ``location_of_fixed_flux``.
301     def __STAB_Aprod(self,p):          :type u0: vector values on the domain (i.e. `escript.Data`) or ``None``
302        gp=util.grad(p)          :return: flux
303        self.__pde_v.setValue(Y=-0.5*gp,          :rtype: `escript.Data`
304                              X=-p*util.kronecker(self.__pde_v.getDomain()),          """
305                              y= p * self.domain.getNormal(),            if self.solver  == self.EVAL:
306                              r=escript.Data())             u = self.__g - util.tensor_mult(self.__permeability, self.perm_scale * (util.grad(pp) + self.__permeability_invXg_ref))
307        u = -self.__pde_v.getSolution()          elif self.solver  == self.POST or self.solver  == self.SMOOTH:
308        self.__pde_p.setValue(Y=util.div(u),              self.__pde_v.setValue(Y= self.__permeability_invXg - (util.grad(pp) + self.__permeability_invXg_ref))
                             X=0.5*(u+util.tensor_mult(self.__permeability,gp)),  
                             y=escript.Data(),  
                             r=escript.Data())  
       
       return  self.__pde_p.getSolution()  
     
    def __SYMSTAB_Aprod(self,p):  
       gp=util.grad(p)  
       self.__pde_v.setValue(Y=0.5*gp ,  
                             X=escript.Data(),  
                             y=escript.Data(),    
                             r=escript.Data())  
       u = -self.__pde_v.getSolution()  
       self.__pde_p.setValue(Y=escript.Data(),  
                             X=0.5*(-u+util.tensor_mult(self.__permeability,gp)),  
                             y=escript.Data(),  
                             r=escript.Data())  
       
       return  self.__pde_p.getSolution()  
         
309    
310  class StokesProblemCartesian(HomogeneousSaddlePointProblem):              if u0 == None:
311                   self.__pde_v.setValue(r=escore.Data())
312                else:
313                   if not isinstance(u0, escore.Data) : u0 = escore.Vector(u0, escore.Solution(self.domain))
314                   self.__pde_v.setValue(r=1./self.perm_scale * u0)
315                u= self.__pde_v.getSolution() * self.perm_scale
316            return u
317          
318    class StokesProblemCartesian(pdt.HomogeneousSaddlePointProblem):
319       """       """
320       solves       solves
321    
# Line 436  class StokesProblemCartesian(Homogeneous Line 333  class StokesProblemCartesian(Homogeneous
333              sp.setTolerance()              sp.setTolerance()
334              sp.initialize(...)              sp.initialize(...)
335              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
336                sp.setStokesEquation(...) # new values for some parameters
337                v1,p1=sp.solve(v,p)
338       """       """
339       def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
340           """           """
# Line 448  class StokesProblemCartesian(Homogeneous Line 347  class StokesProblemCartesian(Homogeneous
347           :param domain: domain of the problem.           :param domain: domain of the problem.
348           :type domain: `Domain`           :type domain: `Domain`
349           """           """
350           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           pdt.HomogeneousSaddlePointProblem.__init__(self,**kwargs)
351           self.domain=domain           self.domain=domain
352           self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_v=lpe.LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
353           self.__pde_v.setSymmetryOn()           self.__pde_v.setSymmetryOn()
354            
355           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=lpe.LinearPDE(domain)
356           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
357           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
358    
359           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=lpe.LinearPDE(domain)
360           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
361       self.__pde_proj.setValue(D=1)           self.__pde_proj.setValue(D=1)
362           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
363    
364       def getSolverOptionsVelocity(self):       def getSolverOptionsVelocity(self):
# Line 468  class StokesProblemCartesian(Homogeneous Line 367  class StokesProblemCartesian(Homogeneous
367            
368       :rtype: `SolverOptions`       :rtype: `SolverOptions`
369       """       """
370       return self.__pde_v.getSolverOptions()           return self.__pde_v.getSolverOptions()
371       def setSolverOptionsVelocity(self, options=None):       def setSolverOptionsVelocity(self, options=None):
372           """           """
373       set the solver options for solving the equation for velocity.       set the solver options for solving the equation for velocity.
# Line 482  class StokesProblemCartesian(Homogeneous Line 381  class StokesProblemCartesian(Homogeneous
381       returns the solver options used  solve the equation for pressure.       returns the solver options used  solve the equation for pressure.
382       :rtype: `SolverOptions`       :rtype: `SolverOptions`
383       """       """
384       return self.__pde_prec.getSolverOptions()           return self.__pde_prec.getSolverOptions()
385       def setSolverOptionsPressure(self, options=None):       def setSolverOptionsPressure(self, options=None):
386           """           """
387       set the solver options for solving the equation for pressure.       set the solver options for solving the equation for pressure.
388       :param options: new solver  options       :param options: new solver  options
389       :type options: `SolverOptions`       :type options: `SolverOptions`
390       """       """
391       self.__pde_prec.setSolverOptions(options)           self.__pde_prec.setSolverOptions(options)
392    
393       def setSolverOptionsDiv(self, options=None):       def setSolverOptionsDiv(self, options=None):
394           """           """
# Line 499  class StokesProblemCartesian(Homogeneous Line 398  class StokesProblemCartesian(Homogeneous
398       :param options: new solver options       :param options: new solver options
399       :type options: `SolverOptions`       :type options: `SolverOptions`
400       """       """
401       self.__pde_proj.setSolverOptions(options)           self.__pde_proj.setSolverOptions(options)
402       def getSolverOptionsDiv(self):       def getSolverOptionsDiv(self):
403           """           """
404       returns the solver options for solving the equation to project the divergence of       returns the solver options for solving the equation to project the divergence of
# Line 507  class StokesProblemCartesian(Homogeneous Line 406  class StokesProblemCartesian(Homogeneous
406            
407       :rtype: `SolverOptions`       :rtype: `SolverOptions`
408       """       """
409       return self.__pde_proj.getSolverOptions()           return self.__pde_proj.getSolverOptions()
410    
411       def updateStokesEquation(self, v, p):       def updateStokesEquation(self, v, p):
412           """           """
413           updates the Stokes equation to consider dependencies from ``v`` and ``p``           updates the Stokes equation to consider dependencies from ``v`` and ``p``
414           :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.           :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values to model parameters.
415           """           """
416           pass           pass
417       def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):       def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
# Line 533  class StokesProblemCartesian(Homogeneous Line 432  class StokesProblemCartesian(Homogeneous
432          if eta !=None:          if eta !=None:
433              k=util.kronecker(self.domain.getDim())              k=util.kronecker(self.domain.getDim())
434              kk=util.outer(k,k)              kk=util.outer(k,k)
435              self.eta=util.interpolate(eta, escript.Function(self.domain))              self.eta=util.interpolate(eta, escore.Function(self.domain))
436          self.__pde_prec.setValue(D=1/self.eta)              self.__pde_prec.setValue(D=1/self.eta)
437              self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))              self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
438          if restoration_factor!=None:          if restoration_factor!=None:
439              n=self.domain.getNormal()              n=self.domain.getNormal()
# Line 545  class StokesProblemCartesian(Homogeneous Line 444  class StokesProblemCartesian(Homogeneous
444          if surface_stress!=None: self.__surface_stress=surface_stress          if surface_stress!=None: self.__surface_stress=surface_stress
445          if stress!=None: self.__stress=stress          if stress!=None: self.__stress=stress
446    
447       def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):       def initialize(self,f=escore.Data(),fixed_u_mask=escore.Data(),eta=1,surface_stress=escore.Data(),stress=escore.Data(), restoration_factor=0):
448          """          """
449          assigns values to the model parameters          assigns values to the model parameters
450    
# Line 558  class StokesProblemCartesian(Homogeneous Line 457  class StokesProblemCartesian(Homogeneous
457          :param surface_stress: normal surface stress          :param surface_stress: normal surface stress
458          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
459          :param stress: initial stress          :param stress: initial stress
460      :type stress: `Tensor` object on `FunctionSpace` `Function` or similar          :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
461          """          """
462          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
463    
# Line 571  class StokesProblemCartesian(Homogeneous Line 470  class StokesProblemCartesian(Homogeneous
470           :rtype: ``float``           :rtype: ``float``
471           """           """
472           self.__pde_proj.setValue(Y=-util.div(v))           self.__pde_proj.setValue(Y=-util.div(v))
473       self.getSolverOptionsDiv().setTolerance(tol)           self.getSolverOptionsDiv().setTolerance(tol)
474       self.getSolverOptionsDiv().setAbsoluteTolerance(0.)           self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
475           out=self.__pde_proj.getSolution()           out=self.__pde_proj.getSolution()
476           return out           return out
477    
# Line 585  class StokesProblemCartesian(Homogeneous Line 484  class StokesProblemCartesian(Homogeneous
484           :return: inner product of element p and Bv=-div(v)           :return: inner product of element p and Bv=-div(v)
485           :rtype: ``float``           :rtype: ``float``
486           """           """
487           return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))           return util.integrate(util.interpolate(p,escore.Function(self.domain))*util.interpolate(Bv, escore.Function(self.domain)))
488    
489       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
490           """           """
# Line 596  class StokesProblemCartesian(Homogeneous Line 495  class StokesProblemCartesian(Homogeneous
495           :return: inner product of p0 and p1           :return: inner product of p0 and p1
496           :rtype: ``float``           :rtype: ``float``
497           """           """
498           s0=util.interpolate(p0, escript.Function(self.domain))           s0=util.interpolate(p0, escore.Function(self.domain))
499           s1=util.interpolate(p1, escript.Function(self.domain))           s1=util.interpolate(p1, escore.Function(self.domain))
500           return util.integrate(s0*s1)           return util.integrate(s0*s1)
501    
502       def norm_v(self,v):       def norm_v(self,v):
# Line 613  class StokesProblemCartesian(Homogeneous Line 512  class StokesProblemCartesian(Homogeneous
512    
513       def getDV(self, p, v, tol):       def getDV(self, p, v, tol):
514           """           """
515           return the value for v for a given p (overwrite)           return the value for v for a given p
516    
517           :param p: a pressure           :param p: a pressure
518           :param v: a initial guess for the value v to return.           :param v: a initial guess for the value v to return.
# Line 621  class StokesProblemCartesian(Homogeneous Line 520  class StokesProblemCartesian(Homogeneous
520           """           """
521           self.updateStokesEquation(v,p)           self.updateStokesEquation(v,p)
522           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
523       self.getSolverOptionsVelocity().setTolerance(tol)           self.getSolverOptionsVelocity().setTolerance(tol)
524       self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)           self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
525           if self.__stress.isEmpty():           if self.__stress.isEmpty():
526              self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))              self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
527           else:           else:
# Line 637  class StokesProblemCartesian(Homogeneous Line 536  class StokesProblemCartesian(Homogeneous
536          :rtype: equal to the type of p          :rtype: equal to the type of p
537          :note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
538          """          """
539          return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))          return util.sqrt(util.integrate(util.interpolate(Bv, escore.Function(self.domain))**2))
540    
541       def solve_AinvBt(self,p, tol):       def solve_AinvBt(self,p, tol):
542           """           """
# Line 647  class StokesProblemCartesian(Homogeneous Line 546  class StokesProblemCartesian(Homogeneous
546           :return: the solution of *Av=B^*p*           :return: the solution of *Av=B^*p*
547           :note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
548           """           """
549           self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))           self.__pde_v.setValue(Y=escore.Data(), y=escore.Data(), X=-p*util.kronecker(self.domain))
550           out=self.__pde_v.getSolution()           out=self.__pde_v.getSolution()
551           return  out           return  out
552    
553       def solve_prec(self,Bv, tol):       def solve_prec(self,Bv, tol):
554           """           """
555           applies preconditioner for for *BA^{-1}B^** to *Bv*           applies preconditioner for for *BA^{-1}B^** to *Bv*
556           with accuracy `self.getSubProblemTolerance()`           with accuracy ``self.getSubProblemTolerance()``
557    
558           :param Bv: velocity increment           :param Bv: velocity increment
559           :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*           :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
560           :note: boundary conditions on p are zero.           :note: boundary conditions on p are zero.
561           """           """
562           self.__pde_prec.setValue(Y=Bv)           self.__pde_prec.setValue(Y=Bv)
563       self.getSolverOptionsPressure().setTolerance(tol)           self.getSolverOptionsPressure().setTolerance(tol)
564       self.getSolverOptionsPressure().setAbsoluteTolerance(0.)           self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
565           out=self.__pde_prec.getSolution()           out=self.__pde_prec.getSolution()
566           return out           return out

Legend:
Removed from v.3515  
changed lines
  Added in v.4507

  ViewVC Help
Powered by ViewVC 1.1.26