/[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 3515 by gross, Thu May 19 08:20:57 2011 UTC revision 3990 by caltinay, Tue Sep 25 05:03:20 2012 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-2012 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-2012 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 escript
37  import util  from . import util
38  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions  from .linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE, SolverOptions
39  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES  from .pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
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              
       self.__pde_v=LinearPDESystem(domain)  
       self.__pde_v.setSymmetryOn()  
       if self.useReduced: self.__pde_v.setReducedOrderOn()  
85        self.__pde_p=LinearSinglePDE(domain)        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=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 xrange(self.domain.getDim()) ] )
103      raise ValueError,"unknown solver %s"%self.solver           #self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale
104    
105          elif self.solver  == self.SMOOTH:
106             self.__pde_v=LinearPDESystem(domain)
107             self.__pde_v.setSymmetryOn()
108             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=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))        self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
113        self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))        self.__g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
114          self.__permeability_invXg=escript.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 = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))        self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
119        self.location_of_fixed_flux = escript.Vector(0, self.__pde_v.getFunctionSpaceForCoefficient("q"))        self.location_of_fixed_flux = escript.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
120        self.setTolerance()        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)
291              u= self.__pde_v.getSolution()  
292      elif self.solver  == self.POST:     def _getFlux(self, pp, u0=None):
293              self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g)-util.grad(p),          """
294                                    X=self.lamb * self.__f * util.kronecker(self.domain.getDim()))          returns the flux for a given pressure ``pp`` where the flux is equal to
295            ``u0`` on locations where ``location_of_fixed_flux`` is positive (see
296            `setValue`). Notice that ``g`` is used, see `setValue`.
297    
298            :param pp: pressure.
299            :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            :type u0: vector values on the domain (i.e. `escript.Data`) or ``None``
302            :return: flux
303            :rtype: `escript.Data`
304            """
305            if self.solver  == self.EVAL:
306               u = self.__g - util.tensor_mult(self.__permeability, self.perm_scale * (util.grad(pp) + self.__permeability_invXg_ref))
307            elif self.solver  == self.POST or self.solver  == self.SMOOTH:
308                self.__pde_v.setValue(Y= self.__permeability_invXg - (util.grad(pp) + self.__permeability_invXg_ref))
309                print
310              if u0 == None:              if u0 == None:
311             self.__pde_v.setValue(r=escript.Data())                 self.__pde_v.setValue(r=escript.Data())
312          else:              else:
313             self.__pde_v.setValue(r=u0)                 if not isinstance(u0, escript.Data) : u0 = escript.Vector(u0, escript.Solution(self.domain))
314              u= self.__pde_v.getSolution()                 self.__pde_v.setValue(r=1./self.perm_scale * u0)
315      elif self.solver  == self.STAB:              u= self.__pde_v.getSolution() * self.perm_scale
316           gp=util.grad(p)          return u
          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  
317                
       
    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  
   
    def __L2(self,v):  
          return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))        
     
    def __norm(self,r):  
          return util.sqrt(self.__inner(r,r))  
           
    def __inner(self,r,s):  
          return util.integrate(util.inner(r,s), escript.Function(self.domain))  
           
    def __STAB_Aprod(self,p):  
       gp=util.grad(p)  
       self.__pde_v.setValue(Y=-0.5*gp,  
                             X=-p*util.kronecker(self.__pde_v.getDomain()),  
                             y= p * self.domain.getNormal(),    
                             r=escript.Data())  
       u = -self.__pde_v.getSolution()  
       self.__pde_p.setValue(Y=util.div(u),  
                             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()  
         
   
318  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
319       """       """
320       solves       solves
# 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 459  class StokesProblemCartesian(Homogeneous Line 358  class StokesProblemCartesian(Homogeneous
358    
359           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=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 534  class StokesProblemCartesian(Homogeneous Line 433  class StokesProblemCartesian(Homogeneous
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, escript.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 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 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 654  class StokesProblemCartesian(Homogeneous Line 553  class StokesProblemCartesian(Homogeneous
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.3990

  ViewVC Help
Powered by ViewVC 1.1.26