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

Legend:
Removed from v.2208  
changed lines
  Added in v.3909

  ViewVC Help
Powered by ViewVC 1.1.26