/[escript]/trunk/escriptcore/py_src/flows.py
ViewVC logotype

Diff of /trunk/escriptcore/py_src/flows.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

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

  ViewVC Help
Powered by ViewVC 1.1.26