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

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

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

revision 2344 by jfenwick, Mon Mar 30 02:13:58 2009 UTC revision 3510 by gross, Fri May 13 06:09:49 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"""
# Line 21  __url__="https://launchpad.net/escript-f Line 22  __url__="https://launchpad.net/escript-f
22  """  """
23  Some models for flow  Some models for flow
24    
25  @var __author__: name of author  :var __author__: name of author
26  @var __copyright__: copyrights  :var __copyright__: copyrights
27  @var __license__: licence agreement  :var __license__: licence agreement
28  @var __url__: url entry point on documentation  :var __url__: url entry point on documentation
29  @var __version__: version  :var __version__: version
30  @var __date__: date of the version  :var __date__: date of the version
31  """  """
32    
33  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
34    
35  from escript import *  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 SIMPLE: simple solver
50      """     :cvar POST: solver using global postprocessing of flux
51       :cvar STAB: solver uses (non-symmetric) stabilization
52      def __init__(self, domain,useReduced=False):     :cvar SYMSTAB: solver uses symmetric stabilization
53          """     """
54          initializes the Darcy flux problem     SIMPLE="SIMPLE"
55          @param domain: domain of the problem     POST="POST"
56          @type domain: L{Domain}     STAB="STAB"
57          """     SYMSTAB="SYMSTAB"
58          self.domain=domain     def __init__(self, domain, useReduced=False, solver="SYMSTAB", 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.SIMPLE`, `DarcyFlow.POST', `DarcyFlow.STAB`, `DarcyFlow.SYMSTAB` ]
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):        self.domain=domain
74          """        self.solver=solver
75          assigns values to model parameters        self.useReduced=useReduced
76          self.verbose=verbose
77          @param f: volumetic sources/sinks        self.scale=1.
78          @type f: scalar value on the domain (e.g. L{Data})        
79          @param g: flux sources/sinks        
80          @type g: vector values on the domain (e.g. L{Data})        self.__pde_v=LinearPDESystem(domain)
81          @param location_of_fixed_pressure: mask for locations where pressure is fixed        self.__pde_v.setSymmetryOn()
82          @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})        if self.useReduced: self.__pde_v.setReducedOrderOn()
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.SIMPLE:
88          @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})       if self.verbose: print "DarcyFlow: simple solver is used."
89             self.__pde_v.setValue(D=util.kronecker(self.domain.getDim()))
90          @note: the values of parameters which are not set by calling C{setValue} are not altered.        elif self.solver  == self.POST:
91          @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)       self.w=w
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:        elif self.solver  == self.STAB:
96             f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))        if self.verbose: print "DarcyFlow: (non-symmetric) stabilization is used."
97             if f.isEmpty():        elif  self.solver  == self.SYMSTAB:
98                 f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))        if self.verbose: print "DarcyFlow: symmetric stabilization is used."
99             else:        else:
100                 if f.getRank()>0: raise ValueError,"illegal rank of f."      raise ValueError,"unknown solver %s"%self.solver
101             self.__f=f        self.__f=escript.Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("X"))
102          if g !=None:        self.__g=escript.Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
103             g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))        self.location_of_fixed_pressure = escript.Scalar(0, self.__pde_p.getFunctionSpaceForCoefficient("q"))
104             if g.isEmpty():        self.location_of_fixed_flux = escript.Vector(0, self.__pde_v.getFunctionSpaceForCoefficient("q"))
105               g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))        self.setTolerance()
106             else:      
107               if not g.getShape()==(self.domain.getDim(),):          
108                 raise ValueError,"illegal shape of g"     def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
109             self.__g=g        """
110          assigns values to model parameters
111          if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)  
112          if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)        :param f: volumetic sources/sinks
113          :type f: scalar value on the domain (e.g. `escript.Data`)
114          if permeability!=None:        :param g: flux sources/sinks
115             perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))        :type g: vector values on the domain (e.g. `escript.Data`)
116             if perm.getRank()==0:        :param location_of_fixed_pressure: mask for locations where pressure is fixed
117                 perm=perm*util.kronecker(self.domain.getDim())        :type location_of_fixed_pressure: scalar value on the domain (e.g. `escript.Data`)
118             elif perm.getRank()==1:        :param location_of_fixed_flux:  mask for locations where flux is fixed.
119                 perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm        :type location_of_fixed_flux: vector values on the domain (e.g. `escript.Data`)
120                 for i in range(self.domain.getDim()): perm[i,i]=perm2[i]        :param permeability: permeability tensor. If scalar ``s`` is given the tensor with ``s`` on the main diagonal is used.
121             elif perm.getRank()==2:        :type permeability: scalar or symmetric tensor values on the domain (e.g. `escript.Data`)
122                pass  
123             else:        :note: the values of parameters which are not set by calling ``setValue`` are not altered.
124                raise ValueError,"illegal rank of permeability."        :note: at any point on the boundary of the domain the pressure
125             self.__permeability=perm               (``location_of_fixed_pressure`` >0) or the normal component of the
126             self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))               flux (``location_of_fixed_flux[i]>0``) if direction of the normal
127                 is along the *x_i* axis.
128      def setTolerance(self,rtol=1e-4):  
129          """        """
130          sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if        if location_of_fixed_pressure!=None:
131               self.location_of_fixed_pressure=util.wherePositive(location_of_fixed_pressure)
132          M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }             self.__pde_p.setValue(q=self.location_of_fixed_pressure)
133          if location_of_fixed_flux!=None:
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}}.            self.location_of_fixed_flux=util.wherePositive(location_of_fixed_flux)
135              self.__pde_v.setValue(q=self.location_of_fixed_flux)
136          @param rtol: relative tolerance for the pressure        
137          @type rtol: non-negative C{float}              
138          """        # pressure  is rescaled by the factor 1/self.scale
139          if rtol<0:        if permeability!=None:
140              raise ValueError,"Relative tolerance needs to be non-negative."      
141          self.__rtol=rtol       perm=util.interpolate(permeability,self.__pde_v.getFunctionSpaceForCoefficient("A"))
142      def getTolerance(self):           V=util.vol(self.domain)
143          """           l=V**(1./self.domain.getDim())
144          returns the relative tolerance          
145         if perm.getRank()==0:
146          @return: current relative tolerance          perm_inv=(1./perm)
147          @rtype: C{float}              self.scale=util.integrate(perm_inv)/V*l
148          """          perm_inv=perm_inv*((1./self.scale)*util.kronecker(self.domain.getDim()))
149          return self.__rtol          perm=perm*(self.scale*util.kronecker(self.domain.getDim()))
150            
151      def setAbsoluteTolerance(self,atol=0.):          
152          """       elif perm.getRank()==2:
153          sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if          perm_inv=util.inverse(perm)
154                self.scale=util.sqrt(util.integrate(util.length(perm_inv)**2)/V)*l
155          M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }          perm_inv*=(1./self.scale)
156            perm=perm*self.scale
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}}.       else:
158            raise ValueError,"illegal rank of permeability."
159          @param atol: absolute tolerance for the pressure          
160          @type atol: non-negative C{float}       self.__permeability=perm
161          """       self.__permeability_inv=perm_inv
162          if atol<0:       if self.verbose: print "DarcyFlow: scaling factor for pressure is %e."%self.scale
163              raise ValueError,"Absolute tolerance needs to be non-negative."      
164          self.__atol=atol       if self.solver  == self.SIMPLE:
165      def getAbsoluteTolerance(self):          self.__pde_p.setValue(A=self.__permeability)
166         """       elif self.solver  == self.POST:
167         returns the absolute tolerance          self.__pde_p.setValue(A=self.__permeability)
168                  k=util.kronecker(self.domain.getDim())
169         @return: current absolute tolerance          self.lamb = self.w*util.length(perm_inv)*l
170         @rtype: C{float}          self.__pde_v.setValue(D=self.__permeability_inv, A=self.lamb*self.domain.getSize()*util.outer(k,k))
171         """       elif self.solver  == self.STAB:
172         return self.__atol          self.__pde_p.setValue(A=0.5*self.__permeability)
173            self.__pde_v.setValue(D=0.5*self.__permeability_inv)
174      def setSubProblemTolerance(self,rtol=None):       elif  self.solver  == self.SYMSTAB:
175           """          self.__pde_p.setValue(A=0.5*self.__permeability)
176           Sets the relative tolerance to solve the subproblem(s). If C{rtol} is not present          self.__pde_v.setValue(D=0.5*self.__permeability_inv)
177           C{self.getTolerance()**2} is used.  
178          if g != None:
179           @param rtol: relative tolerence      g=util.interpolate(g, self.__pde_v.getFunctionSpaceForCoefficient("Y"))
180           @type rtol: positive C{float}      if g.isEmpty():
181           """            g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
182           if rtol == None:      else:
183                if self.getTolerance()<=0.:          if not g.getShape()==(self.domain.getDim(),): raise ValueError,"illegal shape of g"
184                    raise ValueError,"A positive relative tolerance must be set."      self.__g=g
185                self.__sub_tol=max(util.EPSILON**(0.75),self.getTolerance()**2)        if f !=None:
186           else:       f=util.interpolate(f, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
187               if rtol<=0:       if f.isEmpty():      
188                   raise ValueError,"sub-problem tolerance must be positive."            f=Scalar(0,self.__pde_p.getFunctionSpaceForCoefficient("Y"))
189               self.__sub_tol=max(util.EPSILON**(0.75),rtol)       else:
190             if f.getRank()>0: raise ValueError,"illegal rank of f."
191      def getSubProblemTolerance(self):       self.__f=f
192           """     def getSolverOptionsFlux(self):
193           Returns the subproblem reduction factor.        """
194          Returns the solver options used to solve the flux problems
195           @return: subproblem reduction factor        :return: `SolverOptions`
196           @rtype: C{float}        """
197           """        return self.__pde_v.getSolverOptions()
198           return self.__sub_tol        
199       def setSolverOptionsFlux(self, options=None):
200      def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, max_num_corrections=10):        """
201           """        Sets the solver options used to solve the flux problems
202           solves the problem.        If ``options`` is not present, the options are reset to default
203          :param options: `SolverOptions`
204           The iteration is terminated if the residual norm is less then self.getTolerance().        """
205          return self.__pde_v.setSolverOptions(options)
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}).     def getSolverOptionsPressure(self):
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.        """
209           @type p0: scalar value on the domain (e.g. L{Data}).        Returns the solver options used to solve the pressure problems
210           @param verbose: if set some information on iteration progress are printed        :return: `SolverOptions`
211           @type verbose: C{bool}        """
212           @param show_details:  if set information on the subiteration process are printed.        return self.__pde_p.getSolverOptions()
213           @type show_details: C{bool}        
214           @return: flux and pressure     def setSolverOptionsPressure(self, options=None):
215           @rtype: C{tuple} of L{Data}.        """
216          Sets the solver options used to solve the pressure problems
217           @note: The problem is solved as a least squares form        If ``options`` is not present, the options are reset to default
218          
219           M{(I+D^*D)u+Qp=D^*f+g}        :param options: `SolverOptions`
220           M{Q^*u+Q^*Qp=Q^*g}        :note: if the adaption of subtolerance is choosen, the tolerance set by ``options`` will be overwritten before the solver is called.
221          """
222           where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.        return self.__pde_p.setSolverOptions(options)
223           We eliminate the flux form the problem by setting        
224       def setTolerance(self,rtol=1e-4):
225           M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux        """
226          sets the relative tolerance ``rtol`` for the pressure for the stabelized solvers.
227           form the first equation. Inserted into the second equation we get        
228          :param rtol: relative tolerance for the pressure
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        :type rtol: non-negative ``float``
230          """
231           which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step        if rtol<0:
232           PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.       raise ValueError,"Relative tolerance needs to be non-negative."
233           """        self.__rtol=rtol
234           self.verbose=verbose or True        
235           self.show_details= show_details and self.verbose     def getTolerance(self):
236           rtol=self.getTolerance()        """
237           atol=self.getAbsoluteTolerance()        returns the relative tolerance
238           if self.verbose: print "DarcyFlux: initial sub tolerance = %e"%self.getSubProblemTolerance()        :return: current relative tolerance
239          :rtype: ``float``
240           num_corrections=0        """
241           converged=False        return self.__rtol
242           p=p0        
243           norm_r=None     def solve(self,u0,p0, max_iter=100, iter_restart=20):
244           while not converged:        """
245                 v=self.getFlux(p, fixed_flux=u0, show_details=self.show_details)        solves the problem.
246                 Qp=self.__Q(p)        
247                 norm_v=self.__L2(v)        The iteration is terminated if the residual norm is less then self.getTolerance().
248                 norm_Qp=self.__L2(Qp)  
249                 if norm_v == 0.:        :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.
250                    if norm_Qp == 0.:        :type u0: vector value on the domain (e.g. `escript.Data`).
251                       return v,p        :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.
252                    else:        :type p0: scalar value on the domain (e.g. `escript.Data`).
253                      fac=norm_Qp        :param max_iter: maximum number of (outer) iteration steps for the stabilization solvers,
254                 else:        :type max_iter: ``int``
255                    if norm_Qp == 0.:        :param iter_restart: number of steps after which the iteration is restarted. The larger ``iter_restart`` the larger the required memory.
256                      fac=norm_v                             A small value for ``iter_restart`` may require a large number of iteration steps or may even lead to a failure
257                    else:                             of the iteration. ``iter_restart`` is relevant for the stabilization solvers only.
258                      fac=2./(1./norm_v+1./norm_Qp)        :type iter_restart: ``int``
259                 ATOL=(atol+rtol*fac)        :return: flux and pressure
260                 if self.verbose:        :rtype: ``tuple`` of `escript.Data`.
261                      print "DarcyFlux: L2 norm of v = %e."%norm_v  
262                      print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp        """
263                      print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL        # rescale initial guess:
264                 if norm_r == None or norm_r>ATOL:        p0=p0/self.scale
265                     if num_corrections>max_num_corrections:        if self.solver  == self.SIMPLE or self.solver  == self.POST :
266                           raise ValueError,"maximum number of correction steps reached."          self.__pde_p.setValue(X=self.__g ,
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)                                Y=self.__f,
268                     num_corrections+=1                                y= - util.inner(self.domain.getNormal(),u0 * self.location_of_fixed_flux),
269                 else:                                r=p0)
270                     converged=True          p=self.__pde_p.getSolution()
271           return v,p          u = self.getFlux(p, u0)
272  #        elif  self.solver  == self.STAB:
273  #                    u,p = self.__solve_STAB(u0,p0, max_iter, iter_restart)
274  #               r_hat=g-util.interpolate(v,Function(self.domain))-Qp        elif  self.solver  == self.SYMSTAB:
275  #               #===========================================================================      u,p = self.__solve_SYMSTAB(u0,p0, max_iter, iter_restart)
276  #               norm_r_hat=self.__L2(r_hat)      
277  #               norm_v=self.__L2(v)        if self.verbose:
278  #               norm_g=self.__L2(g)          KGp=util.tensor_mult(self.__permeability,util.grad(p))
279  #               norm_gv=self.__L2(g-v)          def_p=self.__g-(u+KGp)
280  #               norm_Qp=self.__L2(Qp)          def_v=self.__f-util.div(u, self.__pde_v.getFunctionSpaceForCoefficient("X"))
281  #               norm_gQp=self.__L2(g-Qp)          print "DarcyFlux: |g-u-K*grad(p)|_2 = %e (|u|_2 = %e)."%(self.__L2(def_p),self.__L2(u))
282  #               fac=min(max(norm_v,norm_gQp),max(norm_Qp,norm_gv))          print "DarcyFlux: |f-div(u)|_2 = %e (|grad(u)|_2 = %e)."%(self.__L2(def_v),self.__L2(util.grad(u)))
283  #               fac=min(norm_v,norm_Qp,norm_gv)        #rescale result
284  #               norm_r_hat_PCG=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))        p=p*self.scale
285  #               print "norm_r_hat = ",norm_r_hat,norm_r_hat_PCG, norm_r_hat_PCG/norm_r_hat        return u,p
286  #               if r!=None:        
287  #                   print "diff = ",self.__L2(r-r_hat)/norm_r_hat     def getFlux(self,p, u0=None):
288  #                   sub_tol=min(rtol/self.__L2(r-r_hat)*norm_r_hat,1.)*self.getSubProblemTolerance()          """
289  #                   self.setSubProblemTolerance(sub_tol)          returns the flux for a given pressure ``p`` where the flux is equal to ``u0``
290  #                   print "subtol_new=",self.getSubProblemTolerance()          on locations where ``location_of_fixed_flux`` is positive (see `setValue`).
291  #               print "norm_v = ",norm_v          Notice that ``g`` and ``f`` are used, see `setValue`.
292  #               print "norm_gv = ",norm_gv  
293  #               print "norm_Qp = ",norm_Qp          :param p: pressure.
294  #               print "norm_gQp = ",norm_gQp          :type p: scalar value on the domain (e.g. `escript.Data`).
295  #               print "norm_g = ",norm_g          :param u0: flux on the locations of the domain marked be ``location_of_fixed_flux``.
296  #               print "max(norm_v,norm_gQp)=",max(norm_v,norm_gQp)          :type u0: vector values on the domain (e.g. `escript.Data`) or ``None``
297  #               print "max(norm_Qp,norm_gv)=",max(norm_Qp,norm_gv)          :return: flux
298  #               if fac == 0:          :rtype: `escript.Data`
299  #                   if self.verbose: print "DarcyFlux: trivial case!"          """
300  #                   return v,p          if self.solver  == self.SIMPLE or self.solver  == self.POST  :
301  #               #===============================================================================              KGp=util.tensor_mult(self.__permeability,util.grad(p))
302  #               # norm_v=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(v),v))              self.__pde_v.setValue(Y=self.__g-KGp, X=escript.Data())
303  #               # norm_Qp=self.__L2(Qp)              if u0 == None:
304  #               norm_r_hat=util.sqrt(self.__inner_PCG(self.__Msolve_PCG(r_hat),r_hat))             self.__pde_v.setValue(r=escript.Data())
305  #               # print "**** norm_v, norm_Qp :",norm_v,norm_Qp          else:
306  #             self.__pde_v.setValue(r=u0)
307  #               ATOL=(atol+rtol*2./(1./norm_v+1./norm_Qp))              u= self.__pde_v.getSolution()
308  #               if self.verbose:      elif self.solver  == self.POST:
309  #                   print "DarcyFlux: residual = %e"%norm_r_hat              self.__pde_v.setValue(Y=util.tensor_mult(self.__permeability_inv,self.__g)-util.grad(p),
310  #                   print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL                                    X=self.lamb * self.__f * util.kronecker(self.domain.getDim()))
311  #               if norm_r_hat <= ATOL:              if u0 == None:
312  #                   print "DarcyFlux: iteration finalized."             self.__pde_v.setValue(r=escript.Data())
313  #                   converged=True          else:
314  #               else:             self.__pde_v.setValue(r=u0)
315  #                   # 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)              u= self.__pde_v.getSolution()
316  #                   # 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)      elif self.solver  == self.STAB:
317  #                   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)           gp=util.grad(p)
318  #               print "norm_r =",norm_r           self.__pde_v.setValue(Y=0.5*(util.tensor_mult(self.__permeability_inv,self.__g)+gp),
319  #         return v,p                                 X= p * util.kronecker(self.domain.getDim()),
320      def __L2(self,v):                                 y= - p * self.domain.getNormal())                          
321           return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))           if u0 == None:
322               self.__pde_v.setValue(r=escript.Data())
323      def __Q(self,p):           else:
324            return util.tensor_mult(self.__permeability,util.grad(p))             self.__pde_v.setValue(r=u0)
325             u= self.__pde_v.getSolution()
326      def __Aprod(self,dp):      elif  self.solver  == self.SYMSTAB:
327            self.__pde_v.setTolerance(self.getSubProblemTolerance())           gp=util.grad(p)
328            if self.show_details: print "DarcyFlux: Applying operator"           self.__pde_v.setValue(Y=0.5*(util.tensor_mult(self.__permeability_inv,self.__g)-gp),
329            Qdp=self.__Q(dp)                                 X= escript.Data() ,
330            self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())                                 y= escript.Data() )                          
331            du=self.__pde_v.getSolution(verbose=self.show_details)           if u0 == None:
332            return Qdp+du             self.__pde_v.setValue(r=escript.Data())
333      def __inner_GMRES(self,r,s):           else:
334           return util.integrate(util.inner(r,s))             self.__pde_v.setValue(r=u0)
335             u= self.__pde_v.getSolution()
336      def __inner_PCG(self,p,r):      return u
337           return util.integrate(util.inner(self.__Q(p), r))        
338        
339      def __Msolve_PCG(self,r):     def __solve_STAB(self, u0, p0, max_iter, iter_restart):
340            self.__pde_p.setTolerance(self.getSubProblemTolerance())            # p0 is used as an initial guess
341            if self.show_details: print "DarcyFlux: Applying preconditioner"        u=self.getFlux(p0, u0)  
342            self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())            self.__pde_p.setValue( Y=self.__f-util.div(u),
343            return self.__pde_p.getSolution(verbose=self.show_details)                                   X=0.5*(self.__g - u - util.tensor_mult(self.__permeability,util.grad(p0)) ),
344                                     y= escript.Data(),
345      def getFlux(self,p=None, fixed_flux=Data(), show_details=False):                                   r=escript.Data())
346          """  
347          returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}        dp=self.__pde_p.getSolution()
348          on locations where C{location_of_fixed_flux} is positive (see L{setValue}).        p=GMRES(dp,
349          Note that C{g} and C{f} are used, see L{setValue}.                self.__STAB_Aprod,
350              p0,
351          @param p: pressure.            self.__inner,
352          @type p: scalar value on the domain (e.g. L{Data}).            atol=self.__norm(p0+dp)*self.getTolerance() ,
353          @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.            rtol=0.,
354          @type fixed_flux: vector values on the domain (e.g. L{Data}).            iter_max=max_iter,
355          @param tol: relative tolerance to be used.            iter_restart=iter_restart,
356          @type tol: positive C{float}.            verbose=self.verbose,P_R=None)
357          @return: flux              
358          @rtype: L{Data}            u=self.getFlux(p, u0)
359          @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}}            return u,p
360                 for the permeability M{k_{ij}}  
361          """     def __solve_SYMSTAB(self, u0, p0, max_iter, iter_restart):
362          self.__pde_v.setTolerance(self.getSubProblemTolerance())            # p0 is used as an initial guess
363          g=self.__g        u=self.getFlux(p0, u0)  
364          f=self.__f            self.__pde_p.setValue( Y= self.__f,
365          self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)                                   X=  0.5*(self.__g + u - util.tensor_mult(self.__permeability,util.grad(p0)) ),
366          if p == None:                                   y=  -  util.inner(self.domain.getNormal(), u),
367             self.__pde_v.setValue(Y=g)                                   r=escript.Data())
368          else:        dp=self.__pde_p.getSolution()
369             self.__pde_v.setValue(Y=g-self.__Q(p))        
370          return self.__pde_v.getSolution(verbose=show_details)        print dp
371              print p0+dp
372              
373          p=GMRES(dp,
374                  self.__SYMSTAB_Aprod,
375              p0,
376              self.__inner,
377              atol=self.__norm(p0+dp)*self.getTolerance() ,
378              rtol=0.,
379              iter_max=max_iter,
380              iter_restart=iter_restart,
381              verbose=self.verbose,P_R=None)
382                
383              u=self.getFlux(p, u0)
384              return u,p
385    
386       def __L2(self,v):
387             return util.sqrt(util.integrate(util.length(util.interpolate(v,escript.Function(self.domain)))**2))      
388      
389       def __norm(self,r):
390             return util.sqrt(self.__inner(r,r))
391            
392       def __inner(self,r,s):
393             return util.integrate(util.inner(r,s), escript.Function(self.domain))
394            
395       def __STAB_Aprod(self,p):
396          gp=util.grad(p)
397          self.__pde_v.setValue(Y=-0.5*gp,
398                                X=-p*util.kronecker(self.__pde_v.getDomain()),
399                                y= p * self.domain.getNormal(),  
400                                r=escript.Data())
401          u = -self.__pde_v.getSolution()
402          self.__pde_p.setValue(Y=util.div(u),
403                                X=0.5*(u+util.tensor_mult(self.__permeability,gp)),
404                                y=escript.Data(),
405                                r=escript.Data())
406        
407          return  self.__pde_p.getSolution()
408      
409       def __SYMSTAB_Aprod(self,p):
410          gp=util.grad(p)
411          self.__pde_v.setValue(Y=0.5*gp ,
412                                X=escript.Data(),
413                                y=escript.Data(),  
414                                r=escript.Data())
415          u = -self.__pde_v.getSolution()
416          self.__pde_p.setValue(Y=escript.Data(),
417                                X=0.5*(-u+util.tensor_mult(self.__permeability,gp)),
418                                y=escript.Data(),
419                                r=escript.Data())
420        
421          return  self.__pde_p.getSolution()
422          
423    
424  class StokesProblemCartesian(HomogeneousSaddlePointProblem):  class StokesProblemCartesian(HomogeneousSaddlePointProblem):
425       """       """
# Line 391  class StokesProblemCartesian(Homogeneous Line 444  class StokesProblemCartesian(Homogeneous
444           """           """
445           initialize the Stokes Problem           initialize the Stokes Problem
446    
447           @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
448           @type domain: L{Domain}           LBB complient, for instance using quadratic and linear approximation on the same element or using linear approximation
449           @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.           with macro elements for the pressure.
450    
451             :param domain: domain of the problem.
452             :type domain: `Domain`
453           """           """
454           HomogeneousSaddlePointProblem.__init__(self,**kwargs)           HomogeneousSaddlePointProblem.__init__(self,**kwargs)
455           self.domain=domain           self.domain=domain
456           self.vol=util.integrate(1.,Function(self.domain))           self.__pde_v=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
457           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_v.setSymmetryOn()
458           self.__pde_u.setSymmetryOn()      
          # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)  
          # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.RILU)  
   
