/[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 2745 by jfenwick, Tue Nov 17 04:23:02 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 426  class Locator: Line 422  class Locator:
422          else:          else:
423             return data             return data
424    
425    
426    def getInfLocator(arg):
427        """
428        Return a Locator for a point with the inf value over all arg.
429        """
430        if not isinstance(arg, escript.Data):
431        raise TypeError,"getInfLocator: Unknown argument type."
432        a_inf=util.inf(arg)
433        loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords
434        x=arg.getFunctionSpace().getX()
435        x_min=x.getTupleForGlobalDataPoint(*loc)
436        return Locator(arg.getFunctionSpace(),x_min)
437    
438    def getSupLocator(arg):
439        """
440        Return a Locator for a point with the sup value over all arg.
441        """
442        if not isinstance(arg, escript.Data):
443        raise TypeError,"getInfLocator: Unknown argument type."
444        a_inf=util.sup(arg)
445        loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords
446        x=arg.getFunctionSpace().getX()
447        x_min=x.getTupleForGlobalDataPoint(*loc)
448        return Locator(arg.getFunctionSpace(),x_min)
449        
450    
451  class SolverSchemeException(Exception):  class SolverSchemeException(Exception):
452     """     """
453     This is a generic exception thrown by solvers.     This is a generic exception thrown by solvers.
# Line 467  def PCG(r, Aprod, x, Msolve, bilinearfor Line 489  def PCG(r, Aprod, x, Msolve, bilinearfor
489     """     """
490     Solver for     Solver for
491    
492     M{Ax=b}     *Ax=b*
493    
494     with a symmetric and positive definite operator A (more details required!).     with a symmetric and positive definite operator A (more details required!).
495     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 497  def PCG(r, Aprod, x, Msolve, bilinearfor
497    
498     The iteration is terminated if     The iteration is terminated if
499    
500     M{|r| <= atol+rtol*|r0|}     *|r| <= atol+rtol*|r0|*
501    
502     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
503    
504     M{|r| = sqrt( bilinearform(Msolve(r),r))}     *|r| = sqrt( bilinearform(Msolve(r),r))*
505    
506     For details on the preconditioned conjugate gradient method see the book:     For details on the preconditioned conjugate gradient method see the book:
507    
# Line 487  def PCG(r, Aprod, x, Msolve, bilinearfor Line 509  def PCG(r, Aprod, x, Msolve, bilinearfor
509     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,
510     C. Romine, and H. van der Vorst}.     C. Romine, and H. van der Vorst}.
511    
512     @param r: initial residual M{r=b-Ax}. C{r} is altered.     :param r: initial residual *r=b-Ax*. ``r`` is altered.
513     @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)
514     @param x: an initial guess for the solution     :param x: an initial guess for the solution
515     @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)
516     @param Aprod: returns the value Ax     :param Aprod: returns the value Ax
517     @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
518                  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
519                  like argument C{r}.                  like argument ``r``.
520     @param Msolve: solves Mx=r     :param Msolve: solves Mx=r
521     @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
522                   argument C{r}. The returned object needs to be of the same                   argument ``r``. The returned object needs to be of the same
523                   type like argument C{x}.                   type like argument ``x``.
524     @param bilinearform: inner product C{<x,r>}     :param bilinearform: inner product ``<x,r>``
525     @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
526                         type like argument C{x} and C{r} is. The returned value                         type like argument ``x`` and ``r`` is. The returned value
527                         is a C{float}.                         is a ``float``.
528     @param atol: absolute tolerance     :param atol: absolute tolerance
529     @type atol: non-negative C{float}     :type atol: non-negative ``float``
530     @param rtol: relative tolerance     :param rtol: relative tolerance
531     @type rtol: non-negative C{float}     :type rtol: non-negative ``float``
532     @param iter_max: maximum number of iteration steps     :param iter_max: maximum number of iteration steps
533     @type iter_max: C{int}     :type iter_max: ``int``
534     @return: the solution approximation and the corresponding residual     :return: the solution approximation and the corresponding residual
535     @rtype: C{tuple}     :rtype: ``tuple``
536     @warning: C{r} and C{x} are altered.     :warning: ``r`` and ``x`` are altered.
537     """     """
538     iter=0     iter=0
539     rhat=Msolve(r)     rhat=Msolve(r)
# Line 564  class Defect(object): Line 586  class Defect(object):
586          """          """
587          Returns the inner product of x0 and x1          Returns the inner product of x0 and x1
588    
589          @param x0: value for x0          :param x0: value for x0
590          @param x1: value for x1          :param x1: value for x1
591          @return: the inner product of x0 and x1          :return: the inner product of x0 and x1
592          @rtype: C{float}          :rtype: ``float``
593          """          """
594          return 0          return 0
595    
596      def norm(self,x):      def norm(self,x):
597          """          """
598          Returns the norm of argument C{x}.          Returns the norm of argument ``x``.
599    
600          @param x: a value          :param x: a value
601          @return: norm of argument x          :return: norm of argument x
602          @rtype: C{float}          :rtype: ``float``
603          @note: by default C{sqrt(self.bilinearform(x,x)} is returned.          :note: by default ``sqrt(self.bilinearform(x,x)`` is returned.
604          """          """
605          s=self.bilinearform(x,x)          s=self.bilinearform(x,x)
606          if s<0: raise NegativeNorm,"negative norm."          if s<0: raise NegativeNorm,"negative norm."
# Line 586  class Defect(object): Line 608  class Defect(object):
608    
609      def eval(self,x):      def eval(self,x):
610          """          """
611          Returns the value F of a given C{x}.          Returns the value F of a given ``x``.
612    
613          @param x: value for which the defect C{F} is evaluated          :param x: value for which the defect ``F`` is evaluated
614          @return: value of the defect at C{x}          :return: value of the defect at ``x``
615          """          """
616          return 0          return 0
617    
618      def __call__(self,x):      def __call__(self,x):
619          return self.eval(x)          return self.eval(x)
620    
621      def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):      def setDerivativeIncrementLength(self,inc=1000.*math.sqrt(util.EPSILON)):
622          """          """
623          Sets the relative length of the increment used to approximate the          Sets the relative length of the increment used to approximate the
624          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
625          direction of v with x as a starting point.          direction of v with x as a starting point.
626    
627          @param inc: relative increment length          :param inc: relative increment length
628          @type inc: positive C{float}          :type inc: positive ``float``
629          """          """
630          if inc<=0: raise ValueError,"positive increment required."          if inc<=0: raise ValueError,"positive increment required."
631          self.__inc=inc          self.__inc=inc
# Line 612  class Defect(object): Line 634  class Defect(object):
634          """          """
635          Returns the relative increment length used to approximate the          Returns the relative increment length used to approximate the
636          derivative of the defect.          derivative of the defect.
637          @return: value of the defect at C{x}          :return: value of the defect at ``x``
638          @rtype: positive C{float}          :rtype: positive ``float``
639          """          """
640          return self.__inc          return self.__inc
641    
642      def derivative(self, F0, x0, v, v_is_normalised=True):      def derivative(self, F0, x0, v, v_is_normalised=True):
643          """          """
644          Returns the directional derivative at C{x0} in the direction of C{v}.          Returns the directional derivative at ``x0`` in the direction of ``v``.
645    
646          @param F0: value of this defect at x0          :param F0: value of this defect at x0
647          @param x0: value at which derivative is calculated          :param x0: value at which derivative is calculated
648          @param v: direction          :param v: direction
649          @param v_is_normalised: True to indicate that C{v} is nomalized          :param v_is_normalised: True to indicate that ``v`` is nomalized
650                                  (self.norm(v)=0)                                  (self.norm(v)=0)
651          @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``
652          @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
653                 used but this method maybe overwritten to use exact evaluation.                 used but this method maybe overwritten to use exact evaluation.
654          """          """
655          normx=self.norm(x0)          normx=self.norm(x0)
# Line 645  class Defect(object): Line 667  class Defect(object):
667          return (F1-F0)/epsnew          return (F1-F0)/epsnew
668    
669  ######################################  ######################################
670  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):
671     """     """
672     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
673     criterion:     criterion:
674    
675     M{norm(F(x) <= atol + rtol * norm(F(x0)}     *norm(F(x) <= atol + rtol * norm(F(x0)*
676    
677     where M{x0} is the initial guess.     where *x0* is the initial guess.
678    
679     @param defect: object defining the function M{F}. C{defect.norm} defines the     :param defect: object defining the function *F*. ``defect.norm`` defines the
680                    M{norm} used in the stopping criterion.                    *norm* used in the stopping criterion.
681     @type defect: L{Defect}     :type defect: `Defect`
682     @param x: initial guess for the solution, C{x} is altered.     :param x: initial guess for the solution, ``x`` is altered.
683     @type x: any object type allowing basic operations such as     :type x: any object type allowing basic operations such as
684              C{numpy.ndarray}, L{Data}              ``numpy.ndarray``, `Data`
685     @param iter_max: maximum number of iteration steps     :param iter_max: maximum number of iteration steps
686     @type iter_max: positive C{int}     :type iter_max: positive ``int``
687     @param sub_iter_max: maximum number of inner iteration steps     :param sub_iter_max: maximum number of inner iteration steps
688     @type sub_iter_max: positive C{int}     :type sub_iter_max: positive ``int``
689     @param atol: absolute tolerance for the solution     :param atol: absolute tolerance for the solution
690     @type atol: positive C{float}     :type atol: positive ``float``
691     @param rtol: relative tolerance for the solution     :param rtol: relative tolerance for the solution
692     @type rtol: positive C{float}     :type rtol: positive ``float``
693     @param gamma: tolerance safety factor for inner iteration     :param gamma: tolerance safety factor for inner iteration
694     @type gamma: positive C{float}, less than 1     :type gamma: positive ``float``, less than 1
695     @param sub_tol_max: upper bound for inner tolerance     :param subtol_max: upper bound for inner tolerance
696     @type sub_tol_max: positive C{float}, less than 1     :type subtol_max: positive ``float``, less than 1
697     @return: an approximation of the solution with the desired accuracy     :return: an approximation of the solution with the desired accuracy
698     @rtype: same type as the initial guess     :rtype: same type as the initial guess
699     """     """
700     lmaxit=iter_max     lmaxit=iter_max
701     if atol<0: raise ValueError,"atol needs to be non-negative."     if atol<0: raise ValueError,"atol needs to be non-negative."
702     if rtol<0: raise ValueError,"rtol needs to be non-negative."     if rtol<0: raise ValueError,"rtol needs to be non-negative."
703     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."
704     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
705     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
706    
707     F=defect(x)     F=defect(x)
708     fnrm=defect.norm(F)     fnrm=defect.norm(F)
709     stop_tol=atol + rtol*fnrm     stop_tol=atol + rtol*fnrm
710     sub_tol=sub_tol_max     subtol=subtol_max
711     if verbose: print "NewtonGMRES: initial residual = %e."%fnrm     if verbose: print "NewtonGMRES: initial residual = %e."%fnrm
712     if verbose: print "             tolerance = %e."%sub_tol     if verbose: print "             tolerance = %e."%subtol
713     iter=1     iter=1
714     #     #
715     # main iteration loop     # main iteration loop
# Line 695  def NewtonGMRES(defect, x, iter_max=100, Line 717  def NewtonGMRES(defect, x, iter_max=100,
717     while not fnrm<=stop_tol:     while not fnrm<=stop_tol:
718              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
719              #              #
720          #   adjust sub_tol_          #   adjust subtol_
721          #          #
722              if iter > 1:              if iter > 1:
723             rat=fnrm/fnrmo             rat=fnrm/fnrmo
724                 sub_tol_old=sub_tol                 subtol_old=subtol
725             sub_tol=gamma*rat**2             subtol=gamma*rat**2
726             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)
727             sub_tol=max(min(sub_tol,sub_tol_max), .5*stop_tol/fnrm)             subtol=max(min(subtol,subtol_max), .5*stop_tol/fnrm)
728          #          #
729          # calculate newton increment xc          # calculate newton increment xc
730              #     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 732  def NewtonGMRES(defect, x, iter_max=100,
732              #     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
733              #              #
734              #              #
735              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
736              try:              try:
737                 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)
738              except MaxIterReached:              except MaxIterReached:
739                 raise MaxIterReached,"maximum number of %s steps reached."%iter_max                 raise MaxIterReached,"maximum number of %s steps reached."%iter_max
740              if sub_iter<0:              if sub_iter<0:
# Line 731  def NewtonGMRES(defect, x, iter_max=100, Line 753  def NewtonGMRES(defect, x, iter_max=100,
753  def __givapp(c,s,vin):  def __givapp(c,s,vin):
754      """      """
755      Applies a sequence of Givens rotations (c,s) recursively to the vector      Applies a sequence of Givens rotations (c,s) recursively to the vector
756      C{vin}      ``vin``
757    
758      @warning: C{vin} is altered.      :warning: ``vin`` is altered.
759      """      """
760      vrot=vin      vrot=vin
761      if isinstance(c,float):      if isinstance(c,float):
# Line 838  def GMRES(r, Aprod, x, bilinearform, ato Line 860  def GMRES(r, Aprod, x, bilinearform, ato
860     """     """
861     Solver for     Solver for
862    
863     M{Ax=b}     *Ax=b*
864    
865     with a general operator A (more details required!).     with a general operator A (more details required!).
866     It uses the generalized minimum residual method (GMRES).     It uses the generalized minimum residual method (GMRES).
867    
868     The iteration is terminated if     The iteration is terminated if
869    
870     M{|r| <= atol+rtol*|r0|}     *|r| <= atol+rtol*|r0|*
871    
872     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
873    
874     M{|r| = sqrt( bilinearform(r,r))}     *|r| = sqrt( bilinearform(r,r))*
875    
876     @param r: initial residual M{r=b-Ax}. C{r} is altered.     :param r: initial residual *r=b-Ax*. ``r`` is altered.
877     @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)
878     @param x: an initial guess for the solution     :param x: an initial guess for the solution
879     @type x: same like C{r}     :type x: same like ``r``
880     @param Aprod: returns the value Ax     :param Aprod: returns the value Ax
881     @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
882                  argument C{x}. The returned object needs to be of the same                  argument ``x``. The returned object needs to be of the same
883                  type like argument C{r}.                  type like argument ``r``.
884     @param bilinearform: inner product C{<x,r>}     :param bilinearform: inner product ``<x,r>``
885     @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
886                         type like argument C{x} and C{r}. The returned value is                         type like argument ``x`` and ``r``. The returned value is
887                         a C{float}.                         a ``float``.
888     @param atol: absolute tolerance     :param atol: absolute tolerance
889     @type atol: non-negative C{float}     :type atol: non-negative ``float``
890     @param rtol: relative tolerance     :param rtol: relative tolerance
891     @type rtol: non-negative C{float}     :type rtol: non-negative ``float``
892     @param iter_max: maximum number of iteration steps     :param iter_max: maximum number of iteration steps
893     @type iter_max: C{int}     :type iter_max: ``int``
894     @param iter_restart: in order to save memory the orthogonalization process     :param iter_restart: in order to save memory the orthogonalization process
895                          is terminated after C{iter_restart} steps and the                          is terminated after ``iter_restart`` steps and the
896                          iteration is restarted.                          iteration is restarted.
897     @type iter_restart: C{int}     :type iter_restart: ``int``
898     @return: the solution approximation and the corresponding residual     :return: the solution approximation and the corresponding residual
899     @rtype: C{tuple}     :rtype: ``tuple``
900     @warning: C{r} and C{x} are altered.     :warning: ``r`` and ``x`` are altered.
901     """     """
902     m=iter_restart     m=iter_restart
903     restarted=False     restarted=False
# Line 1007  def MINRES(r, Aprod, x, Msolve, bilinear Line 1029  def MINRES(r, Aprod, x, Msolve, bilinear
1029      """      """
1030      Solver for      Solver for
1031    
1032      M{Ax=b}      *Ax=b*
1033    
1034      with a symmetric and positive definite operator A (more details required!).      with a symmetric and positive definite operator A (more details required!).
1035      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 1037  def MINRES(r, Aprod, x, Msolve, bilinear
1037    
1038      The iteration is terminated if      The iteration is terminated if
1039    
1040      M{|r| <= atol+rtol*|r0|}      *|r| <= atol+rtol*|r0|*
1041    
1042      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
1043    
1044      M{|r| = sqrt( bilinearform(Msolve(r),r))}      *|r| = sqrt( bilinearform(Msolve(r),r))*
1045    
1046      For details on the preconditioned conjugate gradient method see the book:      For details on the preconditioned conjugate gradient method see the book:
1047    
# Line 1027  def MINRES(r, Aprod, x, Msolve, bilinear Line 1049  def MINRES(r, Aprod, x, Msolve, bilinear
1049      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,
1050      C. Romine, and H. van der Vorst}.      C. Romine, and H. van der Vorst}.
1051    
1052      @param r: initial residual M{r=b-Ax}. C{r} is altered.      :param r: initial residual *r=b-Ax*. ``r`` is altered.
1053      @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)
1054      @param x: an initial guess for the solution      :param x: an initial guess for the solution
1055      @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)
1056      @param Aprod: returns the value Ax      :param Aprod: returns the value Ax
1057      @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
1058                   argument C{x}. The returned object needs to be of the same                   argument ``x``. The returned object needs to be of the same
1059                   type like argument C{r}.                   type like argument ``r``.
1060      @param Msolve: solves Mx=r      :param Msolve: solves Mx=r
1061      @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
1062                    argument C{r}. The returned object needs to be of the same                    argument ``r``. The returned object needs to be of the same
1063                    type like argument C{x}.                    type like argument ``x``.
1064      @param bilinearform: inner product C{<x,r>}      :param bilinearform: inner product ``<x,r>``
1065      @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
1066                          type like argument C{x} and C{r} is. The returned value                          type like argument ``x`` and ``r`` is. The returned value
1067                          is a C{float}.                          is a ``float``.
1068      @param atol: absolute tolerance      :param atol: absolute tolerance
1069      @type atol: non-negative C{float}      :type atol: non-negative ``float``
1070      @param rtol: relative tolerance      :param rtol: relative tolerance
1071      @type rtol: non-negative C{float}      :type rtol: non-negative ``float``
1072      @param iter_max: maximum number of iteration steps      :param iter_max: maximum number of iteration steps
1073      @type iter_max: C{int}      :type iter_max: ``int``
1074      @return: the solution approximation and the corresponding residual      :return: the solution approximation and the corresponding residual
1075      @rtype: C{tuple}      :rtype: ``tuple``
1076      @warning: C{r} and C{x} are altered.      :warning: ``r`` and ``x`` are altered.
1077      """      """
1078      #------------------------------------------------------------------      #------------------------------------------------------------------
1079      # 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 1207  def TFQMR(r, Aprod, x, bilinearform, ato
1207    """    """
1208    Solver for    Solver for
1209    
1210    M{Ax=b}    *Ax=b*
1211    
1212    with a general operator A (more details required!).    with a general operator A (more details required!).
1213    It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).    It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).
1214    
1215    The iteration is terminated if    The iteration is terminated if
1216    
1217    M{|r| <= atol+rtol*|r0|}    *|r| <= atol+rtol*|r0|*
1218    
1219    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
1220    
1221    M{|r| = sqrt( bilinearform(r,r))}    *|r| = sqrt( bilinearform(r,r))*
1222    
1223    @param r: initial residual M{r=b-Ax}. C{r} is altered.    :param r: initial residual *r=b-Ax*. ``r`` is altered.
1224    @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)
1225    @param x: an initial guess for the solution    :param x: an initial guess for the solution
1226    @type x: same like C{r}    :type x: same like ``r``
1227    @param Aprod: returns the value Ax    :param Aprod: returns the value Ax
1228    @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
1229                 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
1230                 like argument C{r}.                 like argument ``r``.
1231    @param bilinearform: inner product C{<x,r>}    :param bilinearform: inner product ``<x,r>``
1232    @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
1233                        type like argument C{x} and C{r}. The returned value is                        type like argument ``x`` and ``r``. The returned value is
1234                        a C{float}.                        a ``float``.
1235    @param atol: absolute tolerance    :param atol: absolute tolerance
1236    @type atol: non-negative C{float}    :type atol: non-negative ``float``
1237    @param rtol: relative tolerance    :param rtol: relative tolerance
1238    @type rtol: non-negative C{float}    :type rtol: non-negative ``float``
1239    @param iter_max: maximum number of iteration steps    :param iter_max: maximum number of iteration steps
1240    @type iter_max: C{int}    :type iter_max: ``int``
1241    @rtype: C{tuple}    :rtype: ``tuple``
1242    @warning: C{r} and C{x} are altered.    :warning: ``r`` and ``x`` are altered.
1243    """    """
1244    u1=0    u1=0
1245    u2=0    u2=0
# Line 1287  def TFQMR(r, Aprod, x, bilinearform, ato Line 1309  def TFQMR(r, Aprod, x, bilinearform, ato
1309    
1310  class ArithmeticTuple(object):  class ArithmeticTuple(object):
1311     """     """
1312     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
1313     ArithmeticTuple and C{a} is a float.     ArithmeticTuple and ``a`` is a float.
1314    
1315     Example of usage::     Example of usage::
1316    
# Line 1302  class ArithmeticTuple(object): Line 1324  class ArithmeticTuple(object):
1324     """     """
1325     def __init__(self,*args):     def __init__(self,*args):
1326         """         """
1327         Initializes object with elements C{args}.         Initializes object with elements ``args``.
1328    
1329         @param args: tuple of objects that support inplace add (x+=y) and         :param args: tuple of objects that support inplace add (x+=y) and
1330                      scaling (x=a*y)                      scaling (x=a*y)
1331         """         """
1332         self.__items=list(args)         self.__items=list(args)
# Line 1313  class ArithmeticTuple(object): Line 1335  class ArithmeticTuple(object):
1335         """         """
1336         Returns the number of items.         Returns the number of items.
1337    
1338         @return: number of items         :return: number of items
1339         @rtype: C{int}         :rtype: ``int``
1340         """         """
1341         return len(self.__items)         return len(self.__items)
1342    
# Line 1322  class ArithmeticTuple(object): Line 1344  class ArithmeticTuple(object):
1344         """         """
1345         Returns item at specified position.         Returns item at specified position.
1346    
1347         @param index: index of item to be returned         :param index: index of item to be returned
1348         @type index: C{int}         :type index: ``int``
1349         @return: item with index C{index}         :return: item with index ``index``
1350         """         """
1351         return self.__items.__getitem__(index)         return self.__items.__getitem__(index)
1352    
1353     def __mul__(self,other):     def __mul__(self,other):
1354         """         """
1355         Scales by C{other} from the right.         Scales by ``other`` from the right.
1356    
1357         @param other: scaling factor         :param other: scaling factor
1358         @type other: C{float}         :type other: ``float``
1359         @return: itemwise self*other         :return: itemwise self*other
1360         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1361         """         """
1362         out=[]         out=[]
1363         try:         try:
# Line 1349  class ArithmeticTuple(object): Line 1371  class ArithmeticTuple(object):
1371    
1372     def __rmul__(self,other):     def __rmul__(self,other):
1373         """         """
1374         Scales by C{other} from the left.         Scales by ``other`` from the left.
1375    
1376         @param other: scaling factor         :param other: scaling factor
1377         @type other: C{float}         :type other: ``float``
1378         @return: itemwise other*self         :return: itemwise other*self
1379         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1380         """         """
1381         out=[]         out=[]
1382         try:         try:
# Line 1368  class ArithmeticTuple(object): Line 1390  class ArithmeticTuple(object):
1390    
1391     def __div__(self,other):     def __div__(self,other):
1392         """         """
1393         Scales by (1/C{other}) from the right.         Scales by (1/``other``) from the right.
1394    
1395         @param other: scaling factor         :param other: scaling factor
1396         @type other: C{float}         :type other: ``float``
1397         @return: itemwise self/other         :return: itemwise self/other
1398         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1399         """         """
1400         return self*(1/other)         return self*(1/other)
1401    
1402     def __rdiv__(self,other):     def __rdiv__(self,other):
1403         """         """
1404         Scales by (1/C{other}) from the left.         Scales by (1/``other``) from the left.
1405    
1406         @param other: scaling factor         :param other: scaling factor
1407         @type other: C{float}         :type other: ``float``
1408         @return: itemwise other/self         :return: itemwise other/self
1409         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1410         """         """
1411         out=[]         out=[]
1412         try:         try:
# Line 1398  class ArithmeticTuple(object): Line 1420  class ArithmeticTuple(object):
1420    
1421     def __iadd__(self,other):     def __iadd__(self,other):
1422         """         """
1423         Inplace addition of C{other} to self.         Inplace addition of ``other`` to self.
1424    
1425         @param other: increment         :param other: increment
1426         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1427         """         """
1428         if len(self) != len(other):         if len(self) != len(other):
1429             raise ValueError,"tuple lengths must match."             raise ValueError,"tuple lengths must match."
# Line 1411  class ArithmeticTuple(object): Line 1433  class ArithmeticTuple(object):
1433    
1434     def __add__(self,other):     def __add__(self,other):
1435         """         """
1436         Adds C{other} to self.         Adds ``other`` to self.
1437    
1438         @param other: increment         :param other: increment
1439         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1440         """         """
1441         out=[]         out=[]
1442         try:         try:
# Line 1428  class ArithmeticTuple(object): Line 1450  class ArithmeticTuple(object):
1450    
1451     def __sub__(self,other):     def __sub__(self,other):
1452         """         """
1453         Subtracts C{other} from self.         Subtracts ``other`` from self.
1454    
1455         @param other: decrement         :param other: decrement
1456         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1457         """         """
1458         out=[]         out=[]
1459         try:         try:
# Line 1445  class ArithmeticTuple(object): Line 1467  class ArithmeticTuple(object):
1467    
1468     def __isub__(self,other):     def __isub__(self,other):
1469         """         """
1470         Inplace subtraction of C{other} from self.         Inplace subtraction of ``other`` from self.
1471    
1472         @param other: decrement         :param other: decrement
1473         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1474         """         """
1475         if len(self) != len(other):         if len(self) != len(other):
1476             raise ValueError,"tuple length must match."             raise ValueError,"tuple length must match."
# Line 1471  class HomogeneousSaddlePointProblem(obje Line 1493  class HomogeneousSaddlePointProblem(obje
1493        This class provides a framework for solving linear homogeneous saddle        This class provides a framework for solving linear homogeneous saddle
1494        point problems of the form::        point problems of the form::
1495    
1496            M{Av+B^*p=f}            *Av+B^*p=f*
1497            M{Bv     =0}            *Bv     =0*
1498    
1499        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
1500        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*.
1501          *A* may depend weakly on *v* and *p*.
1502        """        """
1503        def __init__(self, adaptSubTolerance=True, **kwargs):        def __init__(self, **kwargs):
1504      """      """
1505      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}  
1506      """      """
1507            self.resetControlParameters()
1508          self.setTolerance()          self.setTolerance()
1509          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
1510      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):
1511             """
1512             sets a control parameter
1513    
1514             :param gamma: ``1/(1-gamma)`` controls the perturbation of the converegence rate due to termination errors in the subproblems.
1515             :type gamma: ``float``
1516             :param gamma_min: minimum value for ``gamma``.
1517             :type gamma_min: ``float``
1518             :param chi_max: maximum tolerable converegence rate.
1519             :type chi_max: ``float``
1520             :param omega_div: reduction fact for ``gamma`` if no convergence is detected.
1521             :type omega_div: ``float``
1522             :param omega_conv: raise fact for ``gamma`` if convergence is detected.
1523             :type omega_conv: ``float``
1524             :param rtol_min: minimum relative tolerance used to calculate presssure and velocity increment.
1525             :type rtol_min: ``float``
1526             :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1527             :type rtol_max: ``float``
1528             :param chi: initial convergence measure.
1529             :type chi: ``float``
1530             :param C_p: initial value for constant to adjust pressure tolerance
1531             :type C_p: ``float``
1532             :param C_v: initial value for constant to adjust velocity tolerance
1533             :type C_v: ``float``
1534             :param safety_factor: safety factor for addjustment of pressure and velocity tolerance from stopping criteria
1535             :type safety_factor: ``float``
1536             """
1537             self.setControlParameter(gamma, gamma_min ,chi_max , omega_div , omega_conv, rtol_min , rtol_max, chi,C_p, C_v,safety_factor)
1538    
1539          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):
1540             """
1541             sets a control parameter
1542    
1543             :param gamma: ``1/(1-gamma)`` controls the perturbation of the converegence rate due to termination errors in the subproblems.
1544             :type gamma: ``float``
1545             :param gamma_min: minimum value for ``gamma``.
1546             :type gamma_min: ``float``
1547             :param chi_max: maximum tolerable converegence rate.
1548             :type chi_max: ``float``
1549             :param omega_div: reduction fact for ``gamma`` if no convergence is detected.
1550             :type omega_div: ``float``
1551             :param omega_conv: raise fact for ``gamma`` if convergence is detected.
1552             :type omega_conv: ``float``
1553             :param rtol_min: minimum relative tolerance used to calculate presssure and velocity increment.
1554             :type rtol_min: ``float``
1555             :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1556             :type rtol_max: ``float``
1557             :param chi: initial convergence measure.
1558             :type chi: ``float``
1559             :param C_p: initial value for constant to adjust pressure tolerance
1560             :type C_p: ``float``
1561             :param C_v: initial value for constant to adjust velocity tolerance
1562             :type C_v: ``float``
1563             :param safety_factor: safety factor for addjustment of pressure and velocity tolerance from stopping criteria
1564             :type safety_factor: ``float``
1565             """
1566             if not gamma == None:
1567                if gamma<=0 or gamma>=1:
1568                   raise ValueError,"Initial gamma needs to be positive and less than 1."
1569             else:
1570                gamma=self.__gamma
1571    
1572             if not gamma_min == None:
1573                if gamma_min<=0 or gamma_min>=1:
1574                   raise ValueError,"gamma_min needs to be positive and less than 1."
1575             else:
1576                gamma_min = self.__gamma_min
1577    
1578             if not chi_max == None:
1579                if chi_max<=0 or chi_max>=1:
1580                   raise ValueError,"chi_max needs to be positive and less than 1."
1581             else:
1582                chi_max = self.__chi_max
1583    
1584             if not omega_div == None:
1585                if omega_div<=0 or omega_div >=1:
1586                   raise ValueError,"omega_div needs to be positive and less than 1."
1587             else:
1588                omega_div=self.__omega_div
1589    
1590             if not omega_conv == None:
1591                if omega_conv<1:
1592                   raise ValueError,"omega_conv needs to be greater or equal to 1."
1593             else:
1594                omega_conv=self.__omega_conv
1595    
1596             if not rtol_min == None:
1597                if rtol_min<=0 or rtol_min>=1:
1598                   raise ValueError,"rtol_min needs to be positive and less than 1."
1599             else:
1600                rtol_min=self.__rtol_min
1601    
1602             if not rtol_max == None:
1603                if rtol_max<=0 or rtol_max>=1:
1604                   raise ValueError,"rtol_max needs to be positive and less than 1."
1605             else:
1606                rtol_max=self.__rtol_max
1607    
1608             if not chi == None:
1609                if chi<=0 or chi>1:
1610                   raise ValueError,"chi needs to be positive and less than 1."
1611             else:
1612                chi=self.__chi
1613    
1614             if not C_p == None:
1615                if C_p<1:
1616                   raise ValueError,"C_p need to be greater or equal to 1."
1617             else:
1618                C_p=self.__C_p
1619    
1620             if not C_v == None:
1621                if C_v<1:
1622                   raise ValueError,"C_v need to be greater or equal to 1."
1623             else:
1624                C_v=self.__C_v
1625    
1626             if not safety_factor == None:
1627                if safety_factor<=0 or safety_factor>1:
1628                   raise ValueError,"safety_factor need to be between zero and one."
1629             else:
1630                safety_factor=self.__safety_factor
1631    
1632             if gamma<gamma_min:
1633                   raise ValueError,"gamma = %e needs to be greater or equal gamma_min = %e."%(gamma,gamma_min)
1634             if rtol_max<=rtol_min:
1635                   raise ValueError,"rtol_max = %e needs to be greater rtol_min = %e."%(rtol_max,rtol_min)
1636                
1637             self.__gamma = gamma
1638             self.__gamma_min = gamma_min
1639             self.__chi_max = chi_max
1640             self.__omega_div = omega_div
1641             self.__omega_conv = omega_conv
1642             self.__rtol_min = rtol_min
1643             self.__rtol_max = rtol_max
1644             self.__chi = chi
1645             self.__C_p = C_p
1646             self.__C_v = C_v
1647             self.__safety_factor = safety_factor
1648    
1649        #=============================================================        #=============================================================
1650        def initialize(self):        def initialize(self):
1651          """          """
# Line 1498  class HomogeneousSaddlePointProblem(obje Line 1657  class HomogeneousSaddlePointProblem(obje
1657           """           """
1658           Returns inner product of element p and Bv (overwrite).           Returns inner product of element p and Bv (overwrite).
1659    
1660           @param p: a pressure increment           :param p: a pressure increment
1661           @param v: a residual           :param Bv: a residual
1662           @return: inner product of element p and Bv           :return: inner product of element p and Bv
1663           @rtype: C{float}           :rtype: ``float``
1664           @note: used if PCG is applied.           :note: used if PCG is applied.
1665           """           """
1666           raise NotImplementedError,"no inner product for p and Bv implemented."           raise NotImplementedError,"no inner product for p and Bv implemented."
1667    
# Line 1510  class HomogeneousSaddlePointProblem(obje Line 1669  class HomogeneousSaddlePointProblem(obje
1669           """           """
1670           Returns inner product of p0 and p1 (overwrite).           Returns inner product of p0 and p1 (overwrite).
1671    
1672           @param p0: a pressure           :param p0: a pressure
1673           @param p1: a pressure           :param p1: a pressure
1674           @return: inner product of p0 and p1           :return: inner product of p0 and p1
1675           @rtype: C{float}           :rtype: ``float``
1676           """           """
1677           raise NotImplementedError,"no inner product for p implemented."           raise NotImplementedError,"no inner product for p implemented."
1678        
# Line 1521  class HomogeneousSaddlePointProblem(obje Line 1680  class HomogeneousSaddlePointProblem(obje
1680           """           """
1681           Returns the norm of v (overwrite).           Returns the norm of v (overwrite).
1682    
1683           @param v: a velovity           :param v: a velovity
1684           @return: norm of v           :return: norm of v
1685           @rtype: non-negative C{float}           :rtype: non-negative ``float``
1686           """           """
1687           raise NotImplementedError,"no norm of v implemented."           raise NotImplementedError,"no norm of v implemented."
1688        def getV(self, p, v0):        def getDV(self, p, v, tol):
1689           """           """
1690           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)
1691    
1692           @param p: a pressure           :param p: pressure
1693           @param v0: a initial guess for the value v to return.           :param v: pressure
1694           @return: v given as M{v= A^{-1} (f-B^*p)}           :return: dv given as *dv= A^{-1} (f-A v-B^*p)*
1695             :note: Only *A* may depend on *v* and *p*
1696           """           """
1697           raise NotImplementedError,"no v calculation implemented."           raise NotImplementedError,"no dv calculation implemented."
1698    
1699                    
1700        def Bv(self,v):        def Bv(self,v, tol):
1701          """          """
1702          Returns Bv (overwrite).          Returns Bv with accuracy `tol` (overwrite)
1703    
1704          @rtype: equal to the type of p          :rtype: equal to the type of p
1705          @note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1706          """          """
1707          raise NotImplementedError, "no operator B implemented."          raise NotImplementedError, "no operator B implemented."
1708    
# Line 1550  class HomogeneousSaddlePointProblem(obje Line 1710  class HomogeneousSaddlePointProblem(obje
1710          """          """
1711          Returns the norm of Bv (overwrite).          Returns the norm of Bv (overwrite).
1712    
1713          @rtype: equal to the type of p          :rtype: equal to the type of p
1714          @note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1715          """          """
1716          raise NotImplementedError, "no norm of Bv implemented."          raise NotImplementedError, "no norm of Bv implemented."
1717    
1718        def solve_AinvBt(self,p):        def solve_AinvBt(self,dp, tol):
1719           """           """
1720           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}           Solves *A dv=B^*dp* with accuracy `tol`
          (overwrite).  
