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

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

  ViewVC Help
Powered by ViewVC 1.1.26