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

Legend:
Removed from v.2370  
changed lines
  Added in v.3515

  ViewVC Help
Powered by ViewVC 1.1.26