/[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 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC revision 3990 by caltinay, Tue Sep 25 05:03:20 2012 UTC
# Line 1  Line 1 
1  ########################################################  # -*- coding: utf-8 -*-
2    ##############################################################################
3  #  #
4  # Copyright (c) 2003-2009 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-2008 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 21  __url__="https://launchpad.net/escript-f Line 23  __url__="https://launchpad.net/escript-f
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"
35    
36  from escript import *  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      """     """
43      solves the problem     solves the problem
44      
45      M{u_i+k_{ij}*p_{,j} = g_i}     *u_i+k_{ij}*p_{,j} = g_i*
46      M{u_{i,i} = f}     *u_{i,i} = f*
47      
48      where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,     where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
49      
50      @note: The problem is solved in a least squares formulation.     :cvar EVAL: direct pressure gradient evaluation for flux
51      """     :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                   where *l* is the length scale, *K* is the inverse of the permeability tensor, and *w* is a positive weighting factor.
53      def __init__(self, domain, weight=None, useReduced=False, adaptSubTolerance=True):     :cvar SMOOTH: global smoothing by solving the PDE *K_{ij} u_j= - p_{,j} + K_{ij} g_j*
54          """     """
55          initializes the Darcy flux problem     EVAL="EVAL"
56          @param domain: domain of the problem     SIMPLE="EVAL"
57          @type domain: L{Domain}     POST="POST"
58      @param useReduced: uses reduced oreder on flux and pressure     SMOOTH="SMOOTH"
59      @type useReduced: C{bool}     def __init__(self, domain, useReduced=False, solver="POST", verbose=False, w=1.):
60      @param adaptSubTolerance: switches on automatic subtolerance selection        """
61      @type adaptSubTolerance: C{bool}            initializes the Darcy flux problem.
62          """  
63          self.domain=domain        :param domain: domain of the problem
64          if weight == None:        :type domain: `Domain`
65             s=self.domain.getSize()        :param useReduced: uses reduced oreder on flux and pressure
66             self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2        :type useReduced: ``bool``
67             # self.__l=(3.*util.longestEdge(self.domain))**2        :param solver: solver method
68             # self.__l=(0.1*util.longestEdge(self.domain)*s/util.sup(s))**2        :type solver: in [`DarcyFlow.EVAL`, `DarcyFlow.POST`, `DarcyFlow.SMOOTH` ]
69          else:        :param verbose: if ``True`` some information on the iteration progress are printed.
70             self.__l=weight        :type verbose: ``bool``
71          self.__pde_v=LinearPDESystem(domain)        :param w: weighting factor for `DarcyFlow.POST` solver
72          if useReduced: self.__pde_v.setReducedOrderOn()        :type w: ``float``
73          self.__pde_v.setSymmetryOn()        
74          self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))        """
75          self.__pde_p=LinearSinglePDE(domain)        if not solver in [DarcyFlow.EVAL, DarcyFlow.POST,  DarcyFlow.SMOOTH ] :
76          self.__pde_p.setSymmetryOn()            raise ValueError("unknown solver %d."%solver)
77          if useReduced: self.__pde_p.setReducedOrderOn()  
78          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))        self.domain=domain
79          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))        self.solver=solver
80          self.setTolerance()        self.useReduced=useReduced
81          self.setAbsoluteTolerance()        self.verbose=verbose
82      self.__adaptSubTolerance=adaptSubTolerance        self.l=None
83      self.verbose=False        self.w=None
84      def getSolverOptionsFlux(self):      
85      """        self.__pde_p=LinearSinglePDE(domain)
86      Returns the solver options used to solve the flux problems        self.__pde_p.setSymmetryOn()
87              if self.useReduced: self.__pde_p.setReducedOrderOn()
88      M{(I+D^*D)u=F}  
89              if self.solver  == self.EVAL:
90      @return: L{SolverOptions}           self.__pde_v=None
91      """           if self.verbose: print("DarcyFlow: simple solver is used.")
92      return self.__pde_v.getSolverOptions()  
93      def setSolverOptionsFlux(self, options=None):        elif self.solver  == self.POST:
94      """           if util.inf(w)<0.:
95      Sets the solver options used to solve the flux problems              raise ValueError("Weighting factor must be non-negative.")
96                 if self.verbose: print("DarcyFlow: global postprocessing of flux is used.")
97      M{(I+D^*D)u=F}           self.__pde_v=LinearPDESystem(domain)
98                 self.__pde_v.setSymmetryOn()
99      If C{options} is not present, the options are reset to default           if self.useReduced: self.__pde_v.setReducedOrderOn()
100      @param options: L{SolverOptions}           self.w=w
101      @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.           x=self.domain.getX()
102      """           self.l=min( [util.sup(x[i])-util.inf(x[i]) for i in xrange(self.domain.getDim()) ] )
103      return self.__pde_v.setSolverOptions(options)           #self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale
104      def getSolverOptionsPressure(self):  
105      """        elif self.solver  == self.SMOOTH:
106      Returns the solver options used to solve the pressure problems           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"))
113          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"))
119          self.location_of_fixed_flux = escript.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):
124          """
125          assigns values to model parameters
126    
127          :param f: volumetic sources/sinks
128          :type f: scalar value on the domain (e.g. `escript.Data`)
129          :param g: flux sources/sinks
130          :type g: vector values on the domain (e.g. `escript.Data`)
131          :param location_of_fixed_pressure: mask for locations where pressure is fixed
132          :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
133          :param location_of_fixed_flux:  mask for locations where flux is fixed.
134          :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
135          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
136          :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
137    
138          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
139          :note: at any point on the boundary of the domain the pressure
140                 (``location_of_fixed_pressure`` >0) or the normal component of the
141                 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
142                 is along the *x_i* axis.
143    
144          """
145          if location_of_fixed_pressure!=None:
146               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)
152          if location_of_fixed_flux!=None:
153              self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
154              if not self.__pde_v == None:
155                  self.__pde_v.setValue(q=self.location_of_fixed_flux)
156                
157          if permeability!=None:
158            
159      M{(Q^*Q)p=Q^*G}           perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
160                 self.perm_scale=util.Lsup(util.length(perm))
161      @return: L{SolverOptions}           if self.verbose: print(("DarcyFlow: permeability scaling factor = %e."%self.perm_scale))
162      """           perm=perm*(1./self.perm_scale)
163      return self.__pde_p.getSolverOptions()          
164      def setSolverOptionsPressure(self, options=None):           if perm.getRank()==0:
165      """  
166      Sets the solver options used to solve the pressure problems              perm_inv=(1./perm)
167                    perm_inv=perm_inv*util.kronecker(self.domain.getDim())
168      M{(Q^*Q)p=Q^*G}              perm=perm*util.kronecker(self.domain.getDim())
169                
170      If C{options} is not present, the options are reset to default          
171      @param options: L{SolverOptions}           elif perm.getRank()==2:
172      @note: if the adaption of subtolerance is choosen, the tolerance set by C{options} will be overwritten before the solver is called.              perm_inv=util.inverse(perm)
173      """           else:
174      return self.__pde_p.setSolverOptions(options)              raise ValueError("illegal rank of permeability.")
175            
176      def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):           self.__permeability=perm
177          """           self.__permeability_inv=perm_inv
         assigns values to model parameters  
   
         @param f: volumetic sources/sinks  
         @type f: scalar value on the domain (e.g. L{Data})  
         @param g: flux sources/sinks  
         @type g: vector values on the domain (e.g. L{Data})  
         @param location_of_fixed_pressure: mask for locations where pressure is fixed  
         @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})  
         @param location_of_fixed_flux:  mask for locations where flux is fixed.  
         @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})  
         @param permeability: permeability tensor. If scalar C{s} is given the tensor with  
                              C{s} on the main diagonal is used. If vector C{v} is given the tensor with  
                              C{v} on the main diagonal is used.  
         @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})  
   
         @note: the values of parameters which are not set by calling C{setValue} are not altered.  
         @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)  
                or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal  
                is along the M{x_i} axis.  
         """  
         if f !=None:  
            f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))  
            if f.isEmpty():  
                f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))  
            else:  
                if f.getRank()>0: raise ValueError,"illegal rank of f."  
            self.__f=f  
         if g !=None:  
            g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))  
            if g.isEmpty():  
              g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))  
            else:  
              if not g.getShape()==(self.domain.getDim(),):  
                raise ValueError,"illegal shape of g"  
            self.__g=g  
   
         if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)  
         if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)  
   
         if permeability!=None:  
            perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))  
            if perm.getRank()==0:  
                perm=perm*util.kronecker(self.domain.getDim())  
            elif perm.getRank()==1:  
                perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm  
                for i in range(self.domain.getDim()): perm[i,i]=perm2[i]  
            elif perm.getRank()==2:  
               pass  
            else:  
               raise ValueError,"illegal rank of permeability."  
            self.__permeability=perm  
            self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))  
   
     def setTolerance(self,rtol=1e-4):  
         """  
         sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if  
   
         M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }  
   
         where C{atol} is an absolut tolerance (see L{setAbsoluteTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.  
   
         @param rtol: relative tolerance for the pressure  
         @type rtol: non-negative C{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: C{float}  
         """  
         return self.__rtol  
   
     def setAbsoluteTolerance(self,atol=0.):  
         """  
         sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if  
   
         M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }  
   
         where C{rtol} is an absolut tolerance (see L{setTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.  
   
         @param atol: absolute tolerance for the pressure  
         @type atol: non-negative C{float}  
         """  
         if atol<0:  
             raise ValueError,"Absolute tolerance needs to be non-negative."  
         self.__atol=atol  
     def getAbsoluteTolerance(self):  
        """  
        returns the absolute tolerance  
         
        @return: current absolute tolerance  
        @rtype: C{float}  
        """  
        return self.__atol  
     def getSubProblemTolerance(self):  
     """  
     Returns a suitable subtolerance  
     @type: C{float}  
     """  
     return max(util.EPSILON**(0.75),self.getTolerance()**2)  
     def setSubProblemTolerance(self):  
          """  
          Sets the relative tolerance to solve the subproblem(s) if subtolerance adaption is selected.  
          """  
      if self.__adaptSubTolerance:  
          sub_tol=self.getSubProblemTolerance()  
              self.getSolverOptionsFlux().setTolerance(sub_tol)  
          self.getSolverOptionsFlux().setAbsoluteTolerance(0.)  
          self.getSolverOptionsPressure().setTolerance(sub_tol)  
          self.getSolverOptionsPressure().setAbsoluteTolerance(0.)  
          if self.verbose: print "DarcyFlux: relative subtolerance is set to %e."%sub_tol  
   
     def solve(self,u0,p0, max_iter=100, verbose=False, max_num_corrections=10):  
          """  
          solves the problem.  
   
          The iteration is terminated if the residual norm is less then self.getTolerance().  
   
          @param u0: initial guess for the flux. At locations in the domain marked by C{location_of_fixed_flux} the value of C{u0} is kept unchanged.  
          @type u0: vector value on the domain (e.g. L{Data}).  
          @param p0: initial guess for the pressure. At locations in the domain marked by C{location_of_fixed_pressure} the value of C{p0} is kept unchanged.  
          @type p0: scalar value on the domain (e.g. L{Data}).  
          @param verbose: if set some information on iteration progress are printed  
          @type verbose: C{bool}  
          @return: flux and pressure  
          @rtype: C{tuple} of L{Data}.  
   
          @note: The problem is solved as a least squares form  
   
          M{(I+D^*D)u+Qp=D^*f+g}  
          M{Q^*u+Q^*Qp=Q^*g}  
   
          where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.  
          We eliminate the flux form the problem by setting  
   
          M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux  
   
          form the first equation. Inserted into the second equation we get  
   
          M{Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))} with p=p0  on location_of_fixed_pressure  
   
          which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step  
          PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.  
          """  
          self.verbose=verbose  
          rtol=self.getTolerance()  
          atol=self.getAbsoluteTolerance()  
      self.setSubProblemTolerance()  
