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

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

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

revision 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC revision 2676 by gross, Mon Sep 21 08:06:37 2009 UTC
# Line 11  Line 11 
11  #  #
12  ########################################################  ########################################################
13    
14  __copyright__="""Copyright (c) 2003-2008 by University of Queensland  __copyright__="""Copyright (c) 2003-2009 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 28  Currently includes: Line 28  Currently includes:
28      - TimeIntegrationManager - to handle extrapolation in time      - TimeIntegrationManager - to handle extrapolation in time
29      - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme      - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme
30    
31  @var __author__: name of author  :var __author__: name of author
32  @var __copyright__: copyrights  :var __copyright__: copyrights
33  @var __license__: licence agreement  :var __license__: licence agreement
34  @var __url__: url entry point on documentation  :var __url__: url entry point on documentation
35  @var __version__: version  :var __version__: version
36  @var __date__: date of the version  :var __date__: date of the version
37  """  """
38    
39  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
# Line 45  import numpy Line 45  import numpy
45  import util  import util
46  import math  import math
47    
 ##### Added by Artak  
 # from Numeric import zeros,Int,Float64  
 ###################################  
   
   
48  class TimeIntegrationManager:  class TimeIntegrationManager:
49    """    """
50    A simple mechanism to manage time dependend values.    A simple mechanism to manage time dependend values.
# Line 64  class TimeIntegrationManager: Line 59  class TimeIntegrationManager:
59           tm.checkin(dt,v)           tm.checkin(dt,v)
60           t+=dt           t+=dt
61    
62    @note: currently only p=1 is supported.    :note: currently only p=1 is supported.
63    """    """
64    def __init__(self,*inital_values,**kwargs):    def __init__(self,*inital_values,**kwargs):
65       """       """
66       Sets up the value manager where C{inital_values} are the initial values       Sets up the value manager where ``inital_values`` are the initial values
67       and p is the order used for extrapolation.       and p is the order used for extrapolation.
68       """       """
69       if kwargs.has_key("p"):       if kwargs.has_key("p"):
# Line 113  class TimeIntegrationManager: Line 108  class TimeIntegrationManager:
108    
109    def extrapolate(self,dt):    def extrapolate(self,dt):
110        """        """
111        Extrapolates to C{dt} forward in time.        Extrapolates to ``dt`` forward in time.
112        """        """
113        if self.__order==0:        if self.__order==0:
114           out=self.__v_mem[0]           out=self.__v_mem[0]
# Line 139  class Projector: Line 134  class Projector:
134      """      """
135      Creates a continuous function space projector for a domain.      Creates a continuous function space projector for a domain.
136    
137      @param domain: Domain of the projection.      :param domain: Domain of the projection.
138      @param reduce: Flag to reduce projection order      :param reduce: Flag to reduce projection order
139      @param fast: Flag to use a fast method based on matrix lumping      :param fast: Flag to use a fast method based on matrix lumping
140      """      """
141      self.__pde = linearPDEs.LinearPDE(domain)      self.__pde = linearPDEs.LinearPDE(domain)
142      if fast:      if fast:
# Line 154  class Projector: Line 149  class Projector:
149      """      """
150      Returns the solver options of the PDE solver.      Returns the solver options of the PDE solver.
151            
152      @rtype: L{linearPDEs.SolverOptions}      :rtype: `linearPDEs.SolverOptions`
153      """      """
154    
155    def __call__(self, input_data):    def __call__(self, input_data):
156      """      """
157      Projects C{input_data} onto a continuous function.      Projects ``input_data`` onto a continuous function.
158    
159      @param input_data: the data to be projected      :param input_data: the data to be projected
160      """      """
161      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
162      self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())      self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())
# Line 196  class NoPDE: Line 191  class NoPDE:
191       """       """
192       Solves the following problem for u:       Solves the following problem for u:
193    
194       M{kronecker[i,j]*D[j]*u[j]=Y[i]}       *kronecker[i,j]*D[j]*u[j]=Y[i]*
195    
196       with constraint       with constraint
197    
198       M{u[j]=r[j]}  where M{q[j]>0}       *u[j]=r[j]*  where *q[j]>0*
199    
200       where M{D}, M{Y}, M{r} and M{q} are given functions of rank 1.       where *D*, *Y*, *r* and *q* are given functions of rank 1.
201    
202       In the case of scalars this takes the form       In the case of scalars this takes the form
203    
204       M{D*u=Y}       *D*u=Y*
205    
206       with constraint       with constraint
207    
208       M{u=r} where M{q>0}       *u=r* where *q>0*
209    
210       where M{D}, M{Y}, M{r} and M{q} are given scalar functions.       where *D*, *Y*, *r* and *q* are given scalar functions.
211    
212       The constraint overwrites any other condition.       The constraint overwrites any other condition.
213    
214       @note: This class is similar to the L{linearPDEs.LinearPDE} class with       :note: This class is similar to the `linearPDEs.LinearPDE` class with
215              A=B=C=X=0 but has the intention that all input parameters are given              A=B=C=X=0 but has the intention that all input parameters are given
216              in L{Solution} or L{ReducedSolution}.              in `Solution` or `ReducedSolution`.
217       """       """
218       # The whole thing is a bit strange and I blame Rob Woodcock (CSIRO) for       # The whole thing is a bit strange and I blame Rob Woodcock (CSIRO) for
219       # this.       # this.
# Line 226  class NoPDE: Line 221  class NoPDE:
221           """           """
222           Initializes the problem.           Initializes the problem.
223    
224           @param domain: domain of the PDE           :param domain: domain of the PDE
225           @type domain: L{Domain}           :type domain: `Domain`
226           @param D: coefficient of the solution           :param D: coefficient of the solution
227           @type D: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type D: ``float``, ``int``, ``numpy.ndarray``, `Data`
228           @param Y: right hand side           :param Y: right hand side
229           @type Y: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type Y: ``float``, ``int``, ``numpy.ndarray``, `Data`
230           @param q: location of constraints           :param q: location of constraints
231           @type q: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type q: ``float``, ``int``, ``numpy.ndarray``, `Data`
232           @param r: value of solution at locations of constraints           :param r: value of solution at locations of constraints
233           @type r: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type r: ``float``, ``int``, ``numpy.ndarray``, `Data`
234           """           """
235           self.__domain=domain           self.__domain=domain
236           self.__D=D           self.__D=D
# Line 247  class NoPDE: Line 242  class NoPDE:
242    
243       def setReducedOn(self):       def setReducedOn(self):
244           """           """
245           Sets the L{FunctionSpace} of the solution to L{ReducedSolution}.           Sets the `FunctionSpace` of the solution to `ReducedSolution`.
246           """           """
247           self.__function_space=escript.ReducedSolution(self.__domain)           self.__function_space=escript.ReducedSolution(self.__domain)
248           self.__u=None           self.__u=None
249    
250       def setReducedOff(self):       def setReducedOff(self):
251           """           """
252           Sets the L{FunctionSpace} of the solution to L{Solution}.           Sets the `FunctionSpace` of the solution to `Solution`.
253           """           """
254           self.__function_space=escript.Solution(self.__domain)           self.__function_space=escript.Solution(self.__domain)
255           self.__u=None           self.__u=None
# Line 263  class NoPDE: Line 258  class NoPDE:
258           """           """
259           Assigns values to the parameters.           Assigns values to the parameters.
260    
261           @param D: coefficient of the solution           :param D: coefficient of the solution
262           @type D: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type D: ``float``, ``int``, ``numpy.ndarray``, `Data`
263           @param Y: right hand side           :param Y: right hand side
264           @type Y: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type Y: ``float``, ``int``, ``numpy.ndarray``, `Data`
265           @param q: location of constraints           :param q: location of constraints
266           @type q: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type q: ``float``, ``int``, ``numpy.ndarray``, `Data`
267           @param r: value of solution at locations of constraints           :param r: value of solution at locations of constraints
268           @type r: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type r: ``float``, ``int``, ``numpy.ndarray``, `Data`
269           """           """
270           if not D==None:           if not D==None:
271              self.__D=D              self.__D=D
# Line 289  class NoPDE: Line 284  class NoPDE:
284           """           """
285           Returns the solution.           Returns the solution.
286    
287           @return: the solution of the problem           :return: the solution of the problem
288           @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or           :rtype: `Data` object in the `FunctionSpace` `Solution` or
289                   L{ReducedSolution}                   `ReducedSolution`
290           """           """
291           if self.__u==None:           if self.__u==None:
292              if self.__D==None:              if self.__D==None:
# Line 324  class Locator: Line 319  class Locator:
319         or FunctionSpace for the sample point which is closest to the given         or FunctionSpace for the sample point which is closest to the given
320         point x.         point x.
321    
322         @param where: function space         :param where: function space
323         @type where: L{escript.FunctionSpace}         :type where: `escript.FunctionSpace`
324         @param x: location(s) of the Locator         :param x: location(s) of the Locator
325         @type x: C{numpy.ndarray} or C{list} of C{numpy.ndarray}         :type x: ``numpy.ndarray`` or ``list`` of ``numpy.ndarray``
326         """         """
327         if isinstance(where,escript.FunctionSpace):         if isinstance(where,escript.FunctionSpace):
328            self.__function_space=where            self.__function_space=where
# Line 354  class Locator: Line 349  class Locator:
349         Returns the coordinates of the Locator as a string.         Returns the coordinates of the Locator as a string.
350         """         """
351         x=self.getX()         x=self.getX()
352         if instance(x,list):         if isinstance(x,list):
353            out="["            out="["
354            first=True            first=True
355            for xx in x:            for xx in x:
# Line 401  class Locator: Line 396  class Locator:
396    
397       def getValue(self,data):       def getValue(self,data):
398          """          """
399          Returns the value of C{data} at the Locator if C{data} is a L{Data}          Returns the value of ``data`` at the Locator if ``data`` is a `Data`
400          object otherwise the object is returned.          object otherwise the object is returned.
401          """          """
402          if isinstance(data,escript.Data):          if isinstance(data,escript.Data):
# Line 467  def PCG(r, Aprod, x, Msolve, bilinearfor Line 462  def PCG(r, Aprod, x, Msolve, bilinearfor
462     """     """
463     Solver for     Solver for
464    
465     M{Ax=b}     *Ax=b*
466    
467     with a symmetric and positive definite operator A (more details required!).     with a symmetric and positive definite operator A (more details required!).
468     It uses the conjugate gradient method with preconditioner M providing an     It uses the conjugate gradient method with preconditioner M providing an
# Line 475  def PCG(r, Aprod, x, Msolve, bilinearfor Line 470  def PCG(r, Aprod, x, Msolve, bilinearfor
470    
471     The iteration is terminated if     The iteration is terminated if
472    
473     M{|r| <= atol+rtol*|r0|}     *|r| <= atol+rtol*|r0|*
474    
475     where M{r0} is the initial residual and M{|.|} is the energy norm. In fact     where *r0* is the initial residual and *|.|* is the energy norm. In fact
476    
477     M{|r| = sqrt( bilinearform(Msolve(r),r))}     *|r| = sqrt( bilinearform(Msolve(r),r))*
478    
479     For details on the preconditioned conjugate gradient method see the book:     For details on the preconditioned conjugate gradient method see the book:
480    
# Line 487  def PCG(r, Aprod, x, Msolve, bilinearfor Line 482  def PCG(r, Aprod, x, Msolve, bilinearfor
482     T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,     T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
483     C. Romine, and H. van der Vorst}.     C. Romine, and H. van der Vorst}.
484    
485     @param r: initial residual M{r=b-Ax}. C{r} is altered.     :param r: initial residual *r=b-Ax*. ``r`` is altered.
486     @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)     :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
487     @param x: an initial guess for the solution     :param x: an initial guess for the solution
488     @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)     :type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
489     @param Aprod: returns the value Ax     :param Aprod: returns the value Ax
490     @type Aprod: function C{Aprod(x)} where C{x} is of the same object like     :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
491                  argument C{x}. The returned object needs to be of the same type                  argument ``x``. The returned object needs to be of the same type
492                  like argument C{r}.                  like argument ``r``.
493     @param Msolve: solves Mx=r     :param Msolve: solves Mx=r
494     @type Msolve: function C{Msolve(r)} where C{r} is of the same type like     :type Msolve: function ``Msolve(r)`` where ``r`` is of the same type like
495                   argument C{r}. The returned object needs to be of the same                   argument ``r``. The returned object needs to be of the same
496                   type like argument C{x}.                   type like argument ``x``.
497     @param bilinearform: inner product C{<x,r>}     :param bilinearform: inner product ``<x,r>``
498     @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same     :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
499                         type like argument C{x} and C{r} is. The returned value                         type like argument ``x`` and ``r`` is. The returned value
500                         is a C{float}.                         is a ``float``.
501     @param atol: absolute tolerance     :param atol: absolute tolerance
502     @type atol: non-negative C{float}     :type atol: non-negative ``float``
503     @param rtol: relative tolerance     :param rtol: relative tolerance
504     @type rtol: non-negative C{float}     :type rtol: non-negative ``float``
505     @param iter_max: maximum number of iteration steps     :param iter_max: maximum number of iteration steps
506     @type iter_max: C{int}     :type iter_max: ``int``
507     @return: the solution approximation and the corresponding residual     :return: the solution approximation and the corresponding residual
508     @rtype: C{tuple}     :rtype: ``tuple``
509     @warning: C{r} and C{x} are altered.     :warning: ``r`` and ``x`` are altered.
510     """     """
511     iter=0     iter=0
512     rhat=Msolve(r)     rhat=Msolve(r)
# Line 564  class Defect(object): Line 559  class Defect(object):
559          """          """
560          Returns the inner product of x0 and x1          Returns the inner product of x0 and x1
561    
562          @param x0: value for x0          :param x0: value for x0
563          @param x1: value for x1          :param x1: value for x1
564          @return: the inner product of x0 and x1          :return: the inner product of x0 and x1
565          @rtype: C{float}          :rtype: ``float``
566          """          """
567          return 0          return 0
568    
569      def norm(self,x):      def norm(self,x):
570          """          """
571          Returns the norm of argument C{x}.          Returns the norm of argument ``x``.
572    
573          @param x: a value          :param x: a value
574          @return: norm of argument x          :return: norm of argument x
575          @rtype: C{float}          :rtype: ``float``
576          @note: by default C{sqrt(self.bilinearform(x,x)} is returned.          :note: by default ``sqrt(self.bilinearform(x,x)`` is returned.
577          """          """
578          s=self.bilinearform(x,x)          s=self.bilinearform(x,x)
579          if s<0: raise NegativeNorm,"negative norm."          if s<0: raise NegativeNorm,"negative norm."
# Line 586  class Defect(object): Line 581  class Defect(object):
581    
582      def eval(self,x):      def eval(self,x):
583          """          """
584          Returns the value F of a given C{x}.          Returns the value F of a given ``x``.
585    
586          @param x: value for which the defect C{F} is evaluated          :param x: value for which the defect ``F`` is evaluated
587          @return: value of the defect at C{x}          :return: value of the defect at ``x``
588          """          """
589          return 0          return 0
590    
# Line 602  class Defect(object): Line 597  class Defect(object):
597          derivative of the defect. The increment is inc*norm(x)/norm(v)*v in the          derivative of the defect. The increment is inc*norm(x)/norm(v)*v in the
598          direction of v with x as a starting point.          direction of v with x as a starting point.
599    
600          @param inc: relative increment length          :param inc: relative increment length
601          @type inc: positive C{float}          :type inc: positive ``float``
602          """          """
603          if inc<=0: raise ValueError,"positive increment required."          if inc<=0: raise ValueError,"positive increment required."
604          self.__inc=inc          self.__inc=inc
# Line 612  class Defect(object): Line 607  class Defect(object):
607          """          """
608          Returns the relative increment length used to approximate the          Returns the relative increment length used to approximate the
609          derivative of the defect.          derivative of the defect.
610          @return: value of the defect at C{x}          :return: value of the defect at ``x``
611          @rtype: positive C{float}          :rtype: positive ``float``
612          """          """
613          return self.__inc          return self.__inc
614    
615      def derivative(self, F0, x0, v, v_is_normalised=True):      def derivative(self, F0, x0, v, v_is_normalised=True):
616          """          """
617          Returns the directional derivative at C{x0} in the direction of C{v}.          Returns the directional derivative at ``x0`` in the direction of ``v``.
618    
619          @param F0: value of this defect at x0          :param F0: value of this defect at x0
620          @param x0: value at which derivative is calculated          :param x0: value at which derivative is calculated
621          @param v: direction          :param v: direction
622          @param v_is_normalised: True to indicate that C{v} is nomalized          :param v_is_normalised: True to indicate that ``v`` is nomalized
623                                  (self.norm(v)=0)                                  (self.norm(v)=0)
624          @return: derivative of this defect at x0 in the direction of C{v}          :return: derivative of this defect at x0 in the direction of ``v``
625          @note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is          :note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is
626                 used but this method maybe overwritten to use exact evaluation.                 used but this method maybe overwritten to use exact evaluation.
627          """          """
628          normx=self.norm(x0)          normx=self.norm(x0)
# Line 647  class Defect(object): Line 642  class Defect(object):
642  ######################################  ######################################
643  def NewtonGMRES(defect, x, iter_max=100, sub_iter_max=20, atol=0,rtol=1.e-4, sub_tol_max=0.5, gamma=0.9, verbose=False):  def NewtonGMRES(defect, x, iter_max=100, sub_iter_max=20, atol=0,rtol=1.e-4, sub_tol_max=0.5, gamma=0.9, verbose=False):
644     """     """
645     Solves a non-linear problem M{F(x)=0} for unknown M{x} using the stopping     Solves a non-linear problem *F(x)=0* for unknown *x* using the stopping
646     criterion:     criterion:
647    
648     M{norm(F(x) <= atol + rtol * norm(F(x0)}     *norm(F(x) <= atol + rtol * norm(F(x0)*
649    
650     where M{x0} is the initial guess.     where *x0* is the initial guess.
651    
652     @param defect: object defining the function M{F}. C{defect.norm} defines the     :param defect: object defining the function *F*. ``defect.norm`` defines the
653                    M{norm} used in the stopping criterion.                    *norm* used in the stopping criterion.
654     @type defect: L{Defect}     :type defect: `Defect`
655     @param x: initial guess for the solution, C{x} is altered.     :param x: initial guess for the solution, ``x`` is altered.
656     @type x: any object type allowing basic operations such as     :type x: any object type allowing basic operations such as
657              C{numpy.ndarray}, L{Data}              ``numpy.ndarray``, `Data`
658     @param iter_max: maximum number of iteration steps     :param iter_max: maximum number of iteration steps
659     @type iter_max: positive C{int}     :type iter_max: positive ``int``
660     @param sub_iter_max: maximum number of inner iteration steps     :param sub_iter_max: maximum number of inner iteration steps
661     @type sub_iter_max: positive C{int}     :type sub_iter_max: positive ``int``
662     @param atol: absolute tolerance for the solution     :param atol: absolute tolerance for the solution
663     @type atol: positive C{float}     :type atol: positive ``float``
664     @param rtol: relative tolerance for the solution     :param rtol: relative tolerance for the solution
665     @type rtol: positive C{float}     :type rtol: positive ``float``
666     @param gamma: tolerance safety factor for inner iteration     :param gamma: tolerance safety factor for inner iteration
667     @type gamma: positive C{float}, less than 1     :type gamma: positive ``float``, less than 1
668     @param sub_tol_max: upper bound for inner tolerance     :param sub_tol_max: upper bound for inner tolerance
669     @type sub_tol_max: positive C{float}, less than 1     :type sub_tol_max: positive ``float``, less than 1
670     @return: an approximation of the solution with the desired accuracy     :return: an approximation of the solution with the desired accuracy
671     @rtype: same type as the initial guess     :rtype: same type as the initial guess
672     """     """
673     lmaxit=iter_max     lmaxit=iter_max
674     if atol<0: raise ValueError,"atol needs to be non-negative."     if atol<0: raise ValueError,"atol needs to be non-negative."
# Line 731  def NewtonGMRES(defect, x, iter_max=100, Line 726  def NewtonGMRES(defect, x, iter_max=100,
726  def __givapp(c,s,vin):  def __givapp(c,s,vin):
727      """      """
728      Applies a sequence of Givens rotations (c,s) recursively to the vector      Applies a sequence of Givens rotations (c,s) recursively to the vector
729      C{vin}      ``vin``
730    
731      @warning: C{vin} is altered.      :warning: ``vin`` is altered.
732      """      """
733      vrot=vin      vrot=vin
734      if isinstance(c,float):      if isinstance(c,float):
# Line 838  def GMRES(r, Aprod, x, bilinearform, ato Line 833  def GMRES(r, Aprod, x, bilinearform, ato
833     """     """
834     Solver for     Solver for
835    
836     M{Ax=b}     *Ax=b*
837    
838     with a general operator A (more details required!).     with a general operator A (more details required!).
839     It uses the generalized minimum residual method (GMRES).     It uses the generalized minimum residual method (GMRES).
840    
841     The iteration is terminated if     The iteration is terminated if
842    
843     M{|r| <= atol+rtol*|r0|}     *|r| <= atol+rtol*|r0|*
844    
845     where M{r0} is the initial residual and M{|.|} is the energy norm. In fact     where *r0* is the initial residual and *|.|* is the energy norm. In fact
846    
847     M{|r| = sqrt( bilinearform(r,r))}     *|r| = sqrt( bilinearform(r,r))*
848    
849     @param r: initial residual M{r=b-Ax}. C{r} is altered.     :param r: initial residual *r=b-Ax*. ``r`` is altered.
850     @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)     :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
851     @param x: an initial guess for the solution     :param x: an initial guess for the solution
852     @type x: same like C{r}     :type x: same like ``r``
853     @param Aprod: returns the value Ax     :param Aprod: returns the value Ax
854     @type Aprod: function C{Aprod(x)} where C{x} is of the same object like     :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
855                  argument C{x}. The returned object needs to be of the same                  argument ``x``. The returned object needs to be of the same
856                  type like argument C{r}.                  type like argument ``r``.
857     @param bilinearform: inner product C{<x,r>}     :param bilinearform: inner product ``<x,r>``
858     @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same     :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
859                         type like argument C{x} and C{r}. The returned value is                         type like argument ``x`` and ``r``. The returned value is
860                         a C{float}.                         a ``float``.
861     @param atol: absolute tolerance     :param atol: absolute tolerance
862     @type atol: non-negative C{float}     :type atol: non-negative ``float``
863     @param rtol: relative tolerance     :param rtol: relative tolerance
864     @type rtol: non-negative C{float}     :type rtol: non-negative ``float``
865     @param iter_max: maximum number of iteration steps     :param iter_max: maximum number of iteration steps
866     @type iter_max: C{int}     :type iter_max: ``int``
867     @param iter_restart: in order to save memory the orthogonalization process     :param iter_restart: in order to save memory the orthogonalization process
868                          is terminated after C{iter_restart} steps and the                          is terminated after ``iter_restart`` steps and the
869                          iteration is restarted.                          iteration is restarted.
870     @type iter_restart: C{int}     :type iter_restart: ``int``
871     @return: the solution approximation and the corresponding residual     :return: the solution approximation and the corresponding residual
872     @rtype: C{tuple}     :rtype: ``tuple``
873     @warning: C{r} and C{x} are altered.     :warning: ``r`` and ``x`` are altered.
874     """     """
875     m=iter_restart     m=iter_restart
876     restarted=False     restarted=False
# Line 1007  def MINRES(r, Aprod, x, Msolve, bilinear Line 1002  def MINRES(r, Aprod, x, Msolve, bilinear
1002      """      """
1003      Solver for      Solver for
1004    
1005      M{Ax=b}      *Ax=b*
1006    
1007      with a symmetric and positive definite operator A (more details required!).      with a symmetric and positive definite operator A (more details required!).
1008      It uses the minimum residual method (MINRES) with preconditioner M      It uses the minimum residual method (MINRES) with preconditioner M
# Line 1015  def MINRES(r, Aprod, x, Msolve, bilinear Line 1010  def MINRES(r, Aprod, x, Msolve, bilinear
1010    
1011      The iteration is terminated if      The iteration is terminated if
1012    
1013      M{|r| <= atol+rtol*|r0|}      *|r| <= atol+rtol*|r0|*
1014    
1015      where M{r0} is the initial residual and M{|.|} is the energy norm. In fact      where *r0* is the initial residual and *|.|* is the energy norm. In fact
1016    
1017      M{|r| = sqrt( bilinearform(Msolve(r),r))}      *|r| = sqrt( bilinearform(Msolve(r),r))*
1018    
1019      For details on the preconditioned conjugate gradient method see the book:      For details on the preconditioned conjugate gradient method see the book:
1020    
# Line 1027  def MINRES(r, Aprod, x, Msolve, bilinear Line 1022  def MINRES(r, Aprod, x, Msolve, bilinear
1022      T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,      T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
1023      C. Romine, and H. van der Vorst}.      C. Romine, and H. van der Vorst}.
1024    
1025      @param r: initial residual M{r=b-Ax}. C{r} is altered.      :param r: initial residual *r=b-Ax*. ``r`` is altered.
1026      @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)      :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1027      @param x: an initial guess for the solution      :param x: an initial guess for the solution
1028      @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)      :type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1029      @param Aprod: returns the value Ax      :param Aprod: returns the value Ax
1030      @type Aprod: function C{Aprod(x)} where C{x} is of the same object like      :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
1031                   argument C{x}. The returned object needs to be of the same                   argument ``x``. The returned object needs to be of the same
1032                   type like argument C{r}.                   type like argument ``r``.
1033      @param Msolve: solves Mx=r      :param Msolve: solves Mx=r
1034      @type Msolve: function C{Msolve(r)} where C{r} is of the same type like      :type Msolve: function ``Msolve(r)`` where ``r`` is of the same type like
1035                    argument C{r}. The returned object needs to be of the same                    argument ``r``. The returned object needs to be of the same
1036                    type like argument C{x}.                    type like argument ``x``.
1037      @param bilinearform: inner product C{<x,r>}      :param bilinearform: inner product ``<x,r>``
1038      @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same      :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
1039                          type like argument C{x} and C{r} is. The returned value                          type like argument ``x`` and ``r`` is. The returned value
1040                          is a C{float}.                          is a ``float``.
1041      @param atol: absolute tolerance      :param atol: absolute tolerance
1042      @type atol: non-negative C{float}      :type atol: non-negative ``float``
1043      @param rtol: relative tolerance      :param rtol: relative tolerance
1044      @type rtol: non-negative C{float}      :type rtol: non-negative ``float``
1045      @param iter_max: maximum number of iteration steps      :param iter_max: maximum number of iteration steps
1046      @type iter_max: C{int}      :type iter_max: ``int``
1047      @return: the solution approximation and the corresponding residual      :return: the solution approximation and the corresponding residual
1048      @rtype: C{tuple}      :rtype: ``tuple``
1049      @warning: C{r} and C{x} are altered.      :warning: ``r`` and ``x`` are altered.
1050      """      """
1051      #------------------------------------------------------------------      #------------------------------------------------------------------
1052      # Set up y and v for the first Lanczos vector v1.      # Set up y and v for the first Lanczos vector v1.
# Line 1185  def TFQMR(r, Aprod, x, bilinearform, ato Line 1180  def TFQMR(r, Aprod, x, bilinearform, ato
1180    """    """
1181    Solver for    Solver for
1182    
1183    M{Ax=b}    *Ax=b*
1184    
1185    with a general operator A (more details required!).    with a general operator A (more details required!).
1186    It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).    It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).
1187    
1188    The iteration is terminated if    The iteration is terminated if
1189    
1190    M{|r| <= atol+rtol*|r0|}    *|r| <= atol+rtol*|r0|*
1191    
1192    where M{r0} is the initial residual and M{|.|} is the energy norm. In fact    where *r0* is the initial residual and *|.|* is the energy norm. In fact
1193    
1194    M{|r| = sqrt( bilinearform(r,r))}    *|r| = sqrt( bilinearform(r,r))*
1195    
1196    @param r: initial residual M{r=b-Ax}. C{r} is altered.    :param r: initial residual *r=b-Ax*. ``r`` is altered.
1197    @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)    :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1198    @param x: an initial guess for the solution    :param x: an initial guess for the solution
1199    @type x: same like C{r}    :type x: same like ``r``
1200    @param Aprod: returns the value Ax    :param Aprod: returns the value Ax
1201    @type Aprod: function C{Aprod(x)} where C{x} is of the same object like    :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
1202                 argument C{x}. The returned object needs to be of the same type                 argument ``x``. The returned object needs to be of the same type
1203                 like argument C{r}.                 like argument ``r``.
1204    @param bilinearform: inner product C{<x,r>}    :param bilinearform: inner product ``<x,r>``
1205    @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same    :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
1206                        type like argument C{x} and C{r}. The returned value is                        type like argument ``x`` and ``r``. The returned value is
1207                        a C{float}.                        a ``float``.
1208    @param atol: absolute tolerance    :param atol: absolute tolerance
1209    @type atol: non-negative C{float}    :type atol: non-negative ``float``
1210    @param rtol: relative tolerance    :param rtol: relative tolerance
1211    @type rtol: non-negative C{float}    :type rtol: non-negative ``float``
1212    @param iter_max: maximum number of iteration steps    :param iter_max: maximum number of iteration steps
1213    @type iter_max: C{int}    :type iter_max: ``int``
1214    @rtype: C{tuple}    :rtype: ``tuple``
1215    @warning: C{r} and C{x} are altered.    :warning: ``r`` and ``x`` are altered.
1216    """    """
1217    u1=0    u1=0
1218    u2=0    u2=0
# Line 1287  def TFQMR(r, Aprod, x, bilinearform, ato Line 1282  def TFQMR(r, Aprod, x, bilinearform, ato
1282    
1283  class ArithmeticTuple(object):  class ArithmeticTuple(object):
1284     """     """
1285     Tuple supporting inplace update x+=y and scaling x=a*y where C{x,y} is an     Tuple supporting inplace update x+=y and scaling x=a*y where ``x,y`` is an
1286     ArithmeticTuple and C{a} is a float.     ArithmeticTuple and ``a`` is a float.
1287    
1288     Example of usage::     Example of usage::
1289    
# Line 1302  class ArithmeticTuple(object): Line 1297  class ArithmeticTuple(object):
1297     """     """
1298     def __init__(self,*args):     def __init__(self,*args):
1299         """         """
1300         Initializes object with elements C{args}.         Initializes object with elements ``args``.
1301    
1302         @param args: tuple of objects that support inplace add (x+=y) and         :param args: tuple of objects that support inplace add (x+=y) and
1303                      scaling (x=a*y)                      scaling (x=a*y)
1304         """         """
1305         self.__items=list(args)         self.__items=list(args)
# Line 1313  class ArithmeticTuple(object): Line 1308  class ArithmeticTuple(object):
1308         """         """
1309         Returns the number of items.         Returns the number of items.
1310    
1311         @return: number of items         :return: number of items
1312         @rtype: C{int}         :rtype: ``int``
1313         """         """
1314         return len(self.__items)         return len(self.__items)
1315    
# Line 1322  class ArithmeticTuple(object): Line 1317  class ArithmeticTuple(object):
1317         """         """
1318         Returns item at specified position.         Returns item at specified position.
1319    
1320         @param index: index of item to be returned         :param index: index of item to be returned
1321         @type index: C{int}         :type index: ``int``
1322         @return: item with index C{index}         :return: item with index ``index``
1323         """         """
1324         return self.__items.__getitem__(index)         return self.__items.__getitem__(index)
1325    
1326     def __mul__(self,other):     def __mul__(self,other):
1327         """         """
1328         Scales by C{other} from the right.         Scales by ``other`` from the right.
1329    
1330         @param other: scaling factor         :param other: scaling factor
1331         @type other: C{float}         :type other: ``float``
1332         @return: itemwise self*other         :return: itemwise self*other
1333         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1334         """         """
1335         out=[]         out=[]
1336         try:         try:
# Line 1349  class ArithmeticTuple(object): Line 1344  class ArithmeticTuple(object):
1344    
1345     def __rmul__(self,other):     def __rmul__(self,other):
1346         """         """
1347         Scales by C{other} from the left.         Scales by ``other`` from the left.
1348    
1349         @param other: scaling factor         :param other: scaling factor
1350         @type other: C{float}         :type other: ``float``
1351         @return: itemwise other*self         :return: itemwise other*self
1352         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1353         """         """
1354         out=[]         out=[]
1355         try:         try:
# Line 1368  class ArithmeticTuple(object): Line 1363  class ArithmeticTuple(object):
1363    
1364     def __div__(self,other):     def __div__(self,other):
1365         """         """
1366         Scales by (1/C{other}) from the right.         Scales by (1/``other``) from the right.
1367    
1368         @param other: scaling factor         :param other: scaling factor
1369         @type other: C{float}         :type other: ``float``
1370         @return: itemwise self/other         :return: itemwise self/other
1371         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1372         """         """
1373         return self*(1/other)         return self*(1/other)
1374    
1375     def __rdiv__(self,other):     def __rdiv__(self,other):
1376         """         """
1377         Scales by (1/C{other}) from the left.         Scales by (1/``other``) from the left.
1378    
1379         @param other: scaling factor         :param other: scaling factor
1380         @type other: C{float}         :type other: ``float``
1381         @return: itemwise other/self         :return: itemwise other/self
1382         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1383         """         """
1384         out=[]         out=[]
1385         try:         try:
# Line 1398  class ArithmeticTuple(object): Line 1393  class ArithmeticTuple(object):
1393    
1394     def __iadd__(self,other):     def __iadd__(self,other):
1395         """         """
1396         Inplace addition of C{other} to self.         Inplace addition of ``other`` to self.
1397    
1398         @param other: increment         :param other: increment
1399         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1400         """         """
1401         if len(self) != len(other):         if len(self) != len(other):
1402             raise ValueError,"tuple lengths must match."             raise ValueError,"tuple lengths must match."
# Line 1411  class ArithmeticTuple(object): Line 1406  class ArithmeticTuple(object):
1406    
1407     def __add__(self,other):     def __add__(self,other):
1408         """         """
1409         Adds C{other} to self.         Adds ``other`` to self.
1410    
1411         @param other: increment         :param other: increment
1412         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1413         """         """
1414         out=[]         out=[]
1415         try:         try:
# Line 1428  class ArithmeticTuple(object): Line 1423  class ArithmeticTuple(object):
1423    
1424     def __sub__(self,other):     def __sub__(self,other):
1425         """         """
1426         Subtracts C{other} from self.         Subtracts ``other`` from self.
1427    
1428         @param other: decrement         :param other: decrement
1429         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1430         """         """
1431         out=[]         out=[]
1432         try:         try:
# Line 1445  class ArithmeticTuple(object): Line 1440  class ArithmeticTuple(object):
1440    
1441     def __isub__(self,other):     def __isub__(self,other):
1442         """         """
1443         Inplace subtraction of C{other} from self.         Inplace subtraction of ``other`` from self.
1444    
1445         @param other: decrement         :param other: decrement
1446         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1447         """         """
1448         if len(self) != len(other):         if len(self) != len(other):
1449             raise ValueError,"tuple length must match."             raise ValueError,"tuple length must match."
# Line 1471  class HomogeneousSaddlePointProblem(obje Line 1466  class HomogeneousSaddlePointProblem(obje
1466        This class provides a framework for solving linear homogeneous saddle        This class provides a framework for solving linear homogeneous saddle
1467        point problems of the form::        point problems of the form::
1468    
1469            M{Av+B^*p=f}            *Av+B^*p=f*
1470            M{Bv     =0}            *Bv     =0*
1471    
1472        for the unknowns M{v} and M{p} and given operators M{A} and M{B} and        for the unknowns *v* and *p* and given operators *A* and *B* and
1473        given right hand side M{f}. M{B^*} is the adjoint operator of M{B}.        given right hand side *f*. *B^** is the adjoint operator of *B*.
1474        """        """
1475        def __init__(self, adaptSubTolerance=True, **kwargs):        def __init__(self, adaptSubTolerance=True, **kwargs):
1476      """      """
1477      initializes the saddle point problem      initializes the saddle point problem
1478            
1479      @param adaptSubTolerance: If True the tolerance for subproblem is set automatically.      :param adaptSubTolerance: If True the tolerance for subproblem is set automatically.
1480      @type adaptSubTolerance: C{bool}      :type adaptSubTolerance: ``bool``
1481      """      """
1482          self.setTolerance()          self.setTolerance()
1483          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
# Line 1498  class HomogeneousSaddlePointProblem(obje Line 1493  class HomogeneousSaddlePointProblem(obje
1493           """           """
1494           Returns inner product of element p and Bv (overwrite).           Returns inner product of element p and Bv (overwrite).
1495    
1496           @param p: a pressure increment           :param p: a pressure increment
1497           @param v: a residual           :param Bv: a residual
1498           @return: inner product of element p and Bv           :return: inner product of element p and Bv
1499           @rtype: C{float}           :rtype: ``float``
1500           @note: used if PCG is applied.           :note: used if PCG is applied.
1501           """           """
1502           raise NotImplementedError,"no inner product for p and Bv implemented."           raise NotImplementedError,"no inner product for p and Bv implemented."
1503    
# Line 1510  class HomogeneousSaddlePointProblem(obje Line 1505  class HomogeneousSaddlePointProblem(obje
1505           """           """
1506           Returns inner product of p0 and p1 (overwrite).           Returns inner product of p0 and p1 (overwrite).
1507    
1508           @param p0: a pressure           :param p0: a pressure
1509           @param p1: a pressure           :param p1: a pressure
1510           @return: inner product of p0 and p1           :return: inner product of p0 and p1
1511           @rtype: C{float}           :rtype: ``float``
1512           """           """
1513           raise NotImplementedError,"no inner product for p implemented."           raise NotImplementedError,"no inner product for p implemented."
1514        
# Line 1521  class HomogeneousSaddlePointProblem(obje Line 1516  class HomogeneousSaddlePointProblem(obje
1516           """           """
1517           Returns the norm of v (overwrite).           Returns the norm of v (overwrite).
1518    
1519           @param v: a velovity           :param v: a velovity
1520           @return: norm of v           :return: norm of v
1521           @rtype: non-negative C{float}           :rtype: non-negative ``float``
1522           """           """
1523           raise NotImplementedError,"no norm of v implemented."           raise NotImplementedError,"no norm of v implemented."
1524        def getV(self, p, v0):        def getV(self, p, v0):
1525           """           """
1526           return the value for v for a given p (overwrite)           return the value for v for a given p (overwrite)
1527    
1528           @param p: a pressure           :param p: a pressure
1529           @param v0: a initial guess for the value v to return.           :param v0: a initial guess for the value v to return.
1530           @return: v given as M{v= A^{-1} (f-B^*p)}           :return: v given as *v= A^{-1} (f-B^*p)*
1531           """           """
1532           raise NotImplementedError,"no v calculation implemented."           raise NotImplementedError,"no v calculation implemented."
1533    
# Line 1541  class HomogeneousSaddlePointProblem(obje Line 1536  class HomogeneousSaddlePointProblem(obje
1536          """          """
1537          Returns Bv (overwrite).          Returns Bv (overwrite).
1538    
1539          @rtype: equal to the type of p          :rtype: equal to the type of p
1540          @note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1541          """          """
1542          raise NotImplementedError, "no operator B implemented."          raise NotImplementedError, "no operator B implemented."
1543    
# Line 1550  class HomogeneousSaddlePointProblem(obje Line 1545  class HomogeneousSaddlePointProblem(obje
1545          """          """
1546          Returns the norm of Bv (overwrite).          Returns the norm of Bv (overwrite).
1547    
1548          @rtype: equal to the type of p          :rtype: equal to the type of p
1549          @note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1550          """          """
1551          raise NotImplementedError, "no norm of Bv implemented."          raise NotImplementedError, "no norm of Bv implemented."
1552    
1553        def solve_AinvBt(self,p):        def solve_AinvBt(self,p):
1554           """           """
1555           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}           Solves *Av=B^*p* with accuracy `self.getSubProblemTolerance()`
1556           (overwrite).           (overwrite).
1557    
1558           @param p: a pressure increment           :param p: a pressure increment
1559           @return: the solution of M{Av=B^*p}           :return: the solution of *Av=B^*p*
1560           @note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
1561           """           """
1562           raise NotImplementedError,"no operator A implemented."           raise NotImplementedError,"no operator A implemented."
1563    
1564        def solve_prec(self,Bv):        def solve_prec(self,Bv):
1565           """           """
1566           Provides a preconditioner for M{BA^{-1}B^*} applied to Bv with accuracy           Provides a preconditioner for *(BA^{-1}B^ * )* applied to Bv with accuracy
1567           L{self.getSubProblemTolerance()} (overwrite).           `self.getSubProblemTolerance()` (overwrite).
1568    
1569           @rtype: equal to the type of p           :rtype: equal to the type of p
1570           @note: boundary conditions on p should be zero!           :note: boundary conditions on p should be zero!
1571           """           """
1572           raise NotImplementedError,"no preconditioner for Schur complement implemented."           raise NotImplementedError,"no preconditioner for Schur complement implemented."
1573        def setSubProblemTolerance(self):        def setSubProblemTolerance(self):
1574           """           """
1575       Updates the tolerance for subproblems       Updates the tolerance for subproblems
1576       @note: method is typically the method is overwritten.       :note: method is typically the method is overwritten.
1577           """           """
1578           pass           pass
1579        #=============================================================        #=============================================================
# Line 1600  class HomogeneousSaddlePointProblem(obje Line 1595  class HomogeneousSaddlePointProblem(obje
1595        #=============================================================        #=============================================================
1596        def norm_p(self,p):        def norm_p(self,p):
1597            """            """
1598            calculates the norm of C{p}            calculates the norm of ``p``
1599                        
1600            @param p: a pressure            :param p: a pressure
1601            @return: the norm of C{p} using the inner product for pressure            :return: the norm of ``p`` using the inner product for pressure
1602            @rtype: C{float}            :rtype: ``float``
1603            """            """
1604            f=self.inner_p(p,p)            f=self.inner_p(p,p)
1605            if f<0: raise ValueError,"negative pressure norm."            if f<0: raise ValueError,"negative pressure norm."
# Line 1613  class HomogeneousSaddlePointProblem(obje Line 1608  class HomogeneousSaddlePointProblem(obje
1608        """        """
1609        Returns True if tolerance adaption for subproblem is choosen.        Returns True if tolerance adaption for subproblem is choosen.
1610        """        """
1611            self.__adaptSubTolerance            return self.__adaptSubTolerance
1612                
1613        def solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):        def solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1614           """           """
1615           Solves the saddle point problem using initial guesses v and p.           Solves the saddle point problem using initial guesses v and p.
1616    
1617           @param v: initial guess for velocity           :param v: initial guess for velocity
1618           @param p: initial guess for pressure           :param p: initial guess for pressure
1619           @type v: L{Data}           :type v: `Data`
1620           @type p: L{Data}           :type p: `Data`
1621           @param usePCG: indicates the usage of the PCG rather than GMRES scheme.           :param usePCG: indicates the usage of the PCG rather than GMRES scheme.
1622           @param max_iter: maximum number of iteration steps per correction           :param max_iter: maximum number of iteration steps per correction
1623                            attempt                            attempt
1624           @param verbose: if True, shows information on the progress of the           :param verbose: if True, shows information on the progress of the
1625                           saddlepoint problem solver.                           saddlepoint problem solver.
1626           @param iter_restart: restart the iteration after C{iter_restart} steps           :param iter_restart: restart the iteration after ``iter_restart`` steps
1627                                (only used if useUzaw=False)                                (only used if useUzaw=False)
1628           @type usePCG: C{bool}           :type usePCG: ``bool``
1629           @type max_iter: C{int}           :type max_iter: ``int``
1630           @type verbose: C{bool}           :type verbose: ``bool``
1631           @type iter_restart: C{int}           :type iter_restart: ``int``
1632           @rtype: C{tuple} of L{Data} objects           :rtype: ``tuple`` of `Data` objects
1633           """           """
1634           self.verbose=verbose           self.verbose=verbose
1635           rtol=self.getTolerance()           rtol=self.getTolerance()
# Line 1678  class HomogeneousSaddlePointProblem(obje Line 1673  class HomogeneousSaddlePointProblem(obje
1673           """           """
1674           Sets the relative tolerance for (v,p).           Sets the relative tolerance for (v,p).
1675    
1676           @param tolerance: tolerance to be used           :param tolerance: tolerance to be used
1677           @type tolerance: non-negative C{float}           :type tolerance: non-negative ``float``
1678           """           """
1679           if tolerance<0:           if tolerance<0:
1680               raise ValueError,"tolerance must be positive."               raise ValueError,"tolerance must be positive."
# Line 1689  class HomogeneousSaddlePointProblem(obje Line 1684  class HomogeneousSaddlePointProblem(obje
1684           """           """
1685           Returns the relative tolerance.           Returns the relative tolerance.
1686    
1687           @return: relative tolerance           :return: relative tolerance
1688           @rtype: C{float}           :rtype: ``float``
1689           """           """
1690           return self.__rtol           return self.__rtol
1691    
# Line 1698  class HomogeneousSaddlePointProblem(obje Line 1693  class HomogeneousSaddlePointProblem(obje
1693           """           """
1694           Sets the absolute tolerance.           Sets the absolute tolerance.
1695    
1696           @param tolerance: tolerance to be used           :param tolerance: tolerance to be used
1697           @type tolerance: non-negative C{float}           :type tolerance: non-negative ``float``
1698           """           """
1699           if tolerance<0:           if tolerance<0:
1700               raise ValueError,"tolerance must be non-negative."               raise ValueError,"tolerance must be non-negative."
# Line 1709  class HomogeneousSaddlePointProblem(obje Line 1704  class HomogeneousSaddlePointProblem(obje
1704           """           """
1705           Returns the absolute tolerance.           Returns the absolute tolerance.
1706    
1707           @return: absolute tolerance           :return: absolute tolerance
1708           @rtype: C{float}           :rtype: ``float``
1709           """           """
1710           return self.__atol           return self.__atol
1711    
1712        def getSubProblemTolerance(self):        def getSubProblemTolerance(self):
1713           """           """
1714           Sets the relative tolerance to solve the subproblem(s).           Sets the relative tolerance to solve the subproblem(s).
   
          @param rtol: relative tolerence  
          @type rtol: positive C{float}  
