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

Legend:
Removed from v.2349  
changed lines
  Added in v.3905

  ViewVC Help
Powered by ViewVC 1.1.26