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

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

  ViewVC Help
Powered by ViewVC 1.1.26