/[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 3620 by gross, Wed Oct 5 06:06:14 2011 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 *  import escript
36  import util  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
   
     M{u_i+k_{ij}*p_{,j} = g_i}  
     M{u_{i,i} = f}  
   
     where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,  
   
     @note: The problem is solved in a least squares formulation.  
     """  
   
     def __init__(self, domain,useReduced=False):  
         """  
         initializes the Darcy flux problem  
         @param domain: domain of the problem  
         @type domain: L{Domain}  
         """  
         self.domain=domain  
         self.__pde_v=LinearPDESystem(domain)  
         if useReduced: self.__pde_v.setReducedOrderOn()  
         self.__pde_v.setSymmetryOn()  
         self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain)))  
         self.__pde_p=LinearSinglePDE(domain)  
         self.__pde_p.setSymmetryOn()  
         if useReduced: self.__pde_p.setReducedOrderOn()  
         self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))  
         self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))  
         self.__ATOL= None  
   
     def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):  
         """  
         assigns values to model parameters  
   
         @param f: volumetic sources/sinks  
         @type f: scalar value on the domain (e.g. L{Data})  
         @param g: flux sources/sinks  
         @type g: vector values on the domain (e.g. L{Data})  
         @param location_of_fixed_pressure: mask for locations where pressure is fixed  
         @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})  
         @param location_of_fixed_flux:  mask for locations where flux is fixed.  
         @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})  
         @param permeability: permeability tensor. If scalar C{s} is given the tensor with  
                              C{s} on the main diagonal is used. If vector C{v} is given the tensor with  
                              C{v} on the main diagonal is used.  
         @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})  
   
         @note: the values of parameters which are not set by calling C{setValue} are not altered.  
         @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)  
                or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal  
                is along the M{x_i} axis.  
         """  
         if f !=None:  
            f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))  
            if f.isEmpty():  
                f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))  
            else:  
                if f.getRank()>0: raise ValueError,"illegal rank of f."  
            self.f=f  
         if g !=None:    
            g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))  
            if g.isEmpty():  
              g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))  
            else:  
              if not g.getShape()==(self.domain.getDim(),):  
                raise ValueError,"illegal shape of g"  
            self.__g=g  
   
         if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)  
         if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)  
   
         if permeability!=None:  
            perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))  
            if perm.getRank()==0:  
                perm=perm*util.kronecker(self.domain.getDim())  
            elif perm.getRank()==1:  
                perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm  
                for i in range(self.domain.getDim()): perm[i,i]=perm2[i]  
            elif perm.getRank()==2:  
               pass  
            else:  
               raise ValueError,"illegal rank of permeability."  
            self.__permeability=perm  
            self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))  
   
   
     def 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)  
   
     def getPressure(self,v=None, fixed_pressure=Data(),tol=1.e-8, show_details=False):  
         """  
         returns the pressure for a given flux C{v} where the pressure is equal to C{fixed_pressure}  
         on locations where C{location_of_fixed_pressure} is positive (see L{setValue}).  
         Note that C{g} is used, see L{setValue}.  
           
         @param v: flux.  
         @type v: vector-valued on the domain (e.g. L{Data}).  
         @param fixed_pressure: pressure on the locations of the domain marked be C{location_of_fixed_pressure}.  
         @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  
   
         M{<(Q^*Q)^{-1} r,r> <= ATOL}  
   
         holds. M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}  
   
         @param atol: absolute tolerance for the pressure  
         @type atol: non-negative C{float}  
         @param rtol: relative tolerance for the pressure  
         @type rtol: non-negative C{float}  
         @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  
         L{getPressure} method or use  the value from a previous time step.  
         @type p_ref: scalar value on the domain (e.g. L{Data}).  
         @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  
         L{getFlux} method or use  the value from a previous time step.  
         @type v_ref: vector-valued on the domain (e.g. L{Data}).  
         @return: used absolute tolerance.  
         @rtype: positive C{float}  
         """  
         g=self.__g  
         if not v_ref == None:  
            f1=util.integrate(util.length(util.interpolate(g-v_ref,Function(self.domain)))**2)  
         else:  
            f1=util.integrate(util.length(util.interpolate(g))**2)  
         if not p_ref == None:  
            f2=util.integrate(util.length(util.tensor_mult(self.__permeability,util.grad(p_ref)))**2)  
         else:  
            f2=0  
         self.__ATOL= atol + rtol * util.sqrt(max(f1,f2))  
         if self.__ATOL<=0:  
            raise ValueError,"Positive tolerance (=%e) is expected."%self.__ATOL  
         return self.__ATOL  
           
     def getTolerance(self):  
         """  
         returns the current tolerance.  
