/[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 2288 by gross, Tue Feb 24 06:11:48 2009 UTC revision 3885 by gross, Wed Apr 4 22:12:26 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, GMRES  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.__l=util.longestEdge(self.domain)**2        """
60          self.__pde_v=LinearPDESystem(domain)        initializes the Darcy flux problem
61          if useReduced: self.__pde_v.setReducedOrderOn()        :param domain: domain of the problem
62          self.__pde_v.setSymmetryOn()        :type domain: `Domain`
63          self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))        :param useReduced: uses reduced oreder on flux and pressure
64          self.__pde_p=LinearSinglePDE(domain)        :type useReduced: ``bool``
65          self.__pde_p.setSymmetryOn()        :param solver: solver method
66          if useReduced: self.__pde_p.setReducedOrderOn()        :type solver: in [`DarcyFlow.EVAL`, `DarcyFlow.POST',  `DarcyFlow.SMOOTH' ]
67          self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))        :param verbose: if ``True`` some information on the iteration progress are printed.
68          self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))        :type verbose: ``bool``
69          self.setTolerance()        :param w: weighting factor for `DarcyFlow.POST` solver
70          self.setAbsoluteTolerance()        :type w: ``float``
71          self.setSubProblemTolerance()        
72          """
73      def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):        if not solver in [DarcyFlow.EVAL, DarcyFlow.POST,  DarcyFlow.SMOOTH ] :
74          """            raise ValueError("unknown solver %d."%solver)
75          assigns values to model parameters  
76          self.domain=domain
77          @param f: volumetic sources/sinks        self.solver=solver
78          @type f: scalar value on the domain (e.g. L{Data})        self.useReduced=useReduced
79          @param g: flux sources/sinks        self.verbose=verbose
80          @type g: vector values on the domain (e.g. L{Data})        self.l=None
81          @param location_of_fixed_pressure: mask for locations where pressure is fixed        self.w=None
82          @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})      
83          @param location_of_fixed_flux:  mask for locations where flux is fixed.        self.__pde_p=LinearSinglePDE(domain)
84          @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})        self.__pde_p.setSymmetryOn()
85          @param permeability: permeability tensor. If scalar C{s} is given the tensor with        if self.useReduced: self.__pde_p.setReducedOrderOn()
86                               C{s} on the main diagonal is used. If vector C{v} is given the tensor with  
87                               C{v} on the main diagonal is used.        if self.solver  == self.EVAL:
88          @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})           self.__pde_v=None
89             if self.verbose: print("DarcyFlow: simple solver is used.")
90          @note: the values of parameters which are not set by calling C{setValue} are not altered.  
91          @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)        elif self.solver  == self.POST:
92                 or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal           if util.inf(w)<0.:
93                 is along the M{x_i} axis.              raise ValueError("Weighting factor must be non-negative.")
94          """           if self.verbose: print("DarcyFlow: global postprocessing of flux is used.")
95          if f !=None:           self.__pde_v=LinearPDESystem(domain)
96             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))           self.__pde_v.setSymmetryOn()
97             if f.isEmpty():           if self.useReduced: self.__pde_v.setReducedOrderOn()
98                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))           self.w=w
99             else:           self.l=util.vol(self.domain)**(1./self.domain.getDim()) # length scale
100                 if f.getRank()>0: raise ValueError,"illegal rank of f."  
101             self.__f=f        elif self.solver  == self.SMOOTH:
102          if g !=None:           self.__pde_v=LinearPDESystem(domain)
103             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))           self.__pde_v.setSymmetryOn()
104             if g.isEmpty():           if self.useReduced: self.__pde_v.setReducedOrderOn()
105               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))           if self.verbose: print("DarcyFlow: flux smoothing is used.")
106             else:           self.w=0
107               if not g.getShape()==(self.domain.getDim(),):  
108                 raise ValueError,"illegal shape of g"        self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
109             self.__g=g        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          if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)        self.location_of_fixed_flux = escript.Vector(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
112          if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)        self.perm_scale=1.
113        
114          if permeability!=None:          
115             perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))     def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
116             if perm.getRank()==0:        """
117                 perm=perm*util.kronecker(self.domain.getDim())        assigns values to model parameters
118             elif perm.getRank()==1:  
119                 perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm        :param f: volumetic sources/sinks
120                 for i in range(self.domain.getDim()): perm[i,i]=perm2[i]        :type f: scalar value on the domain (e.g. `escript.Data`)
121             elif perm.getRank()==2:        :param g: flux sources/sinks
122                pass        :type g: vector values on the domain (e.g. `escript.Data`)
123             else:        :param location_of_fixed_pressure: mask for locations where pressure is fixed
124                raise ValueError,"illegal rank of permeability."        :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
125             self.__permeability=perm        :param location_of_fixed_flux:  mask for locations where flux is fixed.
126             self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))        :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      def setTolerance(self,rtol=1e-4):        :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
129          """  
130          sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if        :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          M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }               (``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          where C{atol} is an absolut tolerance (see L{setAbsoluteTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.               is along the *x_i* axis.
135    
136          @param rtol: relative tolerance for the pressure        """
137          @type rtol: non-negative C{float}        if location_of_fixed_pressure!=None:
138          """             self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)
139          if rtol<0:             self.__pde_p.setValue(q=self.location_of_fixed_pressure)
140              raise ValueError,"Relative tolerance needs to be non-negative."        if location_of_fixed_flux!=None:
141          self.__rtol=rtol            self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
142      def getTolerance(self):            if not self.__pde_v == None: self.__pde_v.setValue(q=self.location_of_fixed_flux)
143          """              
144          returns the relative tolerance        if permeability!=None:
145        
146          @return: current relative tolerance           perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
147          @rtype: C{float}           self.perm_scale=util.Lsup(util.length(perm))
148          """           if self.verbose: print(("DarcyFlow: permeability scaling factor = %e."%self.perm_scale))
149          return self.__rtol           perm=perm*(1./self.perm_scale)
150            
151      def setAbsoluteTolerance(self,atol=0.):           if perm.getRank()==0:
152          """  
153          sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if              perm_inv=(1./perm)
154                perm_inv=perm_inv*util.kronecker(self.domain.getDim())
155          M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }              perm=perm*util.kronecker(self.domain.getDim())
156            
157          where C{rtol} is an absolut tolerance (see L{setTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.          
158             elif perm.getRank()==2:
159          @param atol: absolute tolerance for the pressure              perm_inv=util.inverse(perm)
         @type atol: non-negative C{float}  
         """  
         if atol<0:  
             raise ValueError,"Absolute tolerance needs to be non-negative."  
         self.__atol=atol  
     def getAbsoluteTolerance(self):  
        """  
        returns the absolute tolerance  
         
        @return: current absolute tolerance  
        @rtype: C{float}  
        """  
        return self.__atol  
   
     def setSubProblemTolerance(self,rtol=None):  
          """  
          Sets the relative tolerance to solve the subproblem(s). If C{rtol} is not present  
          C{self.getTolerance()**2} is used.  
   
          @param rtol: relative tolerence  
          @type rtol: positive C{float}  
          """  
          if rtol == None:  
               if self.getTolerance()<=0.:  
                   raise ValueError,"A positive relative tolerance must be set."  
               self.__sub_tol=max(util.EPSILON**(0.75),self.getTolerance()**2)  
160           else:           else:
161               if rtol<=0:              raise ValueError("illegal rank of permeability.")
162                   raise ValueError,"sub-problem tolerance must be positive."          
163               self.__sub_tol=max(util.EPSILON**(0.75),rtol)           self.__permeability=perm
164             self.__permeability_inv=perm_inv
165      def getSubProblemTolerance(self):      
166           """           #====================
167           Returns the subproblem reduction factor.           self.__pde_p.setValue(A=self.__permeability)
168             if self.solver  == self.EVAL:
169           @return: subproblem reduction factor                pass # no extra work required
170           @rtype: C{float}           elif self.solver  == self.POST:
171           """                k=util.kronecker(self.domain.getDim())
172           return self.__sub_tol                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      def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, max_num_corrections=10):                self.__pde_v.setValue(D=self.__permeability_inv, A_reduced=self.omega*util.outer(k,k))
175           """           elif self.solver  == self.SMOOTH:
176           solves the problem.              self.__pde_v.setValue(D=self.__permeability_inv)
177    
178           The iteration is terminated if the residual norm is less then self.getTolerance().        if g != None:
179            g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
180           @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.          if g.isEmpty():
181           @type u0: vector value on the domain (e.g. L{Data}).               g=Vector(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
          @param p0: initial guess for the pressure. At locations in the domain marked by C{location_of_fixed_pressure} the value of C{p0} is kept unchanged.  
          @type p0: scalar value on the domain (e.g. L{Data}).  
          @param verbose: if set some information on iteration progress are printed  
          @type verbose: C{bool}  
          @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  
   
          M{Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))} with p=p0  on location_of_fixed_pressure  
   
          which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step  
          PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.  
          """  
          self.verbose=verbose or True  
          self.show_details= show_details and self.verbose  
          rtol=self.getTolerance()  
          atol=self.getAbsoluteTolerance()  
          if self.verbose: print "DarcyFlux: initial sub tolerance = %e"%self.getSubProblemTolerance()  
   
          num_corrections=0  
          converged=False  
          p=p0  
          norm_r=None  
          while not converged:  
                v=self.getFlux(p, fixed_flux=u0, show_details=self.show_details)  
                Qp=self.__Q(p)  
                norm_v=self.__L2(v)  
                norm_Qp=self.__L2(Qp)  
                if norm_v == 0.:  
                   if norm_Qp == 0.:  
                      return v,p  
                   else:  
                     fac=norm_Qp  
                else:  
                   if norm_Qp == 0.:  
                     fac=norm_v  
                   else:  
                     fac=2./(1./norm_v+1./norm_Qp)  
                ATOL=(atol+rtol*fac)  
                if self.verbose:  
                     print "DarcyFlux: L2 norm of v = %e."%norm_v  
                     print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp  
                     print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL  
                if norm_r == None or norm_r>ATOL:  
                    if num_corrections>max_num_corrections:  
                          raise ValueError,"maximum number of correction steps reached."  
                    p,r, norm_r=PCG(self.__g-util.interpolate(v,Function(self.domain))-Qp,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.1*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)  
                    num_corrections+=1  
                else:  
                    converged=True  
          return v,p  
 #  
 #                
 #               r_hat=g-util.interpolate(v,Function(self.domain))-Qp  
 #               #===========================================================================  
 #               norm_r_hat=self.__L2(r_hat)  
 #               norm_v=self.__L2(v)  
 #               norm_g=self.__L2(g)  
 #               norm_gv=self.__L2(g-v)  
 #               norm_Qp=self.__L2(Qp)  
 #               norm_gQp=self.__L2(g-Qp)  
 #               fac=min(max(norm_v,norm_gQp),max(norm_Qp,norm_gv))  
 #               fac=min(norm_v,norm_Qp,norm_gv)  
 #               norm_r_hat_PCG=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))  
 #               print "norm_r_hat = ",norm_r_hat,norm_r_hat_PCG, norm_r_hat_PCG/norm_r_hat  
 #               if r!=None:  
 #                   print "diff = ",self.__L2(r-r_hat)/norm_r_hat  
 #                   sub_tol=min(rtol/self.__L2(r-r_hat)*norm_r_hat,1.)*self.getSubProblemTolerance()  
 #                   self.setSubProblemTolerance(sub_tol)  
 #                   print "subtol_new=",self.getSubProblemTolerance()  
 #               print "norm_v = ",norm_v  
 #               print "norm_gv = ",norm_gv  
 #               print "norm_Qp = ",norm_Qp  
 #               print "norm_gQp = ",norm_gQp  
 #               print "norm_g = ",norm_g  
 #               print "max(norm_v,norm_gQp)=",max(norm_v,norm_gQp)  
 #               print "max(norm_Qp,norm_gv)=",max(norm_Qp,norm_gv)  
 #               if fac == 0:  
 #                   if self.verbose: print "DarcyFlux: trivial case!"  
 #                   return v,p  
 #               #===============================================================================  
 #               # norm_v=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(v),v))  
 #               # norm_Qp=self.__L2(Qp)  
 #               norm_r_hat=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))  
 #               # print "**** norm_v, norm_Qp :",norm_v,norm_Qp  
 #  
 #               ATOL=(atol+rtol*2./(1./norm_v+1./norm_Qp))  
 #               if self.verbose:  
 #                   print "DarcyFlux: residual = %e"%norm_r_hat  
 #                   print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL  
 #               if norm_r_hat <= ATOL:  
 #                   print "DarcyFlux: iteration finalized."  
 #                   converged=True  
 #               else:  
 #                   # p=GMRES(r_hat,self.__Aprod, p, self.__inner_GMRES, atol=ATOL, rtol=0., iter_max=max_iter, iter_restart=20, verbose=self.verbose,P_R=self.__Msolve_PCG)  
 #                   # p,r=PCG(r_hat,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL*min(0.1,norm_r_hat_PCG/norm_r_hat), rtol=0.,iter_max=max_iter, verbose=self.verbose)  
 #                   p,r, norm_r=PCG(r_hat,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.1*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)  
 #               print "norm_r =",norm_r  
 #         return v,p  
     def __L2(self,v):  
          return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))  
   
     def __Q(self,p):  
           return util.tensor_mult(self.__permeability,util.grad(p))  
   
     def __Aprod(self,dp):  
           self.__pde_v.setTolerance(self.getSubProblemTolerance())  
           if self.show_details: print "DarcyFlux: Applying operator"  
           Qdp=self.__Q(dp)  
           self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())  
           du=self.__pde_v.getSolution(verbose=self.show_details)  
           return Qdp+du  
     def __inner_GMRES(self,r,s):  
          return util.integrate(util.inner(r,s))  
   
     def __inner_PCG(self,p,r):  
          return util.integrate(util.inner(self.__Q(p), r))  
   
     def __Msolve_PCG(self,r):  
           self.__pde_p.setTolerance(self.getSubProblemTolerance())  
           if self.show_details: print "DarcyFlux: Applying preconditioner"  
           self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())  
           return self.__pde_p.getSolution(verbose=self.show_details)  
   
     def getFlux(self,p=None, fixed_flux=Data(), 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(self.getSubProblemTolerance())  
         g=self.__g  
         f=self.__f  
         self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)  
         if p == None:  
            self.__pde_v.setValue(Y=g)  