1721    
1722           @param p: a pressure increment           :param dp: a pressure increment
1723           @return: the solution of M{Av=B^*p}           :return: the solution of *A dv=B^*dp*
1724           @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.
1725           """           """
1726           raise NotImplementedError,"no operator A implemented."           raise NotImplementedError,"no operator A implemented."
1727    
1728        def solve_prec(self,Bv):        def solve_prec(self,Bv, tol):
1729           """           """
1730           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).  
1731    
1732           @rtype: equal to the type of p           :rtype: equal to the type of p
1733           @note: boundary conditions on p should be zero!           :note: boundary conditions on p should be zero!
1734           """           """
1735           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  
1736        #=============================================================        #=============================================================
1737        def __Aprod_PCG(self,p):        def __Aprod_PCG(self,dp):
1738            dv=self.solve_AinvBt(p)            dv=self.solve_AinvBt(dp, self.__subtol)
1739            return ArithmeticTuple(dv,self.Bv(dv))            return ArithmeticTuple(dv,self.Bv(dv, self.__subtol))
1740    
1741        def __inner_PCG(self,p,r):        def __inner_PCG(self,p,r):
1742           return self.inner_pBv(p,r[1])           return self.inner_pBv(p,r[1])
1743    
1744        def __Msolve_PCG(self,r):        def __Msolve_PCG(self,r):
1745            return self.solve_prec(r[1])            return self.solve_prec(r[1], self.__subtol)
1746        #=============================================================        #=============================================================
1747        def __Aprod_GMRES(self,p):        def __Aprod_GMRES(self,p):
1748            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)
1749        def __inner_GMRES(self,p0,p1):        def __inner_GMRES(self,p0,p1):
1750           return self.inner_p(p0,p1)           return self.inner_p(p0,p1)
1751    
1752        #=============================================================        #=============================================================
1753        def norm_p(self,p):        def norm_p(self,p):
1754            """            """
1755            calculates the norm of C{p}            calculates the norm of ``p``
1756                        
1757            @param p: a pressure            :param p: a pressure
1758            @return: the norm of C{p} using the inner product for pressure            :return: the norm of ``p`` using the inner product for pressure
1759            @rtype: C{float}            :rtype: ``float``
1760            """            """
1761            f=self.inner_p(p,p)            f=self.inner_p(p,p)
1762            if f<0: raise ValueError,"negative pressure norm."            if f<0: raise ValueError,"negative pressure norm."
1763            return math.sqrt(f)            return math.sqrt(f)
       def adaptSubTolerance(self):  
       """  
       Returns True if tolerance adaption for subproblem is choosen.  
       """  
           self.__adaptSubTolerance  
