/[escript]/trunk/escript/py_src/flows.py
ViewVC logotype

Diff of /trunk/escript/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 3852 by jfenwick, Thu Mar 1 05:34:52 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             self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale
100    
101          elif self.solver  == self.SMOOTH:
102             self.__pde_v=LinearPDESystem(domain)
103             self.__pde_v.setSymmetryOn()
104             if self.useReduced: self.__pde_v.setReducedOrderOn()
105             if self.verbose: print("DarcyFlow: flux smoothing is used.")
106             self.w=0
107    
108          self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
109          self.__g=escript.Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
110          self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
111          self.location_of_fixed_flux = escript.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
112          self.perm_scale=1.
113        
114                    
115          @param v: flux.     def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
116          @type v: vector-valued on the domain (e.g. L{Data}).        """
117          @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  
118    
119          M{<(Q^*Q)^{-1} r,r> <= ATOL}        :param f: volumetic sources/sinks
120          :type f: scalar value on the domain (e.g. `escript.Data`)
121          :param g: flux sources/sinks
122          :type g: vector values on the domain (e.g. `escript.Data`)
123          :param location_of_fixed_pressure: mask for locations where pressure is fixed
124          :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
125          :param location_of_fixed_flux:  mask for locations where flux is fixed.
126          :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
127          :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
128          :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
129    
130          :note: the values of parameters which are not set by calling ``setValue`` are not altered.
131          :note: at any point on the boundary of the domain the pressure
132                 (``location_of_fixed_pressure`` >0) or the normal component of the
133                 flux (``location_of_fixed_flux[i]>0``) if direction of the normal
134                 is along the *x_i* axis.
135    
136          holds. M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}        """
137          if location_of_fixed_pressure!=None:
138               self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)
139               self.__pde_p.setValue(q=self.location_of_fixed_pressure)
140          if location_of_fixed_flux!=None:
141              self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
142              if not self.__pde_v == None: self.__pde_v.setValue(q=self.location_of_fixed_flux)
143                
144          if permeability!=None:
145        
146             perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
147             self.perm_scale=util.Lsup(util.length(perm))
148             if self.verbose: print(("DarcyFlow: permeability scaling factor = %e."%self.perm_scale))
149             perm=perm*(1./self.perm_scale)
150            
151             if perm.getRank()==0:
152    
153          @param atol: absolute tolerance for the pressure              perm_inv=(1./perm)
154          @type atol: non-negative C{float}              perm_inv=perm_inv*util.kronecker(self.domain.getDim())
155          @param rtol: relative tolerance for the pressure              perm=perm*util.kronecker(self.domain.getDim())
156          @type rtol: non-negative C{float}          
157          @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          
158          L{getPressure} method or use  the value from a previous time step.           elif perm.getRank()==2:
159          @type p_ref: scalar value on the domain (e.g. L{Data}).              perm_inv=util.inverse(perm)
160          @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:
161          L{getFlux} method or use  the value from a previous time step.              raise ValueError("illegal rank of permeability.")
162          @type v_ref: vector-valued on the domain (e.g. L{Data}).          
163          @return: used absolute tolerance.           self.__permeability=perm
164          @rtype: positive C{float}           self.__permeability_inv=perm_inv
165          """      
166          g=self.__g           #====================
167          if not v_ref == None:           self.__pde_p.setValue(A=self.__permeability)
168             f1=util.integrate(util.length(util.interpolate(g-v_ref,Function(self.domain)))**2)           if self.solver  == self.EVAL:
169          else:                pass # no extra work required
170             f1=util.integrate(util.length(util.interpolate(g))**2)           elif self.solver  == self.POST:
171          if not p_ref == None:                k=util.kronecker(self.domain.getDim())
172             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()
173                  self.__pde_v.setValue(D=self.__permeability_inv, A=self.omega*util.outer(k,k))
174             elif self.solver  == self.SMOOTH:
175                self.__pde_v.setValue(D=self.__permeability_inv)
176    
177          if g != None:
178            g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
179            if g.isEmpty():
180                 g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
181          else:          else:
182             f2=0               if not g.getShape()==(self.domain.getDim(),): raise ValueError("illegal shape of g")
183          self.__ATOL= atol + rtol * util.sqrt(max(f1,f2))          self.__g=g
184          if self.__ATOL<=0:        if f !=None:
185             raise ValueError,"Positive tolerance (=%e) is expected."%self.__ATOL           f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
186          return self.__ATOL           if f.isEmpty():      
187                         f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
188      def getTolerance(self):           else:
189          """               if f.getRank()>0: raise ValueError("illegal rank of f.")
190          returns the current tolerance.           self.__f=f
     
         @return: used absolute tolerance.  
         @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  