182          else:          else:
183             self.__pde_v.setValue(Y=g-self.__Q(p))               if not g.getShape()==(self.domain.getDim(),): raise ValueError("illegal shape of g")
184          return self.__pde_v.getSolution(verbose=show_details)          self.__g=g
185          if f !=None:
186             f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
187             if f.isEmpty():      
188                 f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
189             else:
190                 if f.getRank()>0: raise ValueError("illegal rank of f.")
191             self.__f=f
192    
193       def getSolverOptionsFlux(self):
194          """
195          Returns the solver options used to solve the flux problems
196          :return: `SolverOptions`
197          """
198          if self.__pde_v == None:
199              return None
200          else:
201              return self.__pde_v.getSolverOptions()
202          
203       def setSolverOptionsFlux(self, options=None):
204          """
205          Sets the solver options used to solve the flux problems
206          If ``options`` is not present, the options are reset to default
207          :param options: `SolverOptions`
208          """
209          if not self.__pde_v == None:
210              self.__pde_v.setSolverOptions(options)
211        
212       def getSolverOptionsPressure(self):
213          """
214          Returns the solver options used to solve the pressure problems
215          :return: `SolverOptions`
216          """
217          return self.__pde_p.getSolverOptions()
218          
219       def setSolverOptionsPressure(self, options=None):
220          """
221          Sets the solver options used to solve the pressure problems
222          If ``options`` is not present, the options are reset to default
223          
224          :param options: `SolverOptions`
225          :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
226          """
227          return self.__pde_p.setSolverOptions(options)
228          
229       def solve(self, u0, p0):
230          """
231          solves the problem.
232          
233          :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.
234          :type u0: vector value on the domain (e.g. `escript.Data`).
235          :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.
236          :type p0: scalar value on the domain (e.g. `escript.Data`).
237          :return: flux and pressure
238          :rtype: ``tuple`` of `escript.Data`.
239    
240          """
241          self.__pde_p.setValue(X=self.__g * 1./self.perm_scale,
242                                Y=self.__f * 1./self.perm_scale,
243                                y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux * 1./self.perm_scale ),
244                                r=p0)
245          p=self.__pde_p.getSolution()
246          u = self.getFlux(p, u0)
247          return u,p
248          
249       def getFlux(self,p, u0=None):
250            """
251            returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
252            on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
253            Notice that ``g`` is used, see `setValue`.
254    
255            :param p: pressure.
256            :type p: scalar value on the domain (e.g. `escript.Data`).
257            :param u0: flux on the locations of the domain marked be ``location_of_fixed_flux``.
258            :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
259            :return: flux
260            :rtype: `escript.Data`
261            """
262            if self.solver  == self.EVAL:
263               u = self.__g-self.perm_scale * util.tensor_mult(self.__permeability,util.grad(p))
264            elif self.solver  == self.POST or self.solver  == self.SMOOTH:
265                self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g * 1./self.perm_scale)-util.grad(p))
266                if u0 == None:
267                   self.__pde_v.setValue(r=escript.Data())
268                else:
269                   if not isinstance(u0, escript.Data) : u0 = escript.Vector(u0, escript.Solution(self.domain))
270                   self.__pde_v.setValue(r=1./self.perm_scale * u0)
271                   u= self.__pde_v.getSolution() * self.perm_scale
272            return u
273          
274  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
275       """       """
276       solves       solves
# Line 386  class StokesProblemCartesian(Homogeneous Line 289  class StokesProblemCartesian(Homogeneous
289              sp.setTolerance()              sp.setTolerance()
290              sp.initialize(...)              sp.initialize(...)
291              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
292                sp.setStokesEquation(...) # new values for some parameters
293                v1,p1=sp.solve(v,p)
294       """       """
295       def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
296           """           """
297           initialize the Stokes Problem           initialize the Stokes Problem
298    
299           @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
300           @type domain: L{Domain}           LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
301           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.           with macro elements for the pressure.
302    
303             :param domain: domain of the problem.
304             :type domain: `Domain`
305           """           """
306           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
307           self.domain=domain           self.domain=domain
308           self.vol=util.integrate(1.,Function(self.domain))           self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
309           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_v.setSymmetryOn()
310           self.__pde_u.setSymmetryOn()      
          # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)  
          # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)  
   