178            
179           num_corrections=0           #====================
180           converged=False           self.__pde_p.setValue(A=self.__permeability)
181           p=p0           if self.solver  == self.EVAL:
182           norm_r=None                pass # no extra work required
183           while not converged:           elif self.solver  == self.POST:
184                 v=self.getFlux(p, fixed_flux=u0)                k=util.kronecker(self.domain.getDim())
185                 Qp=self.__Q(p)                self.omega = self.w*util.length(perm_inv)*self.l*self.domain.getSize()
186                 norm_v=self.__L2(v)                #self.__pde_v.setValue(D=self.__permeability_inv, A=self.omega*util.outer(k,k))
187                 norm_Qp=self.__L2(Qp)                self.__pde_v.setValue(D=self.__permeability_inv, A_reduced=self.omega*util.outer(k,k))
188                 if norm_v == 0.:           elif self.solver  == self.SMOOTH:
189                    if norm_Qp == 0.:              self.__pde_v.setValue(D=self.__permeability_inv)
190                       return v,p  
191                    else:        if g != None:
192                      fac=norm_Qp          g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
193                 else:          if g.isEmpty():
194                    if norm_Qp == 0.:               g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
                     fac=norm_v  
                   else:  
                     fac=2./(1./norm_v+1./norm_Qp)  
                ATOL=(atol+rtol*fac)  
                if self.verbose:  
                     print "DarcyFlux: L2 norm of v = %e."%norm_v  
                     print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp  
                     print "DarcyFlux: L2 defect u = %e."%(util.integrate(util.length(self.__g-util.interpolate(v,Function(self.domain))-Qp)**2)**(0.5),)  
                     print "DarcyFlux: L2 defect div(v) = %e."%(util.integrate((self.__f-util.div(v))**2)**(0.5),)  
                     print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL  
                if norm_r == None or norm_r>ATOL:  
                    if num_corrections>max_num_corrections:  
                          raise ValueError,"maximum number of correction steps reached."  
                    p,r, norm_r=PCG(self.__g-util.interpolate(v,Function(self.domain))-Qp,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.5*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)  
                    num_corrections+=1  
                else:  
                    converged=True  
          return v,p  
     def __L2(self,v):  
          return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))  
   
     def __Q(self,p):  
           return util.tensor_mult(self.__permeability,util.grad(p))  
   
     def __Aprod(self,dp):  
           if self.getSolverOptionsFlux().isVerbose(): print "DarcyFlux: Applying operator"  
           Qdp=self.__Q(dp)  
           self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())  
           du=self.__pde_v.getSolution()  
           # self.__pde_v.getOperator().saveMM("proj.mm")  
           return Qdp+du  
     def __inner_GMRES(self,r,s):  
          return util.integrate(util.inner(r,s))  
   
     def __inner_PCG(self,p,r):  
          return util.integrate(util.inner(self.__Q(p), r))  
   
     def __Msolve_PCG(self,r):  
       if self.getSolverOptionsPressure().isVerbose(): print "DarcyFlux: Applying preconditioner"  
           self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())  
           # self.__pde_p.getOperator().saveMM("prec.mm")  
           return self.__pde_p.getSolution()  
   
     def getFlux(self,p=None, fixed_flux=Data()):  
         """  
         returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}  
         on locations where C{location_of_fixed_flux} is positive (see L{setValue}).  
         Note that C{g} and C{f} are used, see L{setValue}.  
   
         @param p: pressure.  
         @type p: scalar value on the domain (e.g. L{Data}).  
         @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.  
         @type fixed_flux: vector values on the domain (e.g. L{Data}).  
         @param tol: relative tolerance to be used.  
         @type tol: positive C{float}.  
         @return: flux  
         @rtype: L{Data}  
         @note: the method uses the least squares solution M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}}  
                for the permeability M{k_{ij}}  
         """  
     self.setSubProblemTolerance()  
         g=self.__g  
         f=self.__f  
         self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)  
         if p == None:  
            self.__pde_v.setValue(Y=g)  