191    
192           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):
193                  """
194           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
195           PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.        :return: `SolverOptions`
196           """        """
197           self.verbose=verbose        if self.__pde_v == None:
198           self.show_details= show_details and self.verbose            return None
199           self.__pde_v.setTolerance(sub_rtol)        else:
200           self.__pde_p.setTolerance(sub_rtol)            return self.__pde_v.getSolverOptions()
201           ATOL=self.getTolerance()        
202           if self.verbose: print "DarcyFlux: absolute tolerance = %e"%ATOL     def setSolverOptionsFlux(self, options=None):
203           #########################################################################################################################        """
204           #        Sets the solver options used to solve the flux problems
205           #   we solve:        If ``options`` is not present, the options are reset to default
206           #          :param options: `SolverOptions`
207           #      Q^*(I-(I+D^*D)^{-1})Q dp =  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) )        """
208           #        if not self.__pde_v == None:
209           #   residual is            self.__pde_v.setSolverOptions(options)
210           #      
211           #    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):
212           #        """
213           #        with v = (I+D^*D)^{-1} (D^*f+g-Qp) including BC        Returns the solver options used to solve the pressure problems
214           #        :return: `SolverOptions`
215           #    we use (g - Qp, v) to represent the residual. not that        """
216           #        return self.__pde_p.getSolverOptions()
217           #    dr(dp)=( -Q(dp), dv) with dv = - (I+D^*D)^{-1} Q(dp)        
218           #     def setSolverOptionsPressure(self, options=None):
219           #   while the initial residual is        """
220           #        Sets the solver options used to solve the pressure problems
221           #      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
222           #          
223           d0=self.__g-util.tensor_mult(self.__permeability,util.grad(p0))        :param options: `SolverOptions`
224           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.
225           v00=self.__pde_v.getSolution(verbose=show_details)        """
226           if self.verbose: print "DarcyFlux: range of initial flux = ",util.inf(v00), util.sup(v00)        return self.__pde_p.setSolverOptions(options)
227           self.__pde_v.setValue(r=Data())        
228           # start CG     def solve(self, u0, p0):
229           r=ArithmeticTuple(d0, v00)        """
230           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.
231           return r[1],p        
232          :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.
233      def __Aprod_PCG(self,dp):        :type u0: vector value on the domain (e.g. `escript.Data`).
234            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.
235            #  -dr(dp) = (Qdp,du) where du = (I+D^*D)^{-1} (Qdp)        :type p0: scalar value on the domain (e.g. `escript.Data`).
236            mQdp=util.tensor_mult(self.__permeability,util.grad(dp))        :return: flux and pressure
237            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)  
238    
 class StokesProblemCartesian(HomogeneousSaddlePointProblem):  