311           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
312           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
          # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)  
313           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
314    
315       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):           self.__pde_proj=LinearPDE(domain)
316             self.__pde_proj.setReducedOrderOn()
317             self.__pde_proj.setValue(D=1)
318             self.__pde_proj.setSymmetryOn()
319    
320         def getSolverOptionsVelocity(self):
321             """
322         returns the solver options used  solve the equation for velocity.
323        
324         :rtype: `SolverOptions`
325         """
326             return self.__pde_v.getSolverOptions()
327         def setSolverOptionsVelocity(self, options=None):
328             """
329         set the solver options for solving the equation for velocity.
330        
331         :param options: new solver  options
332         :type options: `SolverOptions`
333         """
334             self.__pde_v.setSolverOptions(options)
335         def getSolverOptionsPressure(self):
336             """
337         returns the solver options used  solve the equation for pressure.
338         :rtype: `SolverOptions`
339         """
340             return self.__pde_prec.getSolverOptions()
341         def setSolverOptionsPressure(self, options=None):
342             """
343         set the solver options for solving the equation for pressure.
344         :param options: new solver  options
345         :type options: `SolverOptions`
346         """
347             self.__pde_prec.setSolverOptions(options)
348    
349         def setSolverOptionsDiv(self, options=None):
350             """
351         set the solver options for solving the equation to project the divergence of
352         the velocity onto the function space of presure.
353        
354         :param options: new solver options
355         :type options: `SolverOptions`
356         """
357             self.__pde_proj.setSolverOptions(options)
358         def getSolverOptionsDiv(self):
359             """
360         returns the solver options for solving the equation to project the divergence of
361         the velocity onto the function space of presure.
362        
363         :rtype: `SolverOptions`
364         """
365             return self.__pde_proj.getSolverOptions()
366    
367         def updateStokesEquation(self, v, p):
368             """
369             updates the Stokes equation to consider dependencies from ``v`` and ``p``
370             :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values to model parameters.
371             """
372             pass
373         def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
374            """
375            assigns new values to the model parameters.
376    
377            :param f: external force
378            :type f: `Vector` object in `FunctionSpace` `Function` or similar
379            :param fixed_u_mask: mask of locations with fixed velocity.
380            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
381            :param eta: viscosity
382            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
383            :param surface_stress: normal surface stress
384            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
385            :param stress: initial stress
386        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
387            """
388            if eta !=None:
389                k=util.kronecker(self.domain.getDim())
390                kk=util.outer(k,k)
391                self.eta=util.interpolate(eta, escript.Function(self.domain))
392                self.__pde_prec.setValue(D=1/self.eta)
393                self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
394            if restoration_factor!=None:
395                n=self.domain.getNormal()
396                self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
397            if fixed_u_mask!=None:
398                self.__pde_v.setValue(q=fixed_u_mask)
399            if f!=None: self.__f=f
400            if surface_stress!=None: self.__surface_stress=surface_stress
401            if stress!=None: self.__stress=stress
402    
403         def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
404          """          """
405          assigns values to the model parameters          assigns values to the model parameters
406    
407          @param f: external force          :param f: external force
408          @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar          :type f: `Vector` object in `FunctionSpace` `Function` or similar
409          @param fixed_u_mask: mask of locations with fixed velocity.          :param fixed_u_mask: mask of locations with fixed velocity.
410          @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
411          @param eta: viscosity          :param eta: viscosity
412          @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar          :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
413          @param surface_stress: normal surface stress          :param surface_stress: normal surface stress
414          @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
415          @param stress: initial stress          :param stress: initial stress
416      @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.  
   