1715           """           """
1716           return max(200.*util.EPSILON,self.getTolerance()**2)           return max(200.*util.EPSILON,self.getTolerance()**2)
1717    
# Line 1730  def MaskFromBoundaryTag(domain,*tags): Line 1722  def MaskFromBoundaryTag(domain,*tags):
1722    
1723     Usage: m=MaskFromBoundaryTag(domain, "left", "right")     Usage: m=MaskFromBoundaryTag(domain, "left", "right")
1724    
1725     @param domain: domain to be used     :param domain: domain to be used
1726     @type domain: L{escript.Domain}     :type domain: `escript.Domain`
1727     @param tags: boundary tags     :param tags: boundary tags
1728     @type tags: C{str}     :type tags: ``str``
1729     @return: a mask which marks samples that are touching the boundary tagged     :return: a mask which marks samples that are touching the boundary tagged
1730              by any of the given tags              by any of the given tags
1731     @rtype: L{escript.Data} of rank 0     :rtype: `escript.Data` of rank 0
1732     """     """
1733     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1734     d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))     d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
# Line 1751  def MaskFromTag(domain,*tags): Line 1743  def MaskFromTag(domain,*tags):
1743    
1744     Usage: m=MaskFromTag(domain, "ham")     Usage: m=MaskFromTag(domain, "ham")
1745    
1746     @param domain: domain to be used     :param domain: domain to be used
1747     @type domain: L{escript.Domain}     :type domain: `escript.Domain`
1748     @param tags: boundary tags     :param tags: boundary tags
1749     @type tags: C{str}     :type tags: ``str``
1750     @return: a mask which marks samples that are touching the boundary tagged     :return: a mask which marks samples that are touching the boundary tagged
1751              by any of the given tags              by any of the given tags
1752     @rtype: L{escript.Data} of rank 0     :rtype: `escript.Data` of rank 0
1753     """     """
1754     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1755     d=escript.Scalar(0.,escript.Function(domain))     d=escript.Scalar(0.,escript.Function(domain))

Legend:
Removed from v.2548  
changed lines
  Added in v.2676

  ViewVC Help
Powered by ViewVC 1.1.26