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

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

  ViewVC Help
Powered by ViewVC 1.1.26