417          """          """
418          self.eta=eta          self.setStokesEquation(f,fixed_u_mask, eta, surface_stress, stress, restoration_factor)
         A =self.__pde_u.createCoefficient("A")  
     self.__pde_u.setValue(A=Data())  
         for i in range(self.domain.getDim()):  
         for j in range(self.domain.getDim()):  
             A[i,j,j,i] += 1.  
             A[i,j,i,j] += 1.  
     self.__pde_prec.setValue(D=1/self.eta)  
         self.__pde_u.setValue(A=A*self.eta,q=fixed_u_mask)  
         self.__f=f  
         self.__surface_stress=surface_stress  
         self.__stress=stress  
419    
420       def inner_pBv(self,p,v):       def Bv(self,v,tol):
421           """           """
422           returns inner product of element p and div(v)           returns inner product of element p and div(v)
423    
424           @param p: a pressure increment           :param v: a residual
425           @param v: a residual           :return: inner product of element p and div(v)
426           @return: inner product of element p and div(v)           :rtype: ``float``
427           @rtype: C{float}           """
428             self.__pde_proj.setValue(Y=-util.div(v))
429             self.getSolverOptionsDiv().setTolerance(tol)
430             self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
431             out=self.__pde_proj.getSolution()
432             return out
433    
434         def inner_pBv(self,p,Bv):
435             """
436             returns inner product of element p and Bv=-div(v)
437    
438             :param p: a pressure increment
439             :param Bv: a residual
440             :return: inner product of element p and Bv=-div(v)
441             :rtype: ``float``
442           """           """
443           return util.integrate(-p*util.div(v))           return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
444    
445       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
446           """           """
447           Returns inner product of p0 and p1           Returns inner product of p0 and p1
448    
449           @param p0: a pressure           :param p0: a pressure
450           @param p1: a pressure           :param p1: a pressure
451           @return: inner product of p0 and p1           :return: inner product of p0 and p1
452           @rtype: C{float}           :rtype: ``float``
453           """           """
454           s0=util.interpolate(p0/self.eta,Function(self.domain))           s0=util.interpolate(p0, escript.Function(self.domain))
455           s1=util.interpolate(p1/self.eta,Function(self.domain))           s1=util.interpolate(p1, escript.Function(self.domain))
456           return util.integrate(s0*s1)           return util.integrate(s0*s1)
457    
458       def norm_v(self,v):       def norm_v(self,v):
459           """           """
460           returns the norm of v           returns the norm of v
461    
462           @param v: a velovity           :param v: a velovity
463           @return: norm of v           :return: norm of v
464           @rtype: non-negative C{float}           :rtype: non-negative ``float``
465           """           """
466           return util.sqrt(util.integrate(util.length(util.grad(v))))           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
467    
468       def getV(self, p, v0):  
469         def getDV(self, p, v, tol):
470           """           """
471           return the value for v for a given p (overwrite)           return the value for v for a given p
472    
473           @param p: a pressure           :param p: a pressure
474           @param v0: a initial guess for the value v to return.           :param v: a initial guess for the value v to return.
475           @return: v given as M{v= A^{-1} (f-B^*p)}           :return: dv given as *Adv=(f-Av-B^*p)*
476           """           """
477           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.updateStokesEquation(v,p)
478           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
479             self.getSolverOptionsVelocity().setTolerance(tol)
480             self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
481           if self.__stress.isEmpty():           if self.__stress.isEmpty():
482              self.__pde_u.setValue(X=p*util.kronecker(self.domain))              self.__pde_v.setValue(X=p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
483           else:           else:
484              self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))              self.__pde_v.setValue(X=self.__stress+p*util.kronecker(self.domain)-2*self.eta*util.symmetric(util.grad(v)))
485           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_v.getSolution()
486           return  out           return  out
487    
488         def norm_Bv(self,Bv):
          raise NotImplementedError,"no v calculation implemented."  
   
   
      def norm_Bv(self,v):  
