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

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

  ViewVC Help
Powered by ViewVC 1.1.26