195          else:          else:
196             self.__pde_v.setValue(Y=g-self.__Q(p))               if not g.getShape()==(self.domain.getDim(),): raise ValueError("illegal shape of g")
197          return self.__pde_v.getSolution()          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:
201             f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
202             if f.isEmpty():      
203                 f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
204             else:
205                 if f.getRank()>0: raise ValueError("illegal rank of f.")
206             self.__f=f
207    
208       def getSolverOptionsFlux(self):
209          """
210          Returns the solver options used to solve the flux problems
211          :return: `SolverOptions`
212          """
213          if self.__pde_v == None:
214              return None
215          else:
216              return self.__pde_v.getSolverOptions()
217          
218       def setSolverOptionsFlux(self, options=None):
219          """
220          Sets the solver options used to solve the flux problems
221          If ``options`` is not present, the options are reset to default
222          :param options: `SolverOptions`
223          """
224          if not self.__pde_v == None:
225              self.__pde_v.setSolverOptions(options)
226        
227       def getSolverOptionsPressure(self):
228          """
229          Returns the solver options used to solve the pressure problems
230          :return: `SolverOptions`
231          """
232          return self.__pde_p.getSolverOptions()
233          
234       def setSolverOptionsPressure(self, options=None):
235          """
236          Sets the solver options used to solve the pressure problems
237          If ``options`` is not present, the options are reset to default
238          
239          :param options: `SolverOptions`
240          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
241          """
242          return self.__pde_p.setSolverOptions(options)
243          
244       def solve(self, u0, p0):
245          """
246          solves the problem.
247          
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.
249          :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.
251          :type p0: scalar value on the domain (e.g. `escript.Data`).
252          :return: flux and pressure
253          :rtype: ``tuple`` of `escript.Data`.
254    
255          """
256          p0=util.interpolate(p0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
257          if self.ref_point_id == None:
258              p_ref=0
259          else:
260              p_ref=p0.getTupleForGlobalDataPoint(*self.ref_point_id)[0]
261          p0_hydrostatic=p_ref+util.inner(self.__permeability_invXg_ref, self.__pde_p.getFunctionSpaceForCoefficient("q").getX() - self.ref_point)
262          g_2=self.__g - util.tensor_mult(self.__permeability, self.__permeability_invXg_ref * self.perm_scale)
263          self.__pde_p.setValue(X=g_2 * 1./self.perm_scale,
264                                Y=self.__f * 1./self.perm_scale,
265                                y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux * 1./self.perm_scale ),
266                                r=p0 - p0_hydrostatic)
267          pp=self.__pde_p.getSolution()
268          u = self._getFlux(pp, u0)
269          return u,pp + p0_hydrostatic
270          
271       def getFlux(self,p, u0=None):
272            """
273            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`).
275            Notice that ``g`` is used, see `setValue`.
276    
277            :param p: pressure.
278            :type p: scalar value on the domain (e.g. `escript.Data`).
279            :param u0: flux on the locations of the domain marked be ``location_of_fixed_flux``.
280            :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
281            :return: flux
282            :rtype: `escript.Data`
283            """
284            p=util.interpolate(p, self.__pde_p.getFunctionSpaceForCoefficient("q"))
285            if self.ref_point_id == None:
286                p_ref=0
287            else:
288                p_ref=p.getTupleForGlobalDataPoint(*self.ref_point_id)[0]
289            p_hydrostatic=p_ref+util.inner(self.__permeability_invXg_ref, self.__pde_p.getFunctionSpaceForCoefficient("q").getX() - self.ref_point)
290            return self._getFlux(p-p_hydrostatic, u0)
291    
292       def _getFlux(self, pp, u0=None):
293            """
294            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:
311                   self.__pde_v.setValue(r=escript.Data())
312                else:
313                   if not isinstance(u0, escript.Data) : u0 = escript.Vector(u0, escript.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(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
319       """       """
320       solves       solves
# Line 379  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,adaptSubTolerance=True, **kwargs):       def __init__(self,domain,**kwargs):
340           """           """
341           initialize the Stokes Problem           initialize the Stokes Problem
342    
343           @param domain: domain of the problem. The approximation order needs to be two.           The approximation spaces used for velocity (=Solution(domain)) and pressure (=ReducedSolution(domain)) must be
344           @type domain: L{Domain}           LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
345       @param adaptSubTolerance: If True the tolerance for subproblem is set automatically.           with macro elements for the pressure.
346       @type adaptSubTolerance: C{bool}  
347           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.           :param domain: domain of the problem.
348             :type domain: `Domain`
349           """           """
350           HomogeneousSaddlePointProblem.__init__(self,adaptSubTolerance=adaptSubTolerance,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
351           self.domain=domain           self.domain=domain
352           self.vol=util.integrate(1.,Function(self.domain))           self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
353           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_v.setSymmetryOn()
          self.__pde_u.setSymmetryOn()  
354            
355           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
356           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
# Line 402  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):
365           """           """
366       returns the solver options used  solve the equation for velocity.       returns the solver options used  solve the equation for velocity.
367            
368       @rtype: L{SolverOptions}       :rtype: `SolverOptions`
369       """       """
370       return self.__pde_u.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.
374            
375       @param options: new solver  options       :param options: new solver  options
376       @type options: L{SolverOptions}       :type options: `SolverOptions`
377       """       """
378           self.__pde_u.setSolverOptions(options)           self.__pde_v.setSolverOptions(options)
379       def getSolverOptionsPressure(self):       def getSolverOptionsPressure(self):
380           """           """
381       returns the solver options used  solve the equation for pressure.       returns the solver options used  solve the equation for pressure.
382       @rtype: L{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: L{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           """           """
395       set the solver options for solving the equation to project the divergence of       set the solver options for solving the equation to project the divergence of
396       the velocity onto the function space of presure.       the velocity onto the function space of presure.
397            
398       @param options: new solver options       :param options: new solver options
399       @type options: L{SolverOptions}       :type options: `SolverOptions`
400       """       """
401       self.__pde_prec.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
405       the velocity onto the function space of presure.       the velocity onto the function space of presure.
406            
407       @rtype: L{SolverOptions}       :rtype: `SolverOptions`
408       """       """
409       return self.__pde_prec.getSolverOptions()           return self.__pde_proj.getSolverOptions()
410       def setSubProblemTolerance(self):  
411         def updateStokesEquation(self, v, p):
412           """           """
413       Updates the tolerance for subproblems           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 to model parameters.
415           """           """
416       if self.adaptSubTolerance():           pass
417               sub_tol=self.getSubProblemTolerance()       def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
418           self.getSolverOptionsDiv().setTolerance(sub_tol)          """
419           self.getSolverOptionsDiv().setAbsoluteTolerance(0.)          assigns new values to the model parameters.
420           self.getSolverOptionsPressure().setTolerance(sub_tol)  
421           self.getSolverOptionsPressure().setAbsoluteTolerance(0.)          :param f: external force
422           self.getSolverOptionsVelocity().setTolerance(sub_tol)          :type f: `Vector` object in `FunctionSpace` `Function` or similar
423           self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)          :param fixed_u_mask: mask of locations with fixed velocity.
424                    :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
425            :param eta: viscosity
426            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
427            :param surface_stress: normal surface stress
428            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
429            :param stress: initial stress
430        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
431            """
432            if eta !=None:
433                k=util.kronecker(self.domain.getDim())
434                kk=util.outer(k,k)
435                self.eta=util.interpolate(eta, escript.Function(self.domain))
436                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)))
438            if restoration_factor!=None:
439                n=self.domain.getNormal()
440                self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
441            if fixed_u_mask!=None:
442                self.__pde_v.setValue(q=fixed_u_mask)
443            if f!=None: self.__f=f
444            if surface_stress!=None: self.__surface_stress=surface_stress
445            if stress!=None: self.__stress=stress
446    
447       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):       def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
448          """          """
449          assigns values to the model parameters          assigns values to the model parameters
450    
451          @param f: external force          :param f: external force
452          @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar          :type f: `Vector` object in `FunctionSpace` `Function` or similar
453          @param fixed_u_mask: mask of locations with fixed velocity.          :param fixed_u_mask: mask of locations with fixed velocity.
454          @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar          :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
455          @param eta: viscosity          :param eta: viscosity
456          @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar          :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
457          @param surface_stress: normal surface stress          :param surface_stress: normal surface stress
458          @type eta: L{Vector} object on L{FunctionSpace} L{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: L{Tensor} object on L{FunctionSpace} L{Function} or similar          :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
         @note: All values needs to be set.  
461          """          """
462          self.eta=eta          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
         A =self.__pde_u.createCoefficient("A")  
     self.__pde_u.setValue(A=Data())  
         for i in range(self.domain.getDim()):  
         for j in range(self.domain.getDim()):  
             A[i,j,j,i] += 1.  
             A[i,j,i,j] += 1.  
     self.__pde_prec.setValue(D=1/self.eta)  
         self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask)  
         self.__f=f  
         self.__surface_stress=surface_stress  
         self.__stress=stress  
463    
464       def Bv(self,v):       def Bv(self,v,tol):
465           """           """
466           returns inner product of element p and div(v)           returns inner product of element p and div(v)
467    
468           @param p: a pressure increment           :param v: a residual
469           @param v: a residual           :return: inner product of element p and div(v)
470           @return: inner product of element p and div(v)           :rtype: ``float``
471           @rtype: C{float}           """
472           """           self.__pde_proj.setValue(Y=-util.div(v))
473           self.__pde_proj.setValue(Y=-util.div(v))           self.getSolverOptionsDiv().setTolerance(tol)
474           return self.__pde_proj.getSolution()           self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
475             out=self.__pde_proj.getSolution()
476             return out
477    
478       def inner_pBv(self,p,Bv):       def inner_pBv(self,p,Bv):
479           """           """
480           returns inner product of element p and Bv=-div(v)           returns inner product of element p and Bv=-div(v)
481    
482           @param p: a pressure increment           :param p: a pressure increment
483           @param v: a residual           :param Bv: a residual
484           @return: inner product of element p and Bv=-div(v)           :return: inner product of element p and Bv=-div(v)
485           @rtype: C{float}           :rtype: ``float``
486           """           """
487           return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))           return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
488    
489       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
490           """           """
491           Returns inner product of p0 and p1           Returns inner product of p0 and p1
492    
493           @param p0: a pressure           :param p0: a pressure
494           @param p1: a pressure           :param p1: a pressure
495           @return: inner product of p0 and p1           :return: inner product of p0 and p1
496           @rtype: C{float}           :rtype: ``float``
497           """           """
498           s0=util.interpolate(p0/self.eta,Function(self.domain))           s0=util.interpolate(p0, escript.Function(self.domain))
499           s1=util.interpolate(p1/self.eta,Function(self.domain))           s1=util.interpolate(p1, escript.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):
503           """           """
504           returns the norm of v           returns the norm of v
505    
506           @param v: a velovity           :param v: a velovity
507           @return: norm of v           :return: norm of v
508           @rtype: non-negative C{float}           :rtype: non-negative ``float``
509           """           """
510           return util.sqrt(util.integrate(util.length(util.grad(v))))           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
511    
512    
513       def getV(self, p, v0):       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 v0: a initial guess for the value v to return.           :param v: a initial guess for the value v to return.
519           @return: v given as M{v= A^{-1} (f-B^*p)}           :return: dv given as *Adv=(f-Av-B^*p)*
520           """           """
521           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)           self.updateStokesEquation(v,p)
522             self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
523             self.getSolverOptionsVelocity().setTolerance(tol)
524             self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
525           if self.__stress.isEmpty():           if self.__stress.isEmpty():
526              self.__pde_u.setValue(X=p*util.kronecker(self.domain))              self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
527           else:           else:
528              self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))              self.__pde_v.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
529           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
530           return  out           return  out
531    
532       def norm_Bv(self,Bv):       def norm_Bv(self,Bv):
533          """          """
534          Returns Bv (overwrite).          Returns Bv (overwrite).
535    
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,Function(self.domain))**2))          return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
540    
541       def solve_AinvBt(self,p):       def solve_AinvBt(self,p, tol):
542           """           """
543           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}           Solves *Av=B^*p* with accuracy `tol`
544    
545           @param p: a pressure increment           :param p: a pressure increment
546           @return: the solution of M{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_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))           self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
550           out=self.__pde_u.getSolution()           out=self.__pde_v.getSolution()
551           return  out           return  out
552    
553       def solve_prec(self,Bv):       def solve_prec(self,Bv, tol):
554           """           """
555           applies preconditioner for for M{BA^{-1}B^*} to M{Bv}           applies preconditioner for for *BA^{-1}B^** to *Bv*
556           with accuracy L{self.getSubProblemTolerance()}           with accuracy ``self.getSubProblemTolerance()``
557    
558           @param v: velocity increment           :param Bv: velocity increment
559           @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{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           return self.__pde_prec.getSolution()           self.getSolverOptionsPressure().setTolerance(tol)
564             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
565             out=self.__pde_prec.getSolution()
566             return out

Legend:
Removed from v.2548  
changed lines
  Added in v.3990

  ViewVC Help
Powered by ViewVC 1.1.26