43        
44          @return: used absolute tolerance.     *u_i+k_{ij}*p_{,j} = g_i*
45          @rtype: positive C{float}     *u_{i,i} = f*
46          """    
47          if self.__ATOL==None:     where *p* represents the pressure and *u* the Darcy flux. *k* represents the permeability,
48             raise ValueError,"no tolerance is defined."    
49          return self.__ATOL     :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      def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8):                 where *l* is the length scale, *K* is the inverse of the permeability tensor, and *w* is a positive weighting factor.
52           """     :cvar SMOOTH: global smoothing by solving the PDE *K_{ij} u_j= - p_{,j} + K_{ij} g_j*
53           solves the problem.     """
54       EVAL="EVAL"
55           The iteration is terminated if the residual norm is less then self.getTolerance().     SIMPLE="EVAL"
56       POST="POST"
57           @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.     SMOOTH="SMOOTH"
58           @type u0: vector value on the domain (e.g. L{Data}).     def __init__(self, domain, useReduced=False, solver="EVAL", verbose=False, w=1.):
59           @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.        """
60           @type p0: scalar value on the domain (e.g. L{Data}).        initializes the Darcy flux problem
61           @param sub_rtol: tolerance to be used in the sub iteration. It is recommended that M{sub_rtol<rtol*5.e-3}        :param domain: domain of the problem
62           @type sub_rtol: positive-negative C{float}        :type domain: `Domain`
63           @param verbose: if set some information on iteration progress are printed        :param useReduced: uses reduced oreder on flux and pressure
64           @type verbose: C{bool}        :type useReduced: ``bool``
65           @param show_details:  if set information on the subiteration process are printed.        :param solver: solver method
66           @type show_details: C{bool}        :type solver: in [`DarcyFlow.EVAL`, `DarcyFlow.POST',  `DarcyFlow.SMOOTH' ]
67           @return: flux and pressure        :param verbose: if ``True`` some information on the iteration progress are printed.
68           @rtype: C{tuple} of L{Data}.        :type verbose: ``bool``
69          :param w: weighting factor for `DarcyFlow.POST` solver
70           @note: The problem is solved as a least squares form        :type w: ``float``
71          
72           M{(I+D^*D)u+Qp=D^*f+g}        """
73           M{Q^*u+Q^*Qp=Q^*g}        if not solver in [DarcyFlow.EVAL, DarcyFlow.POST,  DarcyFlow.SMOOTH ] :
74              raise ValueError,"unknown solver %d."%solver
75    
76           where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.        self.domain=domain
77           We eliminate the flux form the problem by setting        self.solver=solver
78          self.useReduced=useReduced
79          self.verbose=verbose
80          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       def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
116          """
117          assigns values to model parameters
118    
119           M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux        :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           form the first equation. Inserted into the second equation we get        """
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           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          perm_inv=(1./perm)
154            perm_inv=perm_inv*util.kronecker(self.domain.getDim())
155            perm=perm*util.kronecker(self.domain.getDim())
156            
157            
158         elif perm.getRank()==2:
159            perm_inv=util.inverse(perm)
160         else:
161            raise ValueError,"illegal rank of permeability."
162                    
163           which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step       self.__permeability=perm
164           PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.       self.__permeability_inv=perm_inv
165           """      
166           self.verbose=verbose           #====================
167           self.show_details= show_details and self.verbose       self.__pde_p.setValue(A=self.__permeability)
168           self.__pde_v.setTolerance(sub_rtol)           if self.solver  == self.EVAL:
169           self.__pde_p.setTolerance(sub_rtol)                pass # no extra work required
170           ATOL=self.getTolerance()           elif self.solver  == self.POST:
171           if self.verbose: print "DarcyFlux: absolute tolerance = %e"%ATOL          k=util.kronecker(self.domain.getDim())
172           #########################################################################################################################          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           #   we solve:           elif self.solver  == self.SMOOTH:
175           #            self.__pde_v.setValue(D=self.__permeability_inv)
176           #      Q^*(I-(I+D^*D)^{-1})Q dp =  Q^* (g-u0-Qp0 - (I+D^*D)^{-1} ( D^*(f-Du0)+g-u0-Qp0) )  
177           #        if g != None:
178           #   residual is      g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
179           #      if g.isEmpty():
180           #    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)            g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
181           #      else:
182           #        with v = (I+D^*D)^{-1} (D^*f+g-Qp) including BC          if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g"
183           #      self.__g=g
184           #    we use (g - Qp, v) to represent the residual. not that        if f !=None:
185           #       f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
186           #    dr(dp)=( -Q(dp), dv) with dv = - (I+D^*D)^{-1} Q(dp)       if f.isEmpty():      
187           #            f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
188           #   while the initial residual is       else:
189           #           if f.getRank()>0: raise ValueError,"illegal rank of f."
190           #      r0=( g - Qp0, v00) with v00=(I+D^*D)^{-1} (D^*f+g-Qp0) including BC       self.__f=f
          #    
          d0=self.__g-util.tensor_mult(self.__permeability,util.grad(p0))  
          self.__pde_v.setValue(Y=d0, X=self.__f*util.kronecker(self.domain), r=u0)  
          v00=self.__pde_v.getSolution(verbose=show_details)  
          if self.verbose: print "DarcyFlux: range of initial flux = ",util.inf(v00), util.sup(v00)  
          self.__pde_v.setValue(r=Data())  
          # start CG  
          r=ArithmeticTuple(d0, v00)  
          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)  
          return r[1],p  
   
     def __Aprod_PCG(self,dp):  
           if self.show_details: print "DarcyFlux: Applying operator"  
           #  -dr(dp) = (Qdp,du) where du = (I+D^*D)^{-1} (Qdp)  
           mQdp=util.tensor_mult(self.__permeability,util.grad(dp))  
           self.__pde_v.setValue(Y=mQdp,X=Data(), r=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)  