1764                
1765        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):
1766           """           """
1767           Solves the saddle point problem using initial guesses v and p.           Solves the saddle point problem using initial guesses v and p.
1768    
1769           @param v: initial guess for velocity           :param v: initial guess for velocity
1770           @param p: initial guess for pressure           :param p: initial guess for pressure
1771           @type v: L{Data}           :type v: `Data`
1772           @type p: L{Data}           :type p: `Data`
1773           @param usePCG: indicates the usage of the PCG rather than GMRES scheme.           :param usePCG: indicates the usage of the PCG rather than GMRES scheme.
1774           @param max_iter: maximum number of iteration steps per correction           :param max_iter: maximum number of iteration steps per correction
1775                            attempt                            attempt
1776           @param verbose: if True, shows information on the progress of the           :param verbose: if True, shows information on the progress of the
1777                           saddlepoint problem solver.                           saddlepoint problem solver.
1778           @param iter_restart: restart the iteration after C{iter_restart} steps           :param iter_restart: restart the iteration after ``iter_restart`` steps
1779                                (only used if useUzaw=False)                                (only used if useUzaw=False)
1780           @type usePCG: C{bool}           :type usePCG: ``bool``
1781           @type max_iter: C{int}           :type max_iter: ``int``
1782           @type verbose: C{bool}           :type verbose: ``bool``
1783           @type iter_restart: C{int}           :type iter_restart: ``int``
1784           @rtype: C{tuple} of L{Data} objects           :rtype: ``tuple`` of `Data` objects
1785           """           """
1786           self.verbose=verbose           self.verbose=verbose
1787           rtol=self.getTolerance()           rtol=self.getTolerance()
1788           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
      if self.adaptSubTolerance(): self.setSubProblemTolerance()  