239        """        """
240        solves        self.__pde_p.setValue(X=self.__g * 1./self.perm_scale,
241                                Y=self.__f * 1./self.perm_scale,
242                                y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux * 1./self.perm_scale ),
243                                r=p0)
244          p=self.__pde_p.getSolution()
245          u = self.getFlux(p, u0)
246          return u,p
247          
248       def getFlux(self,p, u0=None):
249            """
250            returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
251            on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
252            Notice that ``g`` is used, see `setValue`.
253    
254            :param p: pressure.
255            :type p: scalar value on the domain (e.g. `escript.Data`).
256            :param u0: flux on the locations of the domain marked be ``location_of_fixed_flux``.
257            :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
258            :return: flux
259            :rtype: `escript.Data`
260            """
261            if self.solver  == self.EVAL:
262               u = self.__g-self.perm_scale * util.tensor_mult(self.__permeability,util.grad(p))
263            elif self.solver  == self.POST or self.solver  == self.SMOOTH:
264                self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g * 1./self.perm_scale)-util.grad(p))
265                if u0 == None:
266                   self.__pde_v.setValue(r=escript.Data())
267                else:
268                   if not isinstance(u0, escript.Data) : u0 = escript.Vector(u0, escript.Solution(self.domain))
269                   self.__pde_v.setValue(r=1./self.perm_scale * u0)
270                   u= self.__pde_v.getSolution() * self.perm_scale
271            return u
272          
273    class StokesProblemCartesian(HomogeneousSaddlePointProblem):
274         """
275         solves
276    
277            -(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}
278                  u_{i,i}=0                  u_{i,i}=0
# Line 333  class StokesProblemCartesian(Homogeneous Line 280  class StokesProblemCartesian(Homogeneous
280            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
281            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
282    
283        if surface_stress is not given 0 is assumed.       if surface_stress is not given 0 is assumed.
284    
285        typical usage:       typical usage:
286    
287              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
288              sp.setTolerance()              sp.setTolerance()
289              sp.initialize(...)              sp.initialize(...)
290              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
291        """              sp.setStokesEquation(...) # new values for some parameters
292        def __init__(self,domain,**kwargs):              v1,p1=sp.solve(v,p)
293         """
294         def __init__(self,domain,**kwargs):
295           """           """
296           initialize the Stokes Problem           initialize the Stokes Problem
297    
298           @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
299           @type domain: L{Domain}           LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
300           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.           with macro elements for the pressure.
301    
302             :param domain: domain of the problem.
303             :type domain: `Domain`
304           """           """
305           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
306           self.domain=domain           self.domain=domain
307           self.vol=util.integrate(1.,Function(self.domain))           self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
308           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_v.setSymmetryOn()
309           self.__pde_u.setSymmetryOn()      
          # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)  
          # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)  
               
310           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
311           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
          # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)  
312           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
313    
314           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
315           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
316             self.__pde_proj.setValue(D=1)
317           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
          self.__pde_proj.setValue(D=1.)  
318    
319        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):       def getSolverOptionsVelocity(self):
320             """
321         returns the solver options used  solve the equation for velocity.
322        
323         :rtype: `SolverOptions`
324         """
325             return self.__pde_v.getSolverOptions()
326         def setSolverOptionsVelocity(self, options=None):
327             """
328         set the solver options for solving the equation for velocity.
329        
330         :param options: new solver  options
331         :type options: `SolverOptions`
332         """
333             self.__pde_v.setSolverOptions(options)
334         def getSolverOptionsPressure(self):
335             """
336         returns the solver options used  solve the equation for pressure.
337         :rtype: `SolverOptions`
338         """
339             return self.__pde_prec.getSolverOptions()
340         def setSolverOptionsPressure(self, options=None):
341             """
342         set the solver options for solving the equation for pressure.
343         :param options: new solver  options
344         :type options: `SolverOptions`
345         """
346             self.__pde_prec.setSolverOptions(options)
347    
348         def setSolverOptionsDiv(self, options=None):
349             """
350         set the solver options for solving the equation to project the divergence of
351         the velocity onto the function space of presure.
352        
353         :param options: new solver options
354         :type options: `SolverOptions`
355         """
356             self.__pde_proj.setSolverOptions(options)
357         def getSolverOptionsDiv(self):
358             """
359         returns the solver options for solving the equation to project the divergence of
360         the velocity onto the function space of presure.
361        
362         :rtype: `SolverOptions`
363         """
364             return self.__pde_proj.getSolverOptions()
365    
366         def updateStokesEquation(self, v, p):
367             """
368             updates the Stokes equation to consider dependencies from ``v`` and ``p``
369             :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values to model parameters.
370             """
371             pass
372         def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
373            """
374            assigns new values to the model parameters.
375    
376            :param f: external force
377            :type f: `Vector` object in `FunctionSpace` `Function` or similar
378            :param fixed_u_mask: mask of locations with fixed velocity.
379            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
380            :param eta: viscosity
381            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
382            :param surface_stress: normal surface stress
383            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
384            :param stress: initial stress
385        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
386            """
387            if eta !=None:
388                k=util.kronecker(self.domain.getDim())
389                kk=util.outer(k,k)
390                self.eta=util.interpolate(eta, escript.Function(self.domain))
391                self.__pde_prec.setValue(D=1/self.eta)
392                self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
393            if restoration_factor!=None:
394                n=self.domain.getNormal()
395                self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
396            if fixed_u_mask!=None:
397                self.__pde_v.setValue(q=fixed_u_mask)
398            if f!=None: self.__f=f
399            if surface_stress!=None: self.__surface_stress=surface_stress
400            if stress!=None: self.__stress=stress
401    
402         def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
403          """          """
404          assigns values to the model parameters          assigns values to the model parameters
405    
406          @param f: external force          :param f: external force
407          @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar          :type f: `Vector` object in `FunctionSpace` `Function` or similar
408          @param fixed_u_mask: mask of locations with fixed velocity.          :param fixed_u_mask: mask of locations with fixed velocity.
409          @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
410          @param eta: viscosity          :param eta: viscosity
411          @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar          :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
412          @param surface_stress: normal surface stress          :param surface_stress: normal surface stress
413          @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
414          @param stress: initial stress          :param stress: initial stress
415      @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.  
   
416          """          """
417          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  
418    
419        def B(self,v):       def Bv(self,v,tol):
420          """           """
421          returns div(v)           returns inner product of element p and div(v)
         @rtype: equal to the type of p  
422    
423          @note: boundary conditions on p should be zero!           :param v: a residual
424          """           :return: inner product of element p and div(v)
425          if self.show_details: print "apply divergence:"           :rtype: ``float``
426          self.__pde_proj.setValue(Y=-util.div(v))           """
427          self.__pde_proj.setTolerance(self.getSubProblemTolerance())           self.__pde_proj.setValue(Y=-util.div(v))
428          return self.__pde_proj.getSolution(verbose=self.show_details)           self.getSolverOptionsDiv().setTolerance(tol)
429             self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
430             out=self.__pde_proj.getSolution()
431             return out
432    
433        def inner_pBv(self,p,Bv):       def inner_pBv(self,p,Bv):
434           """           """
435           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}  
436    
437           @rtype: equal to the type of p           :param p: a pressure increment
438             :param Bv: a residual
439             :return: inner product of element p and Bv=-div(v)
440             :rtype: ``float``
441           """           """
442           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)  
443    
444        def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
445           """           """
446           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}  
447    
448           @rtype: equal to the type of p           :param p0: a pressure
449             :param p1: a pressure
450             :return: inner product of p0 and p1
451             :rtype: ``float``
452           """           """
453           s0=util.interpolate(p0/self.eta,Function(self.domain))           s0=util.interpolate(p0, escript.Function(self.domain))
454           s1=util.interpolate(p1/self.eta,Function(self.domain))           s1=util.interpolate(p1, escript.Function(self.domain))
455           return util.integrate(s0*s1)           return util.integrate(s0*s1)
456    
457        def inner_v(self,v0,v1):       def norm_v(self,v):
458           """           """
459           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}  
460    
461           @rtype: equal to the type of v           :param v: a velovity
462             :return: norm of v
463             :rtype: non-negative ``float``
464           """           """
465       gv0=util.grad(v0)           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
466       gv1=util.grad(v1)  
          return util.integrate(util.inner(gv0,gv1))  
467    
468        def solve_A(self,u,p):       def getDV(self, p, v, tol):
469           """           """
470           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           return the value for v for a given p
471    
472             :param p: a pressure
473             :param v: a initial guess for the value v to return.
474             :return: dv given as *Adv=(f-Av-B^*p)*
475           """           """
476           if self.show_details: print "solve for velocity:"           self.updateStokesEquation(v,p)
477           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
478             self.getSolverOptionsVelocity().setTolerance(tol)
479             self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
480           if self.__stress.isEmpty():           if self.__stress.isEmpty():
481              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)))
482           else:           else:
483              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)))
484           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_v.getSolution()
485           return  out           return  out
486    
487        def solve_prec(self,p):       def norm_Bv(self,Bv):
488           if self.show_details: print "apply preconditioner:"          """
489           self.__pde_prec.setTolerance(self.getSubProblemTolerance())          Returns Bv (overwrite).
490           self.__pde_prec.setValue(Y=p)  
491           q=self.__pde_prec.getSolution(verbose=self.show_details)          :rtype: equal to the type of p
492           return q          :note: boundary conditions on p should be zero!
493            """
494            return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
495    
496         def solve_AinvBt(self,p, tol):
497             """
498             Solves *Av=B^*p* with accuracy `tol`
499    
500             :param p: a pressure increment
501             :return: the solution of *Av=B^*p*
502             :note: boundary conditions on v should be zero!
503             """
504             self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
505             out=self.__pde_v.getSolution()
506             return  out
507    
508         def solve_prec(self,Bv, tol):
509             """
510             applies preconditioner for for *BA^{-1}B^** to *Bv*
511             with accuracy `self.getSubProblemTolerance()`
512    
513             :param Bv: velocity increment
514             :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
515             :note: boundary conditions on p are zero.
516             """
517             self.__pde_prec.setValue(Y=Bv)
518             self.getSolverOptionsPressure().setTolerance(tol)
519             self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
520             out=self.__pde_prec.getSolution()
521             return out

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

  ViewVC Help
Powered by ViewVC 1.1.26