459           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
460           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
          # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)  
461           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
462    
463       def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1,surface_stress=Data(),stress=Data()):           self.__pde_proj=LinearPDE(domain)
464             self.__pde_proj.setReducedOrderOn()
465         self.__pde_proj.setValue(D=1)
466             self.__pde_proj.setSymmetryOn()
467    
468         def getSolverOptionsVelocity(self):
469             """
470         returns the solver options used  solve the equation for velocity.
471        
472         :rtype: `SolverOptions`
473         """
474         return self.__pde_v.getSolverOptions()
475         def setSolverOptionsVelocity(self, options=None):
476             """
477         set the solver options for solving the equation for velocity.
478        
479         :param options: new solver  options
480         :type options: `SolverOptions`
481         """
482             self.__pde_v.setSolverOptions(options)
483         def getSolverOptionsPressure(self):
484             """
485         returns the solver options used  solve the equation for pressure.
486         :rtype: `SolverOptions`
487         """
488         return self.__pde_prec.getSolverOptions()
489         def setSolverOptionsPressure(self, options=None):
490             """
491         set the solver options for solving the equation for pressure.
492         :param options: new solver  options
493         :type options: `SolverOptions`
494         """
495         self.__pde_prec.setSolverOptions(options)
496    
497         def setSolverOptionsDiv(self, options=None):
498             """
499         set the solver options for solving the equation to project the divergence of
500         the velocity onto the function space of presure.
501        
502         :param options: new solver options
503         :type options: `SolverOptions`
504         """
505         self.__pde_proj.setSolverOptions(options)
506         def getSolverOptionsDiv(self):
507             """
508         returns the solver options for solving the equation to project the divergence of
509         the velocity onto the function space of presure.
510        
511         :rtype: `SolverOptions`
512         """
513         return self.__pde_proj.getSolverOptions()
514    
515         def updateStokesEquation(self, v, p):
516             """
517             updates the Stokes equation to consider dependencies from ``v`` and ``p``
518             :note: This method can be overwritten by a subclass. Use `setStokesEquation` to set new values.
519             """
520             pass
521         def setStokesEquation(self, f=None,fixed_u_mask=None,eta=None,surface_stress=None,stress=None, restoration_factor=None):
522            """
523            assigns new values to the model parameters.
524    
525            :param f: external force
526            :type f: `Vector` object in `FunctionSpace` `Function` or similar
527            :param fixed_u_mask: mask of locations with fixed velocity.
528            :type fixed_u_mask: `Vector` object on `FunctionSpace` `Solution` or similar
529            :param eta: viscosity
530            :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
531            :param surface_stress: normal surface stress
532            :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
533            :param stress: initial stress
534        :type stress: `Tensor` object on `FunctionSpace` `Function` or similar
535            """
536            if eta !=None:
537                k=util.kronecker(self.domain.getDim())
538                kk=util.outer(k,k)
539                self.eta=util.interpolate(eta, escript.Function(self.domain))
540            self.__pde_prec.setValue(D=1/self.eta)
541                self.__pde_v.setValue(A=self.eta*(util.swap_axes(kk,0,3)+util.swap_axes(kk,1,3)))
542            if restoration_factor!=None:
543                n=self.domain.getNormal()
544                self.__pde_v.setValue(d=restoration_factor*util.outer(n,n))
545            if fixed_u_mask!=None:
546                self.__pde_v.setValue(q=fixed_u_mask)
547            if f!=None: self.__f=f
548            if surface_stress!=None: self.__surface_stress=surface_stress
549            if stress!=None: self.__stress=stress
550    
551         def initialize(self,f=escript.Data(),fixed_u_mask=escript.Data(),eta=1,surface_stress=escript.Data(),stress=escript.Data(), restoration_factor=0):
552          """          """
553          assigns values to the model parameters          assigns values to the model parameters
554    
555          @param f: external force          :param f: external force
556          @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar          :type f: `Vector` object in `FunctionSpace` `Function` or similar
557          @param fixed_u_mask: mask of locations with fixed velocity.          :param fixed_u_mask: mask of locations with fixed velocity.
558          @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
559          @param eta: viscosity          :param eta: viscosity
560          @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar          :type eta: `Scalar` object on `FunctionSpace` `Function` or similar
561          @param surface_stress: normal surface stress          :param surface_stress: normal surface stress
562          @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar          :type surface_stress: `Vector` object on `FunctionSpace` `FunctionOnBoundary` or similar
563          @param stress: initial stress          :param stress: initial stress
564      @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.  
   
