/[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 2719 by gross, Wed Oct 14 06:38:03 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 151  class Projector: Line 146  class Projector:
146      self.__pde.setValue(D = 1.)      self.__pde.setValue(D = 1.)
147      return      return
148    def getSolverOptions(self):    def getSolverOptions(self):
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        return self.__pde.getSolverOptions()
155    
156    def __call__(self, input_data):    def __call__(self, input_data):
157      """      """
158      Projects C{input_data} onto a continuous function.      Projects ``input_data`` onto a continuous function.
159    
160      @param input_data: the data to be projected      :param input_data: the data to be projected
161      """      """
162      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
163      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 192  class NoPDE:
192       """       """
193       Solves the following problem for u:       Solves the following problem for u:
194    
195       M{kronecker[i,j]*D[j]*u[j]=Y[i]}       *kronecker[i,j]*D[j]*u[j]=Y[i]*
196    
197       with constraint       with constraint
198    
199       M{u[j]=r[j]}  where M{q[j]>0}       *u[j]=r[j]*  where *q[j]>0*
200    
201       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.
202    
203       In the case of scalars this takes the form       In the case of scalars this takes the form
204    
205       M{D*u=Y}       *D*u=Y*
206    
207       with constraint       with constraint
208    
209       M{u=r} where M{q>0}       *u=r* where *q>0*
210    
211       where M{D}, M{Y}, M{r} and M{q} are given scalar functions.       where *D*, *Y*, *r* and *q* are given scalar functions.
212    
213       The constraint overwrites any other condition.       The constraint overwrites any other condition.
214    
215       @note: This class is similar to the L{linearPDEs.LinearPDE} class with       :note: This class is similar to the `linearPDEs.LinearPDE` class with
216              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
217              in L{Solution} or L{ReducedSolution}.              in `Solution` or `ReducedSolution`.
218       """       """
219       # 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
220       # this.       # this.
# Line 226  class NoPDE: Line 222  class NoPDE:
222           """           """
223           Initializes the problem.           Initializes the problem.
224    
225           @param domain: domain of the PDE           :param domain: domain of the PDE
226           @type domain: L{Domain}           :type domain: `Domain`
227           @param D: coefficient of the solution           :param D: coefficient of the solution
228           @type D: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type D: ``float``, ``int``, ``numpy.ndarray``, `Data`
229           @param Y: right hand side           :param Y: right hand side
230           @type Y: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type Y: ``float``, ``int``, ``numpy.ndarray``, `Data`
231           @param q: location of constraints           :param q: location of constraints
232           @type q: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type q: ``float``, ``int``, ``numpy.ndarray``, `Data`
233           @param r: value of solution at locations of constraints           :param r: value of solution at locations of constraints
234           @type r: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type r: ``float``, ``int``, ``numpy.ndarray``, `Data`
235           """           """
236           self.__domain=domain           self.__domain=domain
237           self.__D=D           self.__D=D
# Line 247  class NoPDE: Line 243  class NoPDE:
243    
244       def setReducedOn(self):       def setReducedOn(self):
245           """           """
246           Sets the L{FunctionSpace} of the solution to L{ReducedSolution}.           Sets the `FunctionSpace` of the solution to `ReducedSolution`.
247           """           """
248           self.__function_space=escript.ReducedSolution(self.__domain)           self.__function_space=escript.ReducedSolution(self.__domain)
249           self.__u=None           self.__u=None
250    
251       def setReducedOff(self):       def setReducedOff(self):
252           """           """
253           Sets the L{FunctionSpace} of the solution to L{Solution}.           Sets the `FunctionSpace` of the solution to `Solution`.
254           """           """
255           self.__function_space=escript.Solution(self.__domain)           self.__function_space=escript.Solution(self.__domain)
256           self.__u=None           self.__u=None
# Line 263  class NoPDE: Line 259  class NoPDE:
259           """           """
260           Assigns values to the parameters.           Assigns values to the parameters.
261    
262           @param D: coefficient of the solution           :param D: coefficient of the solution
263           @type D: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type D: ``float``, ``int``, ``numpy.ndarray``, `Data`
264           @param Y: right hand side           :param Y: right hand side
265           @type Y: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type Y: ``float``, ``int``, ``numpy.ndarray``, `Data`
266           @param q: location of constraints           :param q: location of constraints
267           @type q: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type q: ``float``, ``int``, ``numpy.ndarray``, `Data`
268           @param r: value of solution at locations of constraints           :param r: value of solution at locations of constraints
269           @type r: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type r: ``float``, ``int``, ``numpy.ndarray``, `Data`
270           """           """
271           if not D==None:           if not D==None:
272              self.__D=D              self.__D=D
# Line 289  class NoPDE: Line 285  class NoPDE:
285           """           """
286           Returns the solution.           Returns the solution.
287    
288           @return: the solution of the problem           :return: the solution of the problem
289           @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or           :rtype: `Data` object in the `FunctionSpace` `Solution` or
290                   L{ReducedSolution}                   `ReducedSolution`
291           """           """
292           if self.__u==None:           if self.__u==None:
293              if self.__D==None:              if self.__D==None:
# Line 324  class Locator: Line 320  class Locator:
320         or FunctionSpace for the sample point which is closest to the given         or FunctionSpace for the sample point which is closest to the given
321         point x.         point x.
322    
323         @param where: function space         :param where: function space
324         @type where: L{escript.FunctionSpace}         :type where: `escript.FunctionSpace`
325         @param x: location(s) of the Locator         :param x: location(s) of the Locator
326         @type x: C{numpy.ndarray} or C{list} of C{numpy.ndarray}         :type x: ``numpy.ndarray`` or ``list`` of ``numpy.ndarray``
327         """         """
328         if isinstance(where,escript.FunctionSpace):         if isinstance(where,escript.FunctionSpace):
329            self.__function_space=where            self.__function_space=where
# Line 354  class Locator: Line 350  class Locator:
350         Returns the coordinates of the Locator as a string.         Returns the coordinates of the Locator as a string.
351         """         """
352         x=self.getX()         x=self.getX()
353         if instance(x,list):         if isinstance(x,list):
354            out="["            out="["
355            first=True            first=True
356            for xx in x:            for xx in x:
# Line 401  class Locator: Line 397  class Locator:
397    
398       def getValue(self,data):       def getValue(self,data):
399          """          """
400          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`
401          object otherwise the object is returned.          object otherwise the object is returned.
402          """          """
403          if isinstance(data,escript.Data):          if isinstance(data,escript.Data):
# Line 467  def PCG(r, Aprod, x, Msolve, bilinearfor Line 463  def PCG(r, Aprod, x, Msolve, bilinearfor
463     """     """
464     Solver for     Solver for
465    
466     M{Ax=b}     *Ax=b*
467    
468     with a symmetric and positive definite operator A (more details required!).     with a symmetric and positive definite operator A (more details required!).
469     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 471  def PCG(r, Aprod, x, Msolve, bilinearfor
471    
472     The iteration is terminated if     The iteration is terminated if
473    
474     M{|r| <= atol+rtol*|r0|}     *|r| <= atol+rtol*|r0|*
475    
476     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
477    
478     M{|r| = sqrt( bilinearform(Msolve(r),r))}     *|r| = sqrt( bilinearform(Msolve(r),r))*
479    
480     For details on the preconditioned conjugate gradient method see the book:     For details on the preconditioned conjugate gradient method see the book:
481    
# Line 487  def PCG(r, Aprod, x, Msolve, bilinearfor Line 483  def PCG(r, Aprod, x, Msolve, bilinearfor
483     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,
484     C. Romine, and H. van der Vorst}.     C. Romine, and H. van der Vorst}.
485    
486     @param r: initial residual M{r=b-Ax}. C{r} is altered.     :param r: initial residual *r=b-Ax*. ``r`` is altered.
487     @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)
488     @param x: an initial guess for the solution     :param x: an initial guess for the solution
489     @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)
490     @param Aprod: returns the value Ax     :param Aprod: returns the value Ax
491     @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
492                  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
493                  like argument C{r}.                  like argument ``r``.
494     @param Msolve: solves Mx=r     :param Msolve: solves Mx=r
495     @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
496                   argument C{r}. The returned object needs to be of the same                   argument ``r``. The returned object needs to be of the same
497                   type like argument C{x}.                   type like argument ``x``.
498     @param bilinearform: inner product C{<x,r>}     :param bilinearform: inner product ``<x,r>``
499     @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
500                         type like argument C{x} and C{r} is. The returned value                         type like argument ``x`` and ``r`` is. The returned value
501                         is a C{float}.                         is a ``float``.
502     @param atol: absolute tolerance     :param atol: absolute tolerance
503     @type atol: non-negative C{float}     :type atol: non-negative ``float``
504     @param rtol: relative tolerance     :param rtol: relative tolerance
505     @type rtol: non-negative C{float}     :type rtol: non-negative ``float``
506     @param iter_max: maximum number of iteration steps     :param iter_max: maximum number of iteration steps
507     @type iter_max: C{int}     :type iter_max: ``int``
508     @return: the solution approximation and the corresponding residual     :return: the solution approximation and the corresponding residual
509     @rtype: C{tuple}     :rtype: ``tuple``
510     @warning: C{r} and C{x} are altered.     :warning: ``r`` and ``x`` are altered.
511     """     """
512     iter=0     iter=0
513     rhat=Msolve(r)     rhat=Msolve(r)
# Line 564  class Defect(object): Line 560  class Defect(object):
560          """          """
561          Returns the inner product of x0 and x1          Returns the inner product of x0 and x1
562    
563          @param x0: value for x0          :param x0: value for x0
564          @param x1: value for x1          :param x1: value for x1
565          @return: the inner product of x0 and x1          :return: the inner product of x0 and x1
566          @rtype: C{float}          :rtype: ``float``
567          """          """
568          return 0          return 0
569    
570      def norm(self,x):      def norm(self,x):
571          """          """
572          Returns the norm of argument C{x}.          Returns the norm of argument ``x``.
573    
574          @param x: a value          :param x: a value
575          @return: norm of argument x          :return: norm of argument x
576          @rtype: C{float}          :rtype: ``float``
577          @note: by default C{sqrt(self.bilinearform(x,x)} is returned.          :note: by default ``sqrt(self.bilinearform(x,x)`` is returned.
578          """          """
579          s=self.bilinearform(x,x)          s=self.bilinearform(x,x)
580          if s<0: raise NegativeNorm,"negative norm."          if s<0: raise NegativeNorm,"negative norm."
# Line 586  class Defect(object): Line 582  class Defect(object):
582    
583      def eval(self,x):      def eval(self,x):
584          """          """
585          Returns the value F of a given C{x}.          Returns the value F of a given ``x``.
586    
587          @param x: value for which the defect C{F} is evaluated          :param x: value for which the defect ``F`` is evaluated
588          @return: value of the defect at C{x}          :return: value of the defect at ``x``
589          """          """
590          return 0          return 0
591    
592      def __call__(self,x):      def __call__(self,x):
593          return self.eval(x)          return self.eval(x)
594    
595      def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):      def setDerivativeIncrementLength(self,inc=1000.*math.sqrt(util.EPSILON)):
596          """          """
597          Sets the relative length of the increment used to approximate the          Sets the relative length of the increment used to approximate the
598          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
599          direction of v with x as a starting point.          direction of v with x as a starting point.
600    
601          @param inc: relative increment length          :param inc: relative increment length
602          @type inc: positive C{float}          :type inc: positive ``float``
603          """          """
604          if inc<=0: raise ValueError,"positive increment required."          if inc<=0: raise ValueError,"positive increment required."
605          self.__inc=inc          self.__inc=inc
# Line 612  class Defect(object): Line 608  class Defect(object):
608          """          """
609          Returns the relative increment length used to approximate the          Returns the relative increment length used to approximate the
610          derivative of the defect.          derivative of the defect.
611          @return: value of the defect at C{x}          :return: value of the defect at ``x``
612          @rtype: positive C{float}          :rtype: positive ``float``
613          """          """
614          return self.__inc          return self.__inc
615    
616      def derivative(self, F0, x0, v, v_is_normalised=True):      def derivative(self, F0, x0, v, v_is_normalised=True):
617          """          """
618          Returns the directional derivative at C{x0} in the direction of C{v}.          Returns the directional derivative at ``x0`` in the direction of ``v``.
619    
620          @param F0: value of this defect at x0          :param F0: value of this defect at x0
621          @param x0: value at which derivative is calculated          :param x0: value at which derivative is calculated
622          @param v: direction          :param v: direction
623          @param v_is_normalised: True to indicate that C{v} is nomalized          :param v_is_normalised: True to indicate that ``v`` is nomalized
624                                  (self.norm(v)=0)                                  (self.norm(v)=0)
625          @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``
626          @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
627                 used but this method maybe overwritten to use exact evaluation.                 used but this method maybe overwritten to use exact evaluation.
628          """          """
629          normx=self.norm(x0)          normx=self.norm(x0)
# Line 645  class Defect(object): Line 641  class Defect(object):
641          return (F1-F0)/epsnew          return (F1-F0)/epsnew
642    
643  ######################################  ######################################
644  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, subtol_max=0.5, gamma=0.9, verbose=False):
645     """     """
646     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
647     criterion:     criterion:
648    
649     M{norm(F(x) <= atol + rtol * norm(F(x0)}     *norm(F(x) <= atol + rtol * norm(F(x0)*
650    
651     where M{x0} is the initial guess.     where *x0* is the initial guess.
652    
653     @param defect: object defining the function M{F}. C{defect.norm} defines the     :param defect: object defining the function *F*. ``defect.norm`` defines the
654                    M{norm} used in the stopping criterion.                    *norm* used in the stopping criterion.
655     @type defect: L{Defect}     :type defect: `Defect`
656     @param x: initial guess for the solution, C{x} is altered.     :param x: initial guess for the solution, ``x`` is altered.
657     @type x: any object type allowing basic operations such as     :type x: any object type allowing basic operations such as
658              C{numpy.ndarray}, L{Data}              ``numpy.ndarray``, `Data`
659     @param iter_max: maximum number of iteration steps     :param iter_max: maximum number of iteration steps
660     @type iter_max: positive C{int}     :type iter_max: positive ``int``
661     @param sub_iter_max: maximum number of inner iteration steps     :param sub_iter_max: maximum number of inner iteration steps
662     @type sub_iter_max: positive C{int}     :type sub_iter_max: positive ``int``
663     @param atol: absolute tolerance for the solution     :param atol: absolute tolerance for the solution
664     @type atol: positive C{float}     :type atol: positive ``float``
665     @param rtol: relative tolerance for the solution     :param rtol: relative tolerance for the solution
666     @type rtol: positive C{float}     :type rtol: positive ``float``
667     @param gamma: tolerance safety factor for inner iteration     :param gamma: tolerance safety factor for inner iteration
668     @type gamma: positive C{float}, less than 1     :type gamma: positive ``float``, less than 1
669     @param sub_tol_max: upper bound for inner tolerance     :param subtol_max: upper bound for inner tolerance
670     @type sub_tol_max: positive C{float}, less than 1     :type subtol_max: positive ``float``, less than 1
671     @return: an approximation of the solution with the desired accuracy     :return: an approximation of the solution with the desired accuracy
672     @rtype: same type as the initial guess     :rtype: same type as the initial guess
673     """     """
674     lmaxit=iter_max     lmaxit=iter_max
675     if atol<0: raise ValueError,"atol needs to be non-negative."     if atol<0: raise ValueError,"atol needs to be non-negative."
676     if rtol<0: raise ValueError,"rtol needs to be non-negative."     if rtol<0: raise ValueError,"rtol needs to be non-negative."
677     if rtol+atol<=0: raise ValueError,"rtol or atol needs to be non-negative."     if rtol+atol<=0: raise ValueError,"rtol or atol needs to be non-negative."
678     if gamma<=0 or gamma>=1: raise ValueError,"tolerance safety factor for inner iteration (gamma =%s) needs to be positive and less than 1."%gamma     if gamma<=0 or gamma>=1: raise ValueError,"tolerance safety factor for inner iteration (gamma =%s) needs to be positive and less than 1."%gamma
679     if sub_tol_max<=0 or sub_tol_max>=1: raise ValueError,"upper bound for inner tolerance for inner iteration (sub_tol_max =%s) needs to be positive and less than 1."%sub_tol_max     if subtol_max<=0 or subtol_max>=1: raise ValueError,"upper bound for inner tolerance for inner iteration (subtol_max =%s) needs to be positive and less than 1."%subtol_max
680    
681     F=defect(x)     F=defect(x)
682     fnrm=defect.norm(F)     fnrm=defect.norm(F)
683     stop_tol=atol + rtol*fnrm     stop_tol=atol + rtol*fnrm
684     sub_tol=sub_tol_max     subtol=subtol_max
685     if verbose: print "NewtonGMRES: initial residual = %e."%fnrm     if verbose: print "NewtonGMRES: initial residual = %e."%fnrm
686     if verbose: print "             tolerance = %e."%sub_tol     if verbose: print "             tolerance = %e."%subtol
687     iter=1     iter=1
688     #     #
689     # main iteration loop     # main iteration loop
# Line 695  def NewtonGMRES(defect, x, iter_max=100, Line 691  def NewtonGMRES(defect, x, iter_max=100,
691     while not fnrm<=stop_tol:     while not fnrm<=stop_tol:
692              if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max              if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
693              #              #
694          #   adjust sub_tol_          #   adjust subtol_
695          #          #
696              if iter > 1:              if iter > 1:
697             rat=fnrm/fnrmo             rat=fnrm/fnrmo
698                 sub_tol_old=sub_tol                 subtol_old=subtol
699             sub_tol=gamma*rat**2             subtol=gamma*rat**2
700             if gamma*sub_tol_old**2 > .1: sub_tol=max(sub_tol,gamma*sub_tol_old**2)             if gamma*subtol_old**2 > .1: subtol=max(subtol,gamma*subtol_old**2)
701             sub_tol=max(min(sub_tol,sub_tol_max), .5*stop_tol/fnrm)             subtol=max(min(subtol,subtol_max), .5*stop_tol/fnrm)
702          #          #
703          # calculate newton increment xc          # calculate newton increment xc
704              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown
# Line 710  def NewtonGMRES(defect, x, iter_max=100, Line 706  def NewtonGMRES(defect, x, iter_max=100,
706              #     if  atol is reached sub_iter returns the numer of steps performed to get there              #     if  atol is reached sub_iter returns the numer of steps performed to get there
707              #              #
708              #              #
709              if verbose: print "             subiteration (GMRES) is called with relative tolerance %e."%sub_tol              if verbose: print "             subiteration (GMRES) is called with relative tolerance %e."%subtol
710              try:              try:
711                 xc, sub_iter=__FDGMRES(F, defect, x, sub_tol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)                 xc, sub_iter=__FDGMRES(F, defect, x, subtol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)
712              except MaxIterReached:              except MaxIterReached:
713                 raise MaxIterReached,"maximum number of %s steps reached."%iter_max                 raise MaxIterReached,"maximum number of %s steps reached."%iter_max
714              if sub_iter<0:              if sub_iter<0:
# Line 731  def NewtonGMRES(defect, x, iter_max=100, Line 727  def NewtonGMRES(defect, x, iter_max=100,
727  def __givapp(c,s,vin):  def __givapp(c,s,vin):
728      """      """
729      Applies a sequence of Givens rotations (c,s) recursively to the vector      Applies a sequence of Givens rotations (c,s) recursively to the vector
730      C{vin}      ``vin``
731    
732      @warning: C{vin} is altered.      :warning: ``vin`` is altered.
733      """      """
734      vrot=vin      vrot=vin
735      if isinstance(c,float):      if isinstance(c,float):
# Line 838  def GMRES(r, Aprod, x, bilinearform, ato Line 834  def GMRES(r, Aprod, x, bilinearform, ato
834     """     """
835     Solver for     Solver for
836    
837     M{Ax=b}     *Ax=b*
838    
839     with a general operator A (more details required!).     with a general operator A (more details required!).
840     It uses the generalized minimum residual method (GMRES).     It uses the generalized minimum residual method (GMRES).
841    
842     The iteration is terminated if     The iteration is terminated if
843    
844     M{|r| <= atol+rtol*|r0|}     *|r| <= atol+rtol*|r0|*
845    
846     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
847    
848     M{|r| = sqrt( bilinearform(r,r))}     *|r| = sqrt( bilinearform(r,r))*
849    
850     @param r: initial residual M{r=b-Ax}. C{r} is altered.     :param r: initial residual *r=b-Ax*. ``r`` is altered.
851     @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)
852     @param x: an initial guess for the solution     :param x: an initial guess for the solution
853     @type x: same like C{r}     :type x: same like ``r``
854     @param Aprod: returns the value Ax     :param Aprod: returns the value Ax
855     @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
856                  argument C{x}. The returned object needs to be of the same                  argument ``x``. The returned object needs to be of the same
857                  type like argument C{r}.                  type like argument ``r``.
858     @param bilinearform: inner product C{<x,r>}     :param bilinearform: inner product ``<x,r>``
859     @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
860                         type like argument C{x} and C{r}. The returned value is                         type like argument ``x`` and ``r``. The returned value is
861                         a C{float}.                         a ``float``.
862     @param atol: absolute tolerance     :param atol: absolute tolerance
863     @type atol: non-negative C{float}     :type atol: non-negative ``float``
864     @param rtol: relative tolerance     :param rtol: relative tolerance
865     @type rtol: non-negative C{float}     :type rtol: non-negative ``float``
866     @param iter_max: maximum number of iteration steps     :param iter_max: maximum number of iteration steps
867     @type iter_max: C{int}     :type iter_max: ``int``
868     @param iter_restart: in order to save memory the orthogonalization process     :param iter_restart: in order to save memory the orthogonalization process
869                          is terminated after C{iter_restart} steps and the                          is terminated after ``iter_restart`` steps and the
870                          iteration is restarted.                          iteration is restarted.
871     @type iter_restart: C{int}     :type iter_restart: ``int``
872     @return: the solution approximation and the corresponding residual     :return: the solution approximation and the corresponding residual
873     @rtype: C{tuple}     :rtype: ``tuple``
874     @warning: C{r} and C{x} are altered.     :warning: ``r`` and ``x`` are altered.
875     """     """
876     m=iter_restart     m=iter_restart
877     restarted=False     restarted=False
# Line 1007  def MINRES(r, Aprod, x, Msolve, bilinear Line 1003  def MINRES(r, Aprod, x, Msolve, bilinear
1003      """      """
1004      Solver for      Solver for
1005    
1006      M{Ax=b}      *Ax=b*
1007    
1008      with a symmetric and positive definite operator A (more details required!).      with a symmetric and positive definite operator A (more details required!).
1009      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 1011  def MINRES(r, Aprod, x, Msolve, bilinear
1011    
1012      The iteration is terminated if      The iteration is terminated if
1013    
1014      M{|r| <= atol+rtol*|r0|}      *|r| <= atol+rtol*|r0|*
1015    
1016      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
1017    
1018      M{|r| = sqrt( bilinearform(Msolve(r),r))}      *|r| = sqrt( bilinearform(Msolve(r),r))*
1019    
1020      For details on the preconditioned conjugate gradient method see the book:      For details on the preconditioned conjugate gradient method see the book:
1021    
# Line 1027  def MINRES(r, Aprod, x, Msolve, bilinear Line 1023  def MINRES(r, Aprod, x, Msolve, bilinear
1023      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,
1024      C. Romine, and H. van der Vorst}.      C. Romine, and H. van der Vorst}.
1025    
1026      @param r: initial residual M{r=b-Ax}. C{r} is altered.      :param r: initial residual *r=b-Ax*. ``r`` is altered.
1027      @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)
1028      @param x: an initial guess for the solution      :param x: an initial guess for the solution
1029      @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)
1030      @param Aprod: returns the value Ax      :param Aprod: returns the value Ax
1031      @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
1032                   argument C{x}. The returned object needs to be of the same                   argument ``x``. The returned object needs to be of the same
1033                   type like argument C{r}.                   type like argument ``r``.
1034      @param Msolve: solves Mx=r      :param Msolve: solves Mx=r
1035      @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
1036                    argument C{r}. The returned object needs to be of the same                    argument ``r``. The returned object needs to be of the same
1037                    type like argument C{x}.                    type like argument ``x``.
1038      @param bilinearform: inner product C{<x,r>}      :param bilinearform: inner product ``<x,r>``
1039      @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
1040                          type like argument C{x} and C{r} is. The returned value                          type like argument ``x`` and ``r`` is. The returned value
1041                          is a C{float}.                          is a ``float``.
1042      @param atol: absolute tolerance      :param atol: absolute tolerance
1043      @type atol: non-negative C{float}      :type atol: non-negative ``float``
1044      @param rtol: relative tolerance      :param rtol: relative tolerance
1045      @type rtol: non-negative C{float}      :type rtol: non-negative ``float``
1046      @param iter_max: maximum number of iteration steps      :param iter_max: maximum number of iteration steps
1047      @type iter_max: C{int}      :type iter_max: ``int``
1048      @return: the solution approximation and the corresponding residual      :return: the solution approximation and the corresponding residual
1049      @rtype: C{tuple}      :rtype: ``tuple``
1050      @warning: C{r} and C{x} are altered.      :warning: ``r`` and ``x`` are altered.
1051      """      """
1052      #------------------------------------------------------------------      #------------------------------------------------------------------
1053      # 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 1181  def TFQMR(r, Aprod, x, bilinearform, ato
1181    """    """
1182    Solver for    Solver for
1183    
1184    M{Ax=b}    *Ax=b*
1185    
1186    with a general operator A (more details required!).    with a general operator A (more details required!).
1187    It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).    It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).
1188    
1189    The iteration is terminated if    The iteration is terminated if
1190    
1191    M{|r| <= atol+rtol*|r0|}    *|r| <= atol+rtol*|r0|*
1192    
1193    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
1194    
1195    M{|r| = sqrt( bilinearform(r,r))}    *|r| = sqrt( bilinearform(r,r))*
1196    
1197    @param r: initial residual M{r=b-Ax}. C{r} is altered.    :param r: initial residual *r=b-Ax*. ``r`` is altered.
1198    @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)
1199    @param x: an initial guess for the solution    :param x: an initial guess for the solution
1200    @type x: same like C{r}    :type x: same like ``r``
1201    @param Aprod: returns the value Ax    :param Aprod: returns the value Ax
1202    @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
1203                 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
1204                 like argument C{r}.                 like argument ``r``.
1205    @param bilinearform: inner product C{<x,r>}    :param bilinearform: inner product ``<x,r>``
1206    @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
1207                        type like argument C{x} and C{r}. The returned value is                        type like argument ``x`` and ``r``. The returned value is
1208                        a C{float}.                        a ``float``.
1209    @param atol: absolute tolerance    :param atol: absolute tolerance
1210    @type atol: non-negative C{float}    :type atol: non-negative ``float``
1211    @param rtol: relative tolerance    :param rtol: relative tolerance
1212    @type rtol: non-negative C{float}    :type rtol: non-negative ``float``
1213    @param iter_max: maximum number of iteration steps    :param iter_max: maximum number of iteration steps
1214    @type iter_max: C{int}    :type iter_max: ``int``
1215    @rtype: C{tuple}    :rtype: ``tuple``
1216    @warning: C{r} and C{x} are altered.    :warning: ``r`` and ``x`` are altered.
1217    """    """
1218    u1=0    u1=0
1219    u2=0    u2=0
# Line 1287  def TFQMR(r, Aprod, x, bilinearform, ato Line 1283  def TFQMR(r, Aprod, x, bilinearform, ato
1283    
1284  class ArithmeticTuple(object):  class ArithmeticTuple(object):
1285     """     """
1286     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
1287     ArithmeticTuple and C{a} is a float.     ArithmeticTuple and ``a`` is a float.
1288    
1289     Example of usage::     Example of usage::
1290    
# Line 1302  class ArithmeticTuple(object): Line 1298  class ArithmeticTuple(object):
1298     """     """
1299     def __init__(self,*args):     def __init__(self,*args):
1300         """         """
1301         Initializes object with elements C{args}.         Initializes object with elements ``args``.
1302    
1303         @param args: tuple of objects that support inplace add (x+=y) and         :param args: tuple of objects that support inplace add (x+=y) and
1304                      scaling (x=a*y)                      scaling (x=a*y)
1305         """         """
1306         self.__items=list(args)         self.__items=list(args)
# Line 1313  class ArithmeticTuple(object): Line 1309  class ArithmeticTuple(object):
1309         """         """
1310         Returns the number of items.         Returns the number of items.
1311    
1312         @return: number of items         :return: number of items
1313         @rtype: C{int}         :rtype: ``int``
1314         """         """
1315         return len(self.__items)         return len(self.__items)
1316    
# Line 1322  class ArithmeticTuple(object): Line 1318  class ArithmeticTuple(object):
1318         """         """
1319         Returns item at specified position.         Returns item at specified position.
1320    
1321         @param index: index of item to be returned         :param index: index of item to be returned
1322         @type index: C{int}         :type index: ``int``
1323         @return: item with index C{index}         :return: item with index ``index``
1324         """         """
1325         return self.__items.__getitem__(index)         return self.__items.__getitem__(index)
1326    
1327     def __mul__(self,other):     def __mul__(self,other):
1328         """         """
1329         Scales by C{other} from the right.         Scales by ``other`` from the right.
1330    
1331         @param other: scaling factor         :param other: scaling factor
1332         @type other: C{float}         :type other: ``float``
1333         @return: itemwise self*other         :return: itemwise self*other
1334         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1335         """         """
1336         out=[]         out=[]
1337         try:         try:
# Line 1349  class ArithmeticTuple(object): Line 1345  class ArithmeticTuple(object):
1345    
1346     def __rmul__(self,other):     def __rmul__(self,other):
1347         """         """
1348         Scales by C{other} from the left.         Scales by ``other`` from the left.
1349    
1350         @param other: scaling factor         :param other: scaling factor
1351         @type other: C{float}         :type other: ``float``
1352         @return: itemwise other*self         :return: itemwise other*self
1353         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1354         """         """
1355         out=[]         out=[]
1356         try:         try:
# Line 1368  class ArithmeticTuple(object): Line 1364  class ArithmeticTuple(object):
1364    
1365     def __div__(self,other):     def __div__(self,other):
1366         """         """
1367         Scales by (1/C{other}) from the right.         Scales by (1/``other``) from the right.
1368    
1369         @param other: scaling factor         :param other: scaling factor
1370         @type other: C{float}         :type other: ``float``
1371         @return: itemwise self/other         :return: itemwise self/other
1372         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1373         """         """
1374         return self*(1/other)         return self*(1/other)
1375    
1376     def __rdiv__(self,other):     def __rdiv__(self,other):
1377         """         """
1378         Scales by (1/C{other}) from the left.         Scales by (1/``other``) from the left.
1379    
1380         @param other: scaling factor         :param other: scaling factor
1381         @type other: C{float}         :type other: ``float``
1382         @return: itemwise other/self         :return: itemwise other/self
1383         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1384         """         """
1385         out=[]         out=[]
1386         try:         try:
# Line 1398  class ArithmeticTuple(object): Line 1394  class ArithmeticTuple(object):
1394    
1395     def __iadd__(self,other):     def __iadd__(self,other):
1396         """         """
1397         Inplace addition of C{other} to self.         Inplace addition of ``other`` to self.
1398    
1399         @param other: increment         :param other: increment
1400         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1401         """         """
1402         if len(self) != len(other):         if len(self) != len(other):
1403             raise ValueError,"tuple lengths must match."             raise ValueError,"tuple lengths must match."
# Line 1411  class ArithmeticTuple(object): Line 1407  class ArithmeticTuple(object):
1407    
1408     def __add__(self,other):     def __add__(self,other):
1409         """         """
1410         Adds C{other} to self.         Adds ``other`` to self.
1411    
1412         @param other: increment         :param other: increment
1413         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1414         """         """
1415         out=[]         out=[]
1416         try:         try:
# Line 1428  class ArithmeticTuple(object): Line 1424  class ArithmeticTuple(object):
1424    
1425     def __sub__(self,other):     def __sub__(self,other):
1426         """         """
1427         Subtracts C{other} from self.         Subtracts ``other`` from self.
1428    
1429         @param other: decrement         :param other: decrement
1430         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1431         """         """
1432         out=[]         out=[]
1433         try:         try:
# Line 1445  class ArithmeticTuple(object): Line 1441  class ArithmeticTuple(object):
1441    
1442     def __isub__(self,other):     def __isub__(self,other):
1443         """         """
1444         Inplace subtraction of C{other} from self.         Inplace subtraction of ``other`` from self.
1445    
1446         @param other: decrement         :param other: decrement
1447         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1448         """         """
1449         if len(self) != len(other):         if len(self) != len(other):
1450             raise ValueError,"tuple length must match."             raise ValueError,"tuple length must match."
# Line 1471  class HomogeneousSaddlePointProblem(obje Line 1467  class HomogeneousSaddlePointProblem(obje
1467        This class provides a framework for solving linear homogeneous saddle        This class provides a framework for solving linear homogeneous saddle
1468        point problems of the form::        point problems of the form::
1469    
1470            M{Av+B^*p=f}            *Av+B^*p=f*
1471            M{Bv     =0}            *Bv     =0*
1472    
1473        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
1474        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*.
1475          *A* may depend weakly on *v* and *p*.
1476        """        """
1477        def __init__(self, adaptSubTolerance=True, **kwargs):        def __init__(self, **kwargs):
1478      """      """
1479      initializes the saddle point problem      initializes the saddle point problem
       
     @param adaptSubTolerance: If True the tolerance for subproblem is set automatically.  
     @type adaptSubTolerance: C{bool}  
1480      """      """
1481            self.resetControlParameters()
1482          self.setTolerance()          self.setTolerance()
1483          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
1484      self.__adaptSubTolerance=adaptSubTolerance        def resetControlParameters(self,gamma=0.85, gamma_min=1.e-8,chi_max=0.1, omega_div=0.2, omega_conv=1.1, rtol_min=1.e-7, rtol_max=0.1, chi=1., C_p=1., C_v=1., safety_factor=0.3):
1485             """
1486             sets a control parameter
1487    
1488             :param gamma: ``1/(1-gamma)`` controls the perturbation of the converegence rate due to termination errors in the subproblems.
1489             :type gamma: ``float``
1490             :param gamma_min: minimum value for ``gamma``.
1491             :type gamma_min: ``float``
1492             :param chi_max: maximum tolerable converegence rate.
1493             :type chi_max: ``float``
1494             :param omega_div: reduction fact for ``gamma`` if no convergence is detected.
1495             :type omega_div: ``float``
1496             :param omega_conv: raise fact for ``gamma`` if convergence is detected.
1497             :type omega_conv: ``float``
1498             :param rtol_min: minimum relative tolerance used to calculate presssure and velocity increment.
1499             :type rtol_min: ``float``
1500             :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1501             :type rtol_max: ``float``
1502             :param chi: initial convergence measure.
1503             :type chi: ``float``
1504             :param C_p: initial value for constant to adjust pressure tolerance
1505             :type C_p: ``float``
1506             :param C_v: initial value for constant to adjust velocity tolerance
1507             :type C_v: ``float``
1508             :param safety_factor: safety factor for addjustment of pressure and velocity tolerance from stopping criteria
1509             :type safety_factor: ``float``
1510             """
1511             self.setControlParameter(gamma, gamma_min ,chi_max , omega_div , omega_conv, rtol_min , rtol_max, chi,C_p, C_v,safety_factor)
1512    
1513          def setControlParameter(self,gamma=None, gamma_min=None ,chi_max=None, omega_div=None, omega_conv=None, rtol_min=None, rtol_max=None, chi=None, C_p=None, C_v=None, safety_factor=None):
1514             """
1515             sets a control parameter
1516    
1517             :param gamma: ``1/(1-gamma)`` controls the perturbation of the converegence rate due to termination errors in the subproblems.
1518             :type gamma: ``float``
1519             :param gamma_min: minimum value for ``gamma``.
1520             :type gamma_min: ``float``
1521             :param chi_max: maximum tolerable converegence rate.
1522             :type chi_max: ``float``
1523             :param omega_div: reduction fact for ``gamma`` if no convergence is detected.
1524             :type omega_div: ``float``
1525             :param omega_conv: raise fact for ``gamma`` if convergence is detected.
1526             :type omega_conv: ``float``
1527             :param rtol_min: minimum relative tolerance used to calculate presssure and velocity increment.
1528             :type rtol_min: ``float``
1529             :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1530             :type rtol_max: ``float``
1531             :param chi: initial convergence measure.
1532             :type chi: ``float``
1533             :param C_p: initial value for constant to adjust pressure tolerance
1534             :type C_p: ``float``
1535             :param C_v: initial value for constant to adjust velocity tolerance
1536             :type C_v: ``float``
1537             :param safety_factor: safety factor for addjustment of pressure and velocity tolerance from stopping criteria
1538             :type safety_factor: ``float``
1539             """
1540             if not gamma == None:
1541                if gamma<=0 or gamma>=1:
1542                   raise ValueError,"Initial gamma needs to be positive and less than 1."
1543             else:
1544                gamma=self.__gamma
1545    
1546             if not gamma_min == None:
1547                if gamma_min<=0 or gamma_min>=1:
1548                   raise ValueError,"gamma_min needs to be positive and less than 1."
1549             else:
1550                gamma_min = self.__gamma_min
1551    
1552             if not chi_max == None:
1553                if chi_max<=0 or chi_max>=1:
1554                   raise ValueError,"chi_max needs to be positive and less than 1."
1555             else:
1556                chi_max = self.__chi_max
1557    
1558             if not omega_div == None:
1559                if omega_div<=0 or omega_div >=1:
1560                   raise ValueError,"omega_div needs to be positive and less than 1."
1561             else:
1562                omega_div=self.__omega_div
1563    
1564             if not omega_conv == None:
1565                if omega_conv<1:
1566                   raise ValueError,"omega_conv needs to be greater or equal to 1."
1567             else:
1568                omega_conv=self.__omega_conv
1569    
1570             if not rtol_min == None:
1571                if rtol_min<=0 or rtol_min>=1:
1572                   raise ValueError,"rtol_min needs to be positive and less than 1."
1573             else:
1574                rtol_min=self.__rtol_min
1575    
1576             if not rtol_max == None:
1577                if rtol_max<=0 or rtol_max>=1:
1578                   raise ValueError,"rtol_max needs to be positive and less than 1."
1579             else:
1580                rtol_max=self.__rtol_max
1581    
1582             if not chi == None:
1583                if chi<=0 or chi>1:
1584                   raise ValueError,"chi needs to be positive and less than 1."
1585             else:
1586                chi=self.__chi
1587    
1588             if not C_p == None:
1589                if C_p<1:
1590                   raise ValueError,"C_p need to be greater or equal to 1."
1591             else:
1592                C_p=self.__C_p
1593    
1594             if not C_v == None:
1595                if C_v<1:
1596                   raise ValueError,"C_v need to be greater or equal to 1."
1597             else:
1598                C_v=self.__C_v
1599    
1600             if not safety_factor == None:
1601                if safety_factor<=0 or safety_factor>1:
1602                   raise ValueError,"safety_factor need to be between zero and one."
1603             else:
1604                safety_factor=self.__safety_factor
1605    
1606             if gamma<gamma_min:
1607                   raise ValueError,"gamma = %e needs to be greater or equal gamma_min = %e."%(gamma,gamma_min)
1608             if rtol_max<=rtol_min:
1609                   raise ValueError,"rtol_max = %e needs to be greater rtol_min = %e."%(rtol_max,rtol_min)
1610                
1611             self.__gamma = gamma
1612             self.__gamma_min = gamma_min
1613             self.__chi_max = chi_max
1614             self.__omega_div = omega_div
1615             self.__omega_conv = omega_conv
1616             self.__rtol_min = rtol_min
1617             self.__rtol_max = rtol_max
1618             self.__chi = chi
1619             self.__C_p = C_p
1620             self.__C_v = C_v
1621             self.__safety_factor = safety_factor
1622    
1623        #=============================================================        #=============================================================
1624        def initialize(self):        def initialize(self):
1625          """          """
# Line 1498  class HomogeneousSaddlePointProblem(obje Line 1631  class HomogeneousSaddlePointProblem(obje
1631           """           """
1632           Returns inner product of element p and Bv (overwrite).           Returns inner product of element p and Bv (overwrite).
1633    
1634           @param p: a pressure increment           :param p: a pressure increment
1635           @param v: a residual           :param Bv: a residual
1636           @return: inner product of element p and Bv           :return: inner product of element p and Bv
1637           @rtype: C{float}           :rtype: ``float``
1638           @note: used if PCG is applied.           :note: used if PCG is applied.
1639           """           """
1640           raise NotImplementedError,"no inner product for p and Bv implemented."           raise NotImplementedError,"no inner product for p and Bv implemented."
1641    
# Line 1510  class HomogeneousSaddlePointProblem(obje Line 1643  class HomogeneousSaddlePointProblem(obje
1643           """           """
1644           Returns inner product of p0 and p1 (overwrite).           Returns inner product of p0 and p1 (overwrite).
1645    
1646           @param p0: a pressure           :param p0: a pressure
1647           @param p1: a pressure           :param p1: a pressure
1648           @return: inner product of p0 and p1           :return: inner product of p0 and p1
1649           @rtype: C{float}           :rtype: ``float``
1650           """           """
1651           raise NotImplementedError,"no inner product for p implemented."           raise NotImplementedError,"no inner product for p implemented."
1652        
# Line 1521  class HomogeneousSaddlePointProblem(obje Line 1654  class HomogeneousSaddlePointProblem(obje
1654           """           """
1655           Returns the norm of v (overwrite).           Returns the norm of v (overwrite).
1656    
1657           @param v: a velovity           :param v: a velovity
1658           @return: norm of v           :return: norm of v
1659           @rtype: non-negative C{float}           :rtype: non-negative ``float``
1660           """           """
1661           raise NotImplementedError,"no norm of v implemented."           raise NotImplementedError,"no norm of v implemented."
1662        def getV(self, p, v0):        def getDV(self, p, v, tol):
1663           """           """
1664           return the value for v for a given p (overwrite)           return a correction to the value for a given v and a given p with accuracy `tol` (overwrite)
1665    
1666           @param p: a pressure           :param p: pressure
1667           @param v0: a initial guess for the value v to return.           :param v: pressure
1668           @return: v given as M{v= A^{-1} (f-B^*p)}           :return: dv given as *dv= A^{-1} (f-A v-B^*p)*
1669             :note: Only *A* may depend on *v* and *p*
1670           """           """
1671           raise NotImplementedError,"no v calculation implemented."           raise NotImplementedError,"no dv calculation implemented."
1672    
1673                    
1674        def Bv(self,v):        def Bv(self,v, tol):
1675          """          """
1676          Returns Bv (overwrite).          Returns Bv with accuracy `tol` (overwrite)
1677    
1678          @rtype: equal to the type of p          :rtype: equal to the type of p
1679          @note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1680          """          """
1681          raise NotImplementedError, "no operator B implemented."          raise NotImplementedError, "no operator B implemented."
1682    
# Line 1550  class HomogeneousSaddlePointProblem(obje Line 1684  class HomogeneousSaddlePointProblem(obje
1684          """          """
1685          Returns the norm of Bv (overwrite).          Returns the norm of Bv (overwrite).
1686    
1687          @rtype: equal to the type of p          :rtype: equal to the type of p
1688          @note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1689          """          """
1690          raise NotImplementedError, "no norm of Bv implemented."          raise NotImplementedError, "no norm of Bv implemented."
1691    
1692        def solve_AinvBt(self,p):        def solve_AinvBt(self,dp, tol):
1693           """           """
1694           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}           Solves *A dv=B^*dp* with accuracy `tol`
          (overwrite).  
1695    
1696           @param p: a pressure increment           :param dp: a pressure increment
1697           @return: the solution of M{Av=B^*p}           :return: the solution of *A dv=B^*dp*
1698           @note: boundary conditions on v should be zero!           :note: boundary conditions on dv should be zero! *A* is the operator used in ``getDV`` and must not be altered.
1699           """           """
1700           raise NotImplementedError,"no operator A implemented."           raise NotImplementedError,"no operator A implemented."
1701    
1702        def solve_prec(self,Bv):        def solve_prec(self,Bv, tol):
1703           """           """
1704           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 `tol`
          L{self.getSubProblemTolerance()} (overwrite).  
1705    
1706           @rtype: equal to the type of p           :rtype: equal to the type of p
1707           @note: boundary conditions on p should be zero!           :note: boundary conditions on p should be zero!
1708           """           """
1709           raise NotImplementedError,"no preconditioner for Schur complement implemented."           raise NotImplementedError,"no preconditioner for Schur complement implemented."
       def setSubProblemTolerance(self):  
          """  
      Updates the tolerance for subproblems  
      @note: method is typically the method is overwritten.  
          """  
          pass  
1710        #=============================================================        #=============================================================
1711        def __Aprod_PCG(self,p):        def __Aprod_PCG(self,dp):
1712            dv=self.solve_AinvBt(p)            dv=self.solve_AinvBt(dp, self.__subtol)
1713            return ArithmeticTuple(dv,self.Bv(dv))            return ArithmeticTuple(dv,self.Bv(dv, self.__subtol))
1714    
1715        def __inner_PCG(self,p,r):        def __inner_PCG(self,p,r):
1716           return self.inner_pBv(p,r[1])           return self.inner_pBv(p,r[1])
1717    
1718        def __Msolve_PCG(self,r):        def __Msolve_PCG(self,r):
1719            return self.solve_prec(r[1])            return self.solve_prec(r[1], self.__subtol)
1720        #=============================================================        #=============================================================
1721        def __Aprod_GMRES(self,p):        def __Aprod_GMRES(self,p):
1722            return self.solve_prec(self.Bv(self.solve_AinvBt(p)))            return self.solve_prec(self.Bv(self.solve_AinvBt(p, self.__subtol), self.__subtol), self.__subtol)
1723        def __inner_GMRES(self,p0,p1):        def __inner_GMRES(self,p0,p1):
1724           return self.inner_p(p0,p1)           return self.inner_p(p0,p1)
1725    
1726        #=============================================================        #=============================================================
1727        def norm_p(self,p):        def norm_p(self,p):
1728            """            """
1729            calculates the norm of C{p}            calculates the norm of ``p``
1730                        
1731            @param p: a pressure            :param p: a pressure
1732            @return: the norm of C{p} using the inner product for pressure            :return: the norm of ``p`` using the inner product for pressure
1733            @rtype: C{float}            :rtype: ``float``
1734            """            """
1735            f=self.inner_p(p,p)            f=self.inner_p(p,p)
1736            if f<0: raise ValueError,"negative pressure norm."            if f<0: raise ValueError,"negative pressure norm."
1737            return math.sqrt(f)            return math.sqrt(f)
       def adaptSubTolerance(self):  
       """  
       Returns True if tolerance adaption for subproblem is choosen.  
       """  
           self.__adaptSubTolerance  
1738                
1739        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):
1740           """           """
1741           Solves the saddle point problem using initial guesses v and p.           Solves the saddle point problem using initial guesses v and p.
1742    
1743           @param v: initial guess for velocity           :param v: initial guess for velocity
1744           @param p: initial guess for pressure           :param p: initial guess for pressure
1745           @type v: L{Data}           :type v: `Data`
1746           @type p: L{Data}           :type p: `Data`
1747           @param usePCG: indicates the usage of the PCG rather than GMRES scheme.           :param usePCG: indicates the usage of the PCG rather than GMRES scheme.
1748           @param max_iter: maximum number of iteration steps per correction           :param max_iter: maximum number of iteration steps per correction
1749                            attempt                            attempt
1750           @param verbose: if True, shows information on the progress of the           :param verbose: if True, shows information on the progress of the
1751                           saddlepoint problem solver.                           saddlepoint problem solver.
1752           @param iter_restart: restart the iteration after C{iter_restart} steps           :param iter_restart: restart the iteration after ``iter_restart`` steps
1753                                (only used if useUzaw=False)                                (only used if useUzaw=False)
1754           @type usePCG: C{bool}           :type usePCG: ``bool``
1755           @type max_iter: C{int}           :type max_iter: ``int``
1756           @type verbose: C{bool}           :type verbose: ``bool``
1757           @type iter_restart: C{int}           :type iter_restart: ``int``
1758           @rtype: C{tuple} of L{Data} objects           :rtype: ``tuple`` of `Data` objects
1759           """           """
1760           self.verbose=verbose           self.verbose=verbose
1761           rtol=self.getTolerance()           rtol=self.getTolerance()
1762           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
      if self.adaptSubTolerance(): self.setSubProblemTolerance()  
1763           correction_step=0           correction_step=0
1764           converged=False           converged=False
1765             error=None
1766             chi=None
1767             gamma=self.__gamma
1768             C_p=self.__C_p
1769             C_v=self.__C_v
1770           while not converged:           while not converged:
1771                  if error== None or chi == None:
1772                      gamma_new=gamma/self.__omega_conv
1773                  else:
1774                     if chi < self.__chi_max:
1775                        gamma_new=min(max(gamma*self.__omega_conv,1-chi*error/(self.__safety_factor*ATOL)), 1-chi/self.__chi_max)
1776                     else:
1777                        gamma_new=gamma*self.__omega_div
1778                  gamma=max(gamma_new, self.__gamma_min)
1779                # calculate velocity for current pressure:                # calculate velocity for current pressure:
1780                v=self.getV(p,v)                rtol_v=min(max(gamma/(1.+gamma)/C_v,self.__rtol_min), self.__rtol_max)
1781                Bv=self.Bv(v)                rtol_p=min(max(gamma/C_p, self.__rtol_min), self.__rtol_max)
1782                norm_v=self.norm_v(v)                self.__subtol=rtol_p**2
1783                norm_Bv=self.norm_Bv(Bv)                if self.verbose: print "HomogeneousSaddlePointProblem: step %s: gamma = %e, rtol_v= %e, rtol_p=%e"%(correction_step,gamma,rtol_v,rtol_p)
1784                ATOL=norm_v*rtol+atol                if self.verbose: print "HomogeneousSaddlePointProblem: subtolerance: %e"%self.__subtol
1785                if self.verbose: print "HomogeneousSaddlePointProblem: norm v= %e, norm_Bv= %e, tolerance = %e."%(norm_v, norm_Bv,ATOL)                # calculate velocity for current pressure: A*dv= F-A*v-B*p
1786                  dv1=self.getDV(p,v,rtol_v)
1787                  v1=v+dv1
1788                  Bv1=self.Bv(v1, self.__subtol)
1789                  norm_Bv1=self.norm_Bv(Bv1)
1790                  norm_dv1=self.norm_v(dv1)
1791                  norm_v1=self.norm_v(v1)
1792                  ATOL=norm_v1*rtol+atol
1793                  if self.verbose: print "HomogeneousSaddlePointProblem: step %s: Bv = %e, dv = %e, v=%e"%(correction_step,norm_Bv1, norm_dv1, norm_v1)
1794                if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."                if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."
1795                if norm_Bv <= ATOL:                if max(norm_Bv1, norm_dv1) <= ATOL:
1796                   converged=True                    converged=True
1797                      v=v1
1798                else:                else:
1799                   correction_step+=1                    # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1
1800                   if correction_step>max_correction_steps:                    if usePCG:
1801                        raise CorrectionFailed,"Given up after %d correction steps."%correction_step                      dp,r,a_norm=PCG(ArithmeticTuple(v1,Bv1),self.__Aprod_PCG,0*p,self.__Msolve_PCG,self.__inner_PCG,atol=0, rtol=rtol_p,iter_max=max_iter, verbose=self.verbose)
1802                   dp=self.solve_prec(Bv)                      v2=r[0]
1803                   if usePCG:                      Bv2=r[1]
1804                     norm2=self.inner_pBv(dp,Bv)                    else:
1805                     if norm2<0: raise ValueError,"negative PCG norm."                      dp=GMRES(self.solve_prec(Bv1,self.__subtol),self.__Aprod_GMRES, 0*p, self.__inner_GMRES,atol=0, rtol=rtol_p,iter_max=max_iter, iter_restart=iter_restart, verbose=self.verbose)
1806                     norm2=math.sqrt(norm2)                      dv2=self.solve_AinvBt(dp, self.__subtol)
1807                   else:                      v2=v1-dv2
1808                     norm2=self.norm_p(dp)                      Bv2=self.Bv(v2, self.__subtol)
1809                   ATOL_ITER=ATOL/norm_Bv*norm2*0.5                    #
1810                   if self.verbose: print "HomogeneousSaddlePointProblem: tolerance for solver: %e"%ATOL_ITER                    # convergence indicators:
1811                   if usePCG:                    #
1812                         p,v0,a_norm=PCG(ArithmeticTuple(v,Bv),self.__Aprod_PCG,p,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL_ITER, rtol=0.,iter_max=max_iter, verbose=self.verbose)                    norm_v2=self.norm_v(v2)
1813                   else:                    norm_dv2=self.norm_v(v2-v)
1814                         p=GMRES(dp,self.__Aprod_GMRES, p, self.__inner_GMRES,atol=ATOL_ITER, rtol=0.,iter_max=max_iter, iter_restart=iter_restart, verbose=self.verbose)                    error_new=max(norm_dv2, norm_Bv1)
1815           if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached."                    norm_Bv2=self.norm_Bv(Bv2)
1816                      if self.verbose: print "HomogeneousSaddlePointProblem: step %s (part 2): Bv = %e, dv = %e, v=%e"%(correction_step,norm_Bv2, norm_dv2, norm_v2)
1817                      if error !=None:
1818                          chi_new=error_new/error
1819                          if self.verbose: print "HomogeneousSaddlePointProblem: step %s: convergence rate = %e"%(correction_step,chi_new)
1820                          if chi != None:
1821                              gamma0=max(gamma, 1-chi/chi_new)
1822                              C_p*=gamma0/gamma
1823                              C_v*=gamma0/gamma*(1+gamma)/(1+gamma0)
1824                          chi=chi_new
1825                      error = error_new
1826                      correction_step+=1
1827                      if correction_step>max_correction_steps:
1828                            raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1829                      v,p=v2,p+dp
1830             if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step
1831       return v,p       return v,p
1832    
1833        #========================================================================        #========================================================================
# Line 1678  class HomogeneousSaddlePointProblem(obje Line 1835  class HomogeneousSaddlePointProblem(obje
1835           """           """
1836           Sets the relative tolerance for (v,p).           Sets the relative tolerance for (v,p).
1837    
1838           @param tolerance: tolerance to be used           :param tolerance: tolerance to be used
1839           @type tolerance: non-negative C{float}           :type tolerance: non-negative ``float``
1840           """           """
1841           if tolerance<0:           if tolerance<0:
1842               raise ValueError,"tolerance must be positive."               raise ValueError,"tolerance must be positive."
# Line 1689  class HomogeneousSaddlePointProblem(obje Line 1846  class HomogeneousSaddlePointProblem(obje
1846           """           """
1847           Returns the relative tolerance.           Returns the relative tolerance.
1848    
1849           @return: relative tolerance           :return: relative tolerance
1850           @rtype: C{float}           :rtype: ``float``
1851           """           """
1852           return self.__rtol           return self.__rtol
1853    
# Line 1698  class HomogeneousSaddlePointProblem(obje Line 1855  class HomogeneousSaddlePointProblem(obje
1855           """           """
1856           Sets the absolute tolerance.           Sets the absolute tolerance.
1857    
1858           @param tolerance: tolerance to be used           :param tolerance: tolerance to be used
1859           @type tolerance: non-negative C{float}           :type tolerance: non-negative ``float``
1860           """           """
1861           if tolerance<0:           if tolerance<0:
1862               raise ValueError,"tolerance must be non-negative."               raise ValueError,"tolerance must be non-negative."
# Line 1709  class HomogeneousSaddlePointProblem(obje Line 1866  class HomogeneousSaddlePointProblem(obje
1866           """           """
1867           Returns the absolute tolerance.           Returns the absolute tolerance.
1868    
1869           @return: absolute tolerance           :return: absolute tolerance
1870           @rtype: C{float}           :rtype: ``float``
1871           """           """
1872           return self.__atol           return self.__atol
1873    
       def getSubProblemTolerance(self):  
          """  
          Sets the relative tolerance to solve the subproblem(s).  
   
          @param rtol: relative tolerence  
          @type rtol: positive C{float}  
          """  
          return max(200.*util.EPSILON,self.getTolerance()**2)  
1874    
1875  def MaskFromBoundaryTag(domain,*tags):  def MaskFromBoundaryTag(domain,*tags):
1876     """     """
# Line 1730  def MaskFromBoundaryTag(domain,*tags): Line 1879  def MaskFromBoundaryTag(domain,*tags):
1879    
1880     Usage: m=MaskFromBoundaryTag(domain, "left", "right")     Usage: m=MaskFromBoundaryTag(domain, "left", "right")
1881    
1882     @param domain: domain to be used     :param domain: domain to be used
1883     @type domain: L{escript.Domain}     :type domain: `escript.Domain`
1884     @param tags: boundary tags     :param tags: boundary tags
1885     @type tags: C{str}     :type tags: ``str``
1886     @return: a mask which marks samples that are touching the boundary tagged     :return: a mask which marks samples that are touching the boundary tagged
1887              by any of the given tags              by any of the given tags
1888     @rtype: L{escript.Data} of rank 0     :rtype: `escript.Data` of rank 0
1889     """     """
1890     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1891     d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))     d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
# Line 1751  def MaskFromTag(domain,*tags): Line 1900  def MaskFromTag(domain,*tags):
1900    
1901     Usage: m=MaskFromTag(domain, "ham")     Usage: m=MaskFromTag(domain, "ham")
1902    
1903     @param domain: domain to be used     :param domain: domain to be used
1904     @type domain: L{escript.Domain}     :type domain: `escript.Domain`
1905     @param tags: boundary tags     :param tags: boundary tags
1906     @type tags: C{str}     :type tags: ``str``
1907     @return: a mask which marks samples that are touching the boundary tagged     :return: a mask which marks samples that are touching the boundary tagged
1908              by any of the given tags              by any of the given tags
1909     @rtype: L{escript.Data} of rank 0     :rtype: `escript.Data` of rank 0
1910     """     """
1911     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1912     d=escript.Scalar(0.,escript.Function(domain))     d=escript.Scalar(0.,escript.Function(domain))

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

  ViewVC Help
Powered by ViewVC 1.1.26