489          """          """
490          Returns Bv (overwrite).          Returns Bv (overwrite).
491    
492          @rtype: equal to the type of p          :rtype: equal to the type of p
493          @note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
494          """          """
495          return util.sqrt(util.integrate(util.div(v)**2))          return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
496    
497       def solve_AinvBt(self,p):       def solve_AinvBt(self,p, tol):
498           """           """
499           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}           Solves *Av=B^*p* with accuracy `tol`
500    
501           @param p: a pressure increment           :param p: a pressure increment
502           @return: the solution of M{Av=B^*p}           :return: the solution of *Av=B^*p*
503           @note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
504           """           """
505           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
506           self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))           out=self.__pde_v.getSolution()
          out=self.__pde_u.getSolution(verbose=self.show_details)  
507           return  out           return  out
508    
509       def solve_precB(self,v):       def solve_prec(self,Bv, tol):
510           """           """
511           applies preconditioner for for M{BA^{-1}B^*} to M{Bv}           applies preconditioner for for *BA^{-1}B^** to *Bv*
512           with accuracy L{self.getSubProblemTolerance()} (overwrite).           with accuracy `self.getSubProblemTolerance()`
513    
514           @param v: velocity increment           :param Bv: velocity increment
515           @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}           :return: *p=P(Bv)* where *P^{-1}* is an approximation of *BA^{-1}B^ * )*
516           @note: boundary conditions on p are zero.           :note: boundary conditions on p are zero.
517           """           """
518           self.__pde_prec.setValue(Y=-util.div(v))           self.__pde_prec.setValue(Y=Bv)
519           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           self.getSolverOptionsPressure().setTolerance(tol)
520           return self.__pde_prec.getSolution(verbose=self.show_details)           self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
521             out=self.__pde_prec.getSolution()
522             return out

Legend:
Removed from v.2288  
changed lines
  Added in v.3885

  ViewVC Help
Powered by ViewVC 1.1.26