565          """          """
566          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  
567    
568       def inner_pBv(self,p,v):       def Bv(self,v,tol):
569           """           """
570           returns inner product of element p and div(v)           returns inner product of element p and div(v)
571    
572           @param p: a pressure increment           :param v: a residual
573           @param v: a residual           :return: inner product of element p and div(v)
574           @return: inner product of element p and div(v)           :rtype: ``float``
575           @rtype: C{float}           """
576             self.__pde_proj.setValue(Y=-util.div(v))
577         self.getSolverOptionsDiv().setTolerance(tol)
578         self.getSolverOptionsDiv().setAbsoluteTolerance(0.)
579             out=self.__pde_proj.getSolution()
580             return out
581    
582         def inner_pBv(self,p,Bv):
583             """
584             returns inner product of element p and Bv=-div(v)
585    
586             :param p: a pressure increment
587             :param Bv: a residual
588             :return: inner product of element p and Bv=-div(v)
589             :rtype: ``float``
590           """           """
591           return util.integrate(-p*util.div(v))           return util.integrate(util.interpolate(p,escript.Function(self.domain))*util.interpolate(Bv, escript.Function(self.domain)))
592    
593       def inner_p(self,p0,p1):       def inner_p(self,p0,p1):
594           """           """
595           Returns inner product of p0 and p1           Returns inner product of p0 and p1
596    
597           @param p0: a pressure           :param p0: a pressure
598           @param p1: a pressure           :param p1: a pressure
599           @return: inner product of p0 and p1           :return: inner product of p0 and p1
600           @rtype: C{float}           :rtype: ``float``
601           """           """
602           s0=util.interpolate(p0/self.eta,Function(self.domain))           s0=util.interpolate(p0, escript.Function(self.domain))
603           s1=util.interpolate(p1/self.eta,Function(self.domain))           s1=util.interpolate(p1, escript.Function(self.domain))
604           return util.integrate(s0*s1)           return util.integrate(s0*s1)
605    
606       def norm_v(self,v):       def norm_v(self,v):
607           """           """
608           returns the norm of v           returns the norm of v
609    
610           @param v: a velovity           :param v: a velovity
611           @return: norm of v           :return: norm of v
612           @rtype: non-negative C{float}           :rtype: non-negative ``float``
613           """           """
614           return util.sqrt(util.integrate(util.length(util.grad(v))))           return util.sqrt(util.integrate(util.length(util.grad(v))**2))
615    
616       def getV(self, p, v0):  
617         def getDV(self, p, v, tol):
618           """           """
619           return the value for v for a given p (overwrite)           return the value for v for a given p (overwrite)
620    
621           @param p: a pressure           :param p: a pressure
622           @param v0: a initial guess for the value v to return.           :param v: a initial guess for the value v to return.
623           @return: v given as M{v= A^{-1} (f-B^*p)}           :return: dv given as *Adv=(f-Av-B^*p)*
624           """           """
625           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.updateStokesEquation(v,p)
626           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)           self.__pde_v.setValue(Y=self.__f, y=self.__surface_stress)
627         self.getSolverOptionsVelocity().setTolerance(tol)
628         self.getSolverOptionsVelocity().setAbsoluteTolerance(0.)
629           if self.__stress.isEmpty():           if self.__stress.isEmpty():
630              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)))
631           else:           else:
632              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)))
633           out=self.__pde_u.getSolution(verbose=self.show_details)           out=self.__pde_v.getSolution()
634           return  out           return  out
635    
636         def norm_Bv(self,Bv):
          raise NotImplementedError,"no v calculation implemented."  
   
   
      def norm_Bv(self,v):  