191    
192  class StokesProblemCartesian(HomogeneousSaddlePointProblem):     def getSolverOptionsFlux(self):
193          """
194          Returns the solver options used to solve the flux problems
195          :return: `SolverOptions`
196          """
197          if self.__pde_v == None:
198              return None
199          else:
200              return self.__pde_v.getSolverOptions()
201          
202       def setSolverOptionsFlux(self, options=None):
203          """
204          Sets the solver options used to solve the flux problems
205          If ``options`` is not present, the options are reset to default
206          :param options: `SolverOptions`
207          """
208          if not self.__pde_v == None:
209              self.__pde_v.setSolverOptions(options)
210        
211       def getSolverOptionsPressure(self):
212          """
213          Returns the solver options used to solve the pressure problems
214          :return: `SolverOptions`
215          """
216          return self.__pde_p.getSolverOptions()
217          
218       def setSolverOptionsPressure(self, options=None):
219          """
220          Sets the solver options used to solve the pressure problems
221          If ``options`` is not present, the options are reset to default
222          
223          :param options: `SolverOptions`
224          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
225        """        """
226        solves        return self.__pde_p.setSolverOptions(options)
227          
228       def solve(self, u0, p0):
229          """
230          solves the problem.
231          
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          :type u0: vector value on the domain (e.g. `escript.Data`).
234          :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          :type p0: scalar value on the domain (e.g. `escript.Data`).
236          :return: flux and pressure
237          :rtype: ``tuple`` of `escript.Data`.
238    
239          """
240          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               self.__pde_v.setValue(r=u0/self.perm_scale)
269                u= self.__pde_v.getSolution() * self.perm_scale
270        return u
271          
272    class StokesProblemCartesian(HomogeneousSaddlePointProblem):
273         """
274         solves
275    
276            -(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}
277                  u_{i,i}=0                  u_{i,i}=0
# Line 333  class StokesProblemCartesian(Homogeneous Line 279  class StokesProblemCartesian(Homogeneous
279            u=0 where  fixed_u_mask>0            u=0 where  fixed_u_mask>0
280            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
281    
282        if surface_stress is not given 0 is assumed.       if surface_stress is not given 0 is assumed.
283    
284        typical usage:       typical usage:
285    
286              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
287              sp.setTolerance()              sp.setTolerance()
288              sp.initialize(...)              sp.initialize(...)
289              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
290        """              sp.setStokesEquation(...) # new values for some parameters
291        def __init__(self,domain,**kwargs):              v1,p1=sp.solve(v,p)
292         """
293         def __init__(self,domain,**kwargs):
294           """           """
295           initialize the Stokes Problem           initialize the Stokes Problem
296    
297           @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
298           @type domain: L{Domain}           LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
299           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.           with macro elements for the pressure.
300    
301             :param domain: domain of the problem.
302             :type domain: `Domain`
303           """           """
304           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
305           self.domain=domain           self.domain=domain
306           self.vol=util.integrate(1.,Function(self.domain))           self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
307           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_v.setSymmetryOn()
308           self.__pde_u.setSymmetryOn()      
          # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)  
          # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)  
               
309           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
310           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
          # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)  
311           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
312    
313           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
314           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
315         self.__pde_proj.setValue(D=1)
316           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
          self.__pde_proj.setValue(D=1.)  
317    
318        def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):       def getSolverOptionsVelocity(self):
319             """
320         returns the solver options used  solve the equation for velocity.
321        
322         :rtype: `SolverOptions`
323         """
324         return self.__pde_v.getSolverOptions()
325         def setSolverOptionsVelocity(self, options=None):
326             """
327         set the solver options for solving the equation for velocity.
328        
329         :param options: new solver  options
330         :type options: `SolverOptions`
331         """
332             self.__pde_v.setSolverOptions(options)
333         def getSolverOptionsPressure(self):
334             """
335         returns the solver options used  solve the equation for pressure.
336         :rtype: `SolverOptions`
337         """
338         return self.__pde_prec.getSolverOptions()
339         def setSolverOptionsPressure(self, options=None):
340             """
341         set the solver options for solving the equation for pressure.
342         :param options: new solver  options
343         :type options: `SolverOptions`
344         """
345         self.__pde_prec.setSolverOptions(options)
346    
347         def setSolverOptionsDiv(self, options=None):
348             """
349         set the solver options for solving the equation to project the divergence of
350         the velocity onto the function space of presure.
351        
352         :param options: new solver options
353         :type options: `SolverOptions`
354         """
355         self.__pde_proj.setSolverOptions(options)
356         def getSolverOptionsDiv(self):
357             """
358         returns the solver options for solving the equation to project the divergence of
359         the velocity onto the function space of presure.
360        
361         :rtype: `SolverOptions`
362         """
363         return self.__pde_proj.getSolverOptions()
364    
365         def updateStokesEquation(self, v, p):
366             """
367             updates the Stokes equation to consider dependencies from ``v`` and ``p``
368             :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values to model parameters.
369             """
370             pass
371         def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
372            """
373            assigns new values to the model parameters.
374    
375            :param f: external force
376            :type f: `Vector` object in `FunctionSpace` `Function` or similar
377            :param fixed_u_mask: mask of locations with fixed velocity.
378            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
379            :param eta: viscosity
380            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
381            :param surface_stress: normal surface stress
382            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
383            :param stress: initial stress
384        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
385            """
386            if eta !=None:
387                k=util.kronecker(self.domain.getDim())
388                kk=util.outer(k,k)
389                self.eta=util.interpolate(eta, escript.Function(self.domain))
390            self.__pde_prec.setValue(D=1/self.eta)
391                self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
392            if restoration_factor!=None:
393                n=self.domain.getNormal()
394                self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
395            if fixed_u_mask!=None:
396                self.__pde_v.setValue(q=fixed_u_mask)
397            if f!=None: self.__f=f
398            if surface_stress!=None: self.__surface_stress=surface_stress
399            if stress!=None: self.__stress=stress
400    
401         def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
402          """          """
403          assigns values to the model parameters          assigns values to the model parameters
404    
405          @param f: external force          :param f: external force
406          @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar          :type f: `Vector` object in `FunctionSpace` `Function` or similar
407          @param fixed_u_mask: mask of locations with fixed velocity.          :param fixed_u_mask: mask of locations with fixed velocity.
408          @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
409          @param eta: viscosity          :param eta: viscosity
410          @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar          :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
411          @param surface_stress: normal surface stress          :param surface_stress: normal surface stress
412          @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
413          @param stress: initial stress          :param stress: initial stress
414      @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.  
   
415          """          """
416          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  
417    
418        def B(self,v):       def Bv(self,v,tol):
419          """           """
420          returns div(v)           returns inner product of element p and div(v)
         @rtype: equal to the type of p  
421    
422          @note: boundary conditions on p should be zero!           :param v: a residual
423          """           :return: inner product of element p and div(v)
424          if self.show_details: print "apply divergence:"           :rtype: ``float``
425          self.__pde_proj.setValue(Y=-util.div(v))           """
426          self.__pde_proj.setTolerance(self.getSubProblemTolerance())           self.__pde_proj.setValue(Y=-util.div(v))
427          return self.__pde_proj.getSolution(verbose=self.show_details)       self.getSolverOptionsDiv().setTolerance(tol)
428         self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
429             out=self.__pde_proj.getSolution()
430             return out
431    
432        def inner_pBv(self,p,Bv):       def inner_pBv(self,p,Bv):
433           """           """
434           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}  
435    
436           @rtype: equal to the type of p           :param p: a pressure increment
437             :param Bv: a residual
438             :return: inner product of element p and Bv=-div(v)
439             :rtype: ``float``
440           """           """
441           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)  
442    
443        def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
444           """           """
445           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}  
446    
447           @rtype: equal to the type of p           :param p0: a pressure
448             :param p1: a pressure
449             :return: inner product of p0 and p1
450             :rtype: ``float``
451           """           """
452           s0=util.interpolate(p0/self.eta,Function(self.domain))           s0=util.interpolate(p0, escript.Function(self.domain))
453           s1=util.interpolate(p1/self.eta,Function(self.domain))           s1=util.interpolate(p1, escript.Function(self.domain))
454           return util.integrate(s0*s1)           return util.integrate(s0*s1)
455    
456        def inner_v(self,v0,v1):       def norm_v(self,v):
457           """           """
458           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}  
459    
460           @rtype: equal to the type of v           :param v: a velovity
461             :return: norm of v
462             :rtype: non-negative ``float``
463           """           """
464       gv0=util.grad(v0)           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
465       gv1=util.grad(v1)  
          return util.integrate(util.inner(gv0,gv1))  
466    
467        def solve_A(self,u,p):       def getDV(self, p, v, tol):
468           """           """
469           solves Av=f-Au-B^*p (v=0 on fixed_u_mask)           return the value for v for a given p
470    
471             :param p: a pressure
472             :param v: a initial guess for the value v to return.
473             :return: dv given as *Adv=(f-Av-B^*p)*
474           """           """
475           if self.show_details: print "solve for velocity:"           self.updateStokesEquation(v,p)
476           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
477         self.getSolverOptionsVelocity().setTolerance(tol)
478         self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
479           if self.__stress.isEmpty():           if self.__stress.isEmpty():
480              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)))
481           else:           else:
482              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)))
483           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_v.getSolution()
484           return  out           return  out
485    
486        def solve_prec(self,p):       def norm_Bv(self,Bv):
487           if self.show_details: print "apply preconditioner:"          """
488           self.__pde_prec.setTolerance(self.getSubProblemTolerance())          Returns Bv (overwrite).
489           self.__pde_prec.setValue(Y=p)  
490           q=self.__pde_prec.getSolution(verbose=self.show_details)          :rtype: equal to the type of p
491           return q          :note: boundary conditions on p should be zero!
492            """
493            return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
494    
495         def solve_AinvBt(self,p, tol):
496             """
497             Solves *Av=B^*p* with accuracy `tol`
498    
499             :param p: a pressure increment
500             :return: the solution of *Av=B^*p*
501             :note: boundary conditions on v should be zero!
502             """
503             self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
504             out=self.__pde_v.getSolution()
505             return  out
506    
507         def solve_prec(self,Bv, tol):
508             """
509             applies preconditioner for for *BA^{-1}B^** to *Bv*
510             with accuracy `self.getSubProblemTolerance()`
511    
512             :param Bv: velocity increment
513             :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
514             :note: boundary conditions on p are zero.
515             """
516             self.__pde_prec.setValue(Y=Bv)
517         self.getSolverOptionsPressure().setTolerance(tol)
518         self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
519             out=self.__pde_prec.getSolution()
520             return out

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

  ViewVC Help
Powered by ViewVC 1.1.26