1789           correction_step=0           correction_step=0
1790           converged=False           converged=False
1791             error=None
1792             chi=None
1793             gamma=self.__gamma
1794             C_p=self.__C_p
1795             C_v=self.__C_v
1796           while not converged:           while not converged:
1797                  if error== None or chi == None:
1798                      gamma_new=gamma/self.__omega_conv
1799                  else:
1800                     if chi < self.__chi_max:
1801                        gamma_new=min(max(gamma*self.__omega_conv,1-chi*error/(self.__safety_factor*ATOL)), 1-chi/self.__chi_max)
1802                     else:
1803                        gamma_new=gamma*self.__omega_div
1804                  gamma=max(gamma_new, self.__gamma_min)
1805                # calculate velocity for current pressure:                # calculate velocity for current pressure:
1806                v=self.getV(p,v)                rtol_v=min(max(gamma/(1.+gamma)/C_v,self.__rtol_min), self.__rtol_max)
1807                Bv=self.Bv(v)                rtol_p=min(max(gamma/C_p, self.__rtol_min), self.__rtol_max)
1808                norm_v=self.norm_v(v)                self.__subtol=rtol_p**2
1809                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)
1810                ATOL=norm_v*rtol+atol                if self.verbose: print "HomogeneousSaddlePointProblem: subtolerance: %e"%self.__subtol
1811                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
1812                  dv1=self.getDV(p,v,rtol_v)
1813                  v1=v+dv1
1814                  Bv1=self.Bv(v1, self.__subtol)
1815                  norm_Bv1=self.norm_Bv(Bv1)
1816                  norm_dv1=self.norm_v(dv1)
1817                  norm_v1=self.norm_v(v1)
1818                  ATOL=norm_v1*rtol+atol
1819                  if self.verbose: print "HomogeneousSaddlePointProblem: step %s: Bv = %e, dv = %e, v=%e"%(correction_step,norm_Bv1, norm_dv1, norm_v1)
1820                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."
1821                if norm_Bv <= ATOL:                if max(norm_Bv1, norm_dv1) <= ATOL:
1822                   converged=True                    converged=True
1823                      v=v1
1824                else:                else:
1825                   correction_step+=1                    # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1
1826                   if correction_step>max_correction_steps:                    if usePCG:
1827                        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)
1828                   dp=self.solve_prec(Bv)                      v2=r[0]
1829                   if usePCG:                      Bv2=r[1]
1830                     norm2=self.inner_pBv(dp,Bv)                    else:
1831                     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)
1832                     norm2=math.sqrt(norm2)                      dv2=self.solve_AinvBt(dp, self.__subtol)
1833                   else:                      v2=v1-dv2
1834                     norm2=self.norm_p(dp)                      Bv2=self.Bv(v2, self.__subtol)
1835                   ATOL_ITER=ATOL/norm_Bv*norm2*0.5                    #
1836                   if self.verbose: print "HomogeneousSaddlePointProblem: tolerance for solver: %e"%ATOL_ITER                    # convergence indicators:
1837                   if usePCG:                    #
1838                         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)
1839                   else:                    norm_dv2=self.norm_v(v2-v)
1840                         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)
1841           if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached."                    norm_Bv2=self.norm_Bv(Bv2)
1842                      if self.verbose: print "HomogeneousSaddlePointProblem: step %s (part 2): Bv = %e, dv = %e, v=%e"%(correction_step,norm_Bv2, norm_dv2, norm_v2)
1843                      if error !=None:
1844                          chi_new=error_new/error
1845                          if self.verbose: print "HomogeneousSaddlePointProblem: step %s: convergence rate = %e"%(correction_step,chi_new)
1846                          if chi != None:
1847                              gamma0=max(gamma, 1-chi/chi_new)
1848                              C_p*=gamma0/gamma
1849                              C_v*=gamma0/gamma*(1+gamma)/(1+gamma0)
1850                          chi=chi_new
1851                      error = error_new
1852                      correction_step+=1
1853                      if correction_step>max_correction_steps:
1854                            raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1855                      v,p=v2,p+dp
1856             if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step
1857       return v,p       return v,p
1858    
1859        #========================================================================        #========================================================================
# Line 1678  class HomogeneousSaddlePointProblem(obje Line 1861  class HomogeneousSaddlePointProblem(obje
1861           """           """
1862           Sets the relative tolerance for (v,p).           Sets the relative tolerance for (v,p).
1863    
1864           @param tolerance: tolerance to be used           :param tolerance: tolerance to be used
1865           @type tolerance: non-negative C{float}           :type tolerance: non-negative ``float``
1866           """           """
1867           if tolerance<0:           if tolerance<0:
1868               raise ValueError,"tolerance must be positive."               raise ValueError,"tolerance must be positive."
# Line 1689  class HomogeneousSaddlePointProblem(obje Line 1872  class HomogeneousSaddlePointProblem(obje
1872           """           """
1873           Returns the relative tolerance.           Returns the relative tolerance.
1874    
1875           @return: relative tolerance           :return: relative tolerance
1876           @rtype: C{float}           :rtype: ``float``
1877           """           """
1878           return self.__rtol           return self.__rtol
1879    
# Line 1698  class HomogeneousSaddlePointProblem(obje Line 1881  class HomogeneousSaddlePointProblem(obje
1881           """           """
1882           Sets the absolute tolerance.           Sets the absolute tolerance.
1883    
1884           @param tolerance: tolerance to be used           :param tolerance: tolerance to be used
1885           @type tolerance: non-negative C{float}           :type tolerance: non-negative ``float``
1886           """           """
1887           if tolerance<0:           if tolerance<0:
1888               raise ValueError,"tolerance must be non-negative."               raise ValueError,"tolerance must be non-negative."
# Line 1709  class HomogeneousSaddlePointProblem(obje Line 1892  class HomogeneousSaddlePointProblem(obje
1892           """           """
1893           Returns the absolute tolerance.           Returns the absolute tolerance.
1894    
1895           @return: absolute tolerance           :return: absolute tolerance
1896           @rtype: C{float}           :rtype: ``float``
1897           """           """
1898           return self.__atol           return self.__atol
1899    
       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)  
1900    
1901  def MaskFromBoundaryTag(domain,*tags):  def MaskFromBoundaryTag(domain,*tags):
1902     """     """
# Line 1730  def MaskFromBoundaryTag(domain,*tags): Line 1905  def MaskFromBoundaryTag(domain,*tags):
1905    
1906     Usage: m=MaskFromBoundaryTag(domain, "left", "right")     Usage: m=MaskFromBoundaryTag(domain, "left", "right")
1907    
1908     @param domain: domain to be used     :param domain: domain to be used
1909     @type domain: L{escript.Domain}     :type domain: `escript.Domain`
1910     @param tags: boundary tags     :param tags: boundary tags
1911     @type tags: C{str}     :type tags: ``str``
1912     @return: a mask which marks samples that are touching the boundary tagged     :return: a mask which marks samples that are touching the boundary tagged
1913              by any of the given tags              by any of the given tags
1914     @rtype: L{escript.Data} of rank 0     :rtype: `escript.Data` of rank 0
1915     """     """
1916     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1917     d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))     d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
# Line 1751  def MaskFromTag(domain,*tags): Line 1926  def MaskFromTag(domain,*tags):
1926    
1927     Usage: m=MaskFromTag(domain, "ham")     Usage: m=MaskFromTag(domain, "ham")
1928    
1929     @param domain: domain to be used     :param domain: domain to be used
1930     @type domain: L{escript.Domain}     :type domain: `escript.Domain`
1931     @param tags: boundary tags     :param tags: boundary tags
1932     @type tags: C{str}     :type tags: ``str``
1933     @return: a mask which marks samples that are touching the boundary tagged     :return: a mask which marks samples that are touching the boundary tagged
1934              by any of the given tags              by any of the given tags
1935     @rtype: L{escript.Data} of rank 0     :rtype: `escript.Data` of rank 0
1936     """     """
1937     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1938     d=escript.Scalar(0.,escript.Function(domain))     d=escript.Scalar(0.,escript.Function(domain))

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

  ViewVC Help
Powered by ViewVC 1.1.26