637          """          """
638          Returns Bv (overwrite).          Returns Bv (overwrite).
639    
640          @rtype: equal to the type of p          :rtype: equal to the type of p
641          @note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
642          """          """
643          return util.sqrt(util.integrate(util.div(v)**2))          return util.sqrt(util.integrate(util.interpolate(Bv, escript.Function(self.domain))**2))
644    
645       def solve_AinvBt(self,p):       def solve_AinvBt(self,p, tol):
646           """           """
647           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}           Solves *Av=B^*p* with accuracy `tol`
648    
649           @param p: a pressure increment           :param p: a pressure increment
650           @return: the solution of M{Av=B^*p}           :return: the solution of *Av=B^*p*
651           @note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
652           """           """
653           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_v.setValue(Y=escript.Data(), y=escript.Data(), X=-p*util.kronecker(self.domain))
654           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)  
655           return  out           return  out
656    
657       def solve_precB(self,v):       def solve_prec(self,Bv, tol):
658           """           """
659           applies preconditioner for for M{BA^{-1}B^*} to M{Bv}           applies preconditioner for for *BA^{-1}B^** to *Bv*
660           with accuracy L{self.getSubProblemTolerance()} (overwrite).           with accuracy `self.getSubProblemTolerance()`
661    
662           @param v: velocity increment           :param Bv: velocity increment
663           @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^ * )*
664           @note: boundary conditions on p are zero.           :note: boundary conditions on p are zero.
665           """           """
666           self.__pde_prec.setValue(Y=-util.div(v))           self.__pde_prec.setValue(Y=Bv)
667           self.__pde_prec.setTolerance(self.getSubProblemTolerance())       self.getSolverOptionsPressure().setTolerance(tol)
668           return self.__pde_prec.getSolution(verbose=self.show_details)       self.getSolverOptionsPressure().setAbsoluteTolerance(0.)
669             out=self.__pde_prec.getSolution()
670             return out

Legend:
Removed from v.2344  
changed lines
  Added in v.3510

  ViewVC Help
Powered by ViewVC 1.1.26