/[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 2881 by jfenwick, Thu Jan 28 02:03:15 2010 UTC
# Line 1  Line 1 
1    
2  ########################################################  ########################################################
3  #  #
4  # Copyright (c) 2003-2009 by University of Queensland  # Copyright (c) 2003-2010 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
6  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
7  #  #
# Line 11  Line 11 
11  #  #
12  ########################################################  ########################################################
13    
14  __copyright__="""Copyright (c) 2003-2008 by University of Queensland  __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
16  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
17  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
# Line 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 getValue(self, input_data):
157        """
158        Projects ``input_data`` onto a continuous function.
159    
160        :param input_data: the data to be projected
161        """
162        return self(input_data)
163    
164    def __call__(self, input_data):    def __call__(self, input_data):
165      """      """
166      Projects C{input_data} onto a continuous function.      Projects ``input_data`` onto a continuous function.
167    
168      @param input_data: the data to be projected      :param input_data: the data to be projected
169      """      """
170      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
171      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 200  class NoPDE:
200       """       """
201       Solves the following problem for u:       Solves the following problem for u:
202    
203       M{kronecker[i,j]*D[j]*u[j]=Y[i]}       *kronecker[i,j]*D[j]*u[j]=Y[i]*
204    
205       with constraint       with constraint
206    
207       M{u[j]=r[j]}  where M{q[j]>0}       *u[j]=r[j]*  where *q[j]>0*
208    
209       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.
210    
211       In the case of scalars this takes the form       In the case of scalars this takes the form
212    
213       M{D*u=Y}       *D*u=Y*
214    
215       with constraint       with constraint
216    
217       M{u=r} where M{q>0}       *u=r* where *q>0*
218    
219       where M{D}, M{Y}, M{r} and M{q} are given scalar functions.       where *D*, *Y*, *r* and *q* are given scalar functions.
220    
221       The constraint overwrites any other condition.       The constraint overwrites any other condition.
222    
223       @note: This class is similar to the L{linearPDEs.LinearPDE} class with       :note: This class is similar to the `linearPDEs.LinearPDE` class with
224              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
225              in L{Solution} or L{ReducedSolution}.              in `Solution` or `ReducedSolution`.
226       """       """
227       # 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
228       # this.       # this.
# Line 226  class NoPDE: Line 230  class NoPDE:
230           """           """
231           Initializes the problem.           Initializes the problem.
232    
233           @param domain: domain of the PDE           :param domain: domain of the PDE
234           @type domain: L{Domain}           :type domain: `Domain`
235           @param D: coefficient of the solution           :param D: coefficient of the solution
236           @type D: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type D: ``float``, ``int``, ``numpy.ndarray``, `Data`
237           @param Y: right hand side           :param Y: right hand side
238           @type Y: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type Y: ``float``, ``int``, ``numpy.ndarray``, `Data`
239           @param q: location of constraints           :param q: location of constraints
240           @type q: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type q: ``float``, ``int``, ``numpy.ndarray``, `Data`
241           @param r: value of solution at locations of constraints           :param r: value of solution at locations of constraints
242           @type r: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type r: ``float``, ``int``, ``numpy.ndarray``, `Data`
243           """           """
244           self.__domain=domain           self.__domain=domain
245           self.__D=D           self.__D=D
# Line 247  class NoPDE: Line 251  class NoPDE:
251    
252       def setReducedOn(self):       def setReducedOn(self):
253           """           """
254           Sets the L{FunctionSpace} of the solution to L{ReducedSolution}.           Sets the `FunctionSpace` of the solution to `ReducedSolution`.
255           """           """
256           self.__function_space=escript.ReducedSolution(self.__domain)           self.__function_space=escript.ReducedSolution(self.__domain)
257           self.__u=None           self.__u=None
258    
259       def setReducedOff(self):       def setReducedOff(self):
260           """           """
261           Sets the L{FunctionSpace} of the solution to L{Solution}.           Sets the `FunctionSpace` of the solution to `Solution`.
262           """           """
263           self.__function_space=escript.Solution(self.__domain)           self.__function_space=escript.Solution(self.__domain)
264           self.__u=None           self.__u=None
# Line 263  class NoPDE: Line 267  class NoPDE:
267           """           """
268           Assigns values to the parameters.           Assigns values to the parameters.
269    
270           @param D: coefficient of the solution           :param D: coefficient of the solution
271           @type D: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type D: ``float``, ``int``, ``numpy.ndarray``, `Data`
272           @param Y: right hand side           :param Y: right hand side
273           @type Y: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type Y: ``float``, ``int``, ``numpy.ndarray``, `Data`
274           @param q: location of constraints           :param q: location of constraints
275           @type q: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type q: ``float``, ``int``, ``numpy.ndarray``, `Data`
276           @param r: value of solution at locations of constraints           :param r: value of solution at locations of constraints
277           @type r: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type r: ``float``, ``int``, ``numpy.ndarray``, `Data`
278           """           """
279           if not D==None:           if not D==None:
280              self.__D=D              self.__D=D
# Line 289  class NoPDE: Line 293  class NoPDE:
293           """           """
294           Returns the solution.           Returns the solution.
295    
296           @return: the solution of the problem           :return: the solution of the problem
297           @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or           :rtype: `Data` object in the `FunctionSpace` `Solution` or
298                   L{ReducedSolution}                   `ReducedSolution`
299           """           """
300           if self.__u==None:           if self.__u==None:
301              if self.__D==None:              if self.__D==None:
# Line 324  class Locator: Line 328  class Locator:
328         or FunctionSpace for the sample point which is closest to the given         or FunctionSpace for the sample point which is closest to the given
329         point x.         point x.
330    
331         @param where: function space         :param where: function space
332         @type where: L{escript.FunctionSpace}         :type where: `escript.FunctionSpace`
333         @param x: location(s) of the Locator         :param x: location(s) of the Locator
334         @type x: C{numpy.ndarray} or C{list} of C{numpy.ndarray}         :type x: ``numpy.ndarray`` or ``list`` of ``numpy.ndarray``
335         """         """
336         if isinstance(where,escript.FunctionSpace):         if isinstance(where,escript.FunctionSpace):
337            self.__function_space=where            self.__function_space=where
# Line 354  class Locator: Line 358  class Locator:
358         Returns the coordinates of the Locator as a string.         Returns the coordinates of the Locator as a string.
359         """         """
360         x=self.getX()         x=self.getX()
361         if instance(x,list):         if isinstance(x,list):
362            out="["            out="["
363            first=True            first=True
364            for xx in x:            for xx in x:
# Line 401  class Locator: Line 405  class Locator:
405    
406       def getValue(self,data):       def getValue(self,data):
407          """          """
408          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`
409          object otherwise the object is returned.          object otherwise the object is returned.
410          """          """
411          if isinstance(data,escript.Data):          if isinstance(data,escript.Data):
# Line 426  class Locator: Line 430  class Locator:
430          else:          else:
431             return data             return data
432    
433    
434    def getInfLocator(arg):
435        """
436        Return a Locator for a point with the inf value over all arg.
437        """
438        if not isinstance(arg, escript.Data):
439        raise TypeError,"getInfLocator: Unknown argument type."
440        a_inf=util.inf(arg)
441        loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords
442        x=arg.getFunctionSpace().getX()
443        x_min=x.getTupleForGlobalDataPoint(*loc)
444        return Locator(arg.getFunctionSpace(),x_min)
445    
446    def getSupLocator(arg):
447        """
448        Return a Locator for a point with the sup value over all arg.
449        """
450        if not isinstance(arg, escript.Data):
451        raise TypeError,"getInfLocator: Unknown argument type."
452        a_inf=util.sup(arg)
453        loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords
454        x=arg.getFunctionSpace().getX()
455        x_min=x.getTupleForGlobalDataPoint(*loc)
456        return Locator(arg.getFunctionSpace(),x_min)
457        
458    
459  class SolverSchemeException(Exception):  class SolverSchemeException(Exception):
460     """     """
461     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 497  def PCG(r, Aprod, x, Msolve, bilinearfor
497     """     """
498     Solver for     Solver for
499    
500     M{Ax=b}     *Ax=b*
501    
502     with a symmetric and positive definite operator A (more details required!).     with a symmetric and positive definite operator A (more details required!).
503     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 505  def PCG(r, Aprod, x, Msolve, bilinearfor
505    
506     The iteration is terminated if     The iteration is terminated if
507    
508     M{|r| <= atol+rtol*|r0|}     *|r| <= atol+rtol*|r0|*
509    
510     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
511    
512     M{|r| = sqrt( bilinearform(Msolve(r),r))}     *|r| = sqrt( bilinearform(Msolve(r),r))*
513    
514     For details on the preconditioned conjugate gradient method see the book:     For details on the preconditioned conjugate gradient method see the book:
515    
# Line 487  def PCG(r, Aprod, x, Msolve, bilinearfor Line 517  def PCG(r, Aprod, x, Msolve, bilinearfor
517     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,
518     C. Romine, and H. van der Vorst}.     C. Romine, and H. van der Vorst}.
519    
520     @param r: initial residual M{r=b-Ax}. C{r} is altered.     :param r: initial residual *r=b-Ax*. ``r`` is altered.
521     @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)
522     @param x: an initial guess for the solution     :param x: an initial guess for the solution
523     @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)
524     @param Aprod: returns the value Ax     :param Aprod: returns the value Ax
525     @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
526                  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
527                  like argument C{r}.                  like argument ``r``.
528     @param Msolve: solves Mx=r     :param Msolve: solves Mx=r
529     @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
530                   argument C{r}. The returned object needs to be of the same                   argument ``r``. The returned object needs to be of the same
531                   type like argument C{x}.                   type like argument ``x``.
532     @param bilinearform: inner product C{<x,r>}     :param bilinearform: inner product ``<x,r>``
533     @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
534                         type like argument C{x} and C{r} is. The returned value                         type like argument ``x`` and ``r`` is. The returned value
535                         is a C{float}.                         is a ``float``.
536     @param atol: absolute tolerance     :param atol: absolute tolerance
537     @type atol: non-negative C{float}     :type atol: non-negative ``float``
538     @param rtol: relative tolerance     :param rtol: relative tolerance
539     @type rtol: non-negative C{float}     :type rtol: non-negative ``float``
540     @param iter_max: maximum number of iteration steps     :param iter_max: maximum number of iteration steps
541     @type iter_max: C{int}     :type iter_max: ``int``
542     @return: the solution approximation and the corresponding residual     :return: the solution approximation and the corresponding residual
543     @rtype: C{tuple}     :rtype: ``tuple``
544     @warning: C{r} and C{x} are altered.     :warning: ``r`` and ``x`` are altered.
545     """     """
546     iter=0     iter=0
547     rhat=Msolve(r)     rhat=Msolve(r)
# Line 564  class Defect(object): Line 594  class Defect(object):
594          """          """
595          Returns the inner product of x0 and x1          Returns the inner product of x0 and x1
596    
597          @param x0: value for x0          :param x0: value for x0
598          @param x1: value for x1          :param x1: value for x1
599          @return: the inner product of x0 and x1          :return: the inner product of x0 and x1
600          @rtype: C{float}          :rtype: ``float``
601          """          """
602          return 0          return 0
603    
604      def norm(self,x):      def norm(self,x):
605          """          """
606          Returns the norm of argument C{x}.          Returns the norm of argument ``x``.
607    
608          @param x: a value          :param x: a value
609          @return: norm of argument x          :return: norm of argument x
610          @rtype: C{float}          :rtype: ``float``
611          @note: by default C{sqrt(self.bilinearform(x,x)} is returned.          :note: by default ``sqrt(self.bilinearform(x,x)`` is returned.
612          """          """
613          s=self.bilinearform(x,x)          s=self.bilinearform(x,x)
614          if s<0: raise NegativeNorm,"negative norm."          if s<0: raise NegativeNorm,"negative norm."
# Line 586  class Defect(object): Line 616  class Defect(object):
616    
617      def eval(self,x):      def eval(self,x):
618          """          """
619          Returns the value F of a given C{x}.          Returns the value F of a given ``x``.
620    
621          @param x: value for which the defect C{F} is evaluated          :param x: value for which the defect ``F`` is evaluated
622          @return: value of the defect at C{x}          :return: value of the defect at ``x``
623          """          """
624          return 0          return 0
625    
626      def __call__(self,x):      def __call__(self,x):
627          return self.eval(x)          return self.eval(x)
628    
629      def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):      def setDerivativeIncrementLength(self,inc=1000.*math.sqrt(util.EPSILON)):
630          """          """
631          Sets the relative length of the increment used to approximate the          Sets the relative length of the increment used to approximate the
632          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
633          direction of v with x as a starting point.          direction of v with x as a starting point.
634    
635          @param inc: relative increment length          :param inc: relative increment length
636          @type inc: positive C{float}          :type inc: positive ``float``
637          """          """
638          if inc<=0: raise ValueError,"positive increment required."          if inc<=0: raise ValueError,"positive increment required."
639          self.__inc=inc          self.__inc=inc
# Line 612  class Defect(object): Line 642  class Defect(object):
642          """          """
643          Returns the relative increment length used to approximate the          Returns the relative increment length used to approximate the
644          derivative of the defect.          derivative of the defect.
645          @return: value of the defect at C{x}          :return: value of the defect at ``x``
646          @rtype: positive C{float}          :rtype: positive ``float``
647          """          """
648          return self.__inc          return self.__inc
649    
650      def derivative(self, F0, x0, v, v_is_normalised=True):      def derivative(self, F0, x0, v, v_is_normalised=True):
651          """          """
652          Returns the directional derivative at C{x0} in the direction of C{v}.          Returns the directional derivative at ``x0`` in the direction of ``v``.
653    
654          @param F0: value of this defect at x0          :param F0: value of this defect at x0
655          @param x0: value at which derivative is calculated          :param x0: value at which derivative is calculated
656          @param v: direction          :param v: direction
657          @param v_is_normalised: True to indicate that C{v} is nomalized          :param v_is_normalised: True to indicate that ``v`` is nomalized
658                                  (self.norm(v)=0)                                  (self.norm(v)=0)
659          @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``
660          @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
661                 used but this method maybe overwritten to use exact evaluation.                 used but this method maybe overwritten to use exact evaluation.
662          """          """
663          normx=self.norm(x0)          normx=self.norm(x0)
# Line 645  class Defect(object): Line 675  class Defect(object):
675          return (F1-F0)/epsnew          return (F1-F0)/epsnew
676    
677  ######################################  ######################################
678  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):
679     """     """
680     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
681     criterion:     criterion:
682    
683     M{norm(F(x) <= atol + rtol * norm(F(x0)}     *norm(F(x) <= atol + rtol * norm(F(x0)*
684    
685     where M{x0} is the initial guess.     where *x0* is the initial guess.
686    
687     @param defect: object defining the function M{F}. C{defect.norm} defines the     :param defect: object defining the function *F*. ``defect.norm`` defines the
688                    M{norm} used in the stopping criterion.                    *norm* used in the stopping criterion.
689     @type defect: L{Defect}     :type defect: `Defect`
690     @param x: initial guess for the solution, C{x} is altered.     :param x: initial guess for the solution, ``x`` is altered.
691     @type x: any object type allowing basic operations such as     :type x: any object type allowing basic operations such as
692              C{numpy.ndarray}, L{Data}              ``numpy.ndarray``, `Data`
693     @param iter_max: maximum number of iteration steps     :param iter_max: maximum number of iteration steps
694     @type iter_max: positive C{int}     :type iter_max: positive ``int``
695     @param sub_iter_max: maximum number of inner iteration steps     :param sub_iter_max: maximum number of inner iteration steps
696     @type sub_iter_max: positive C{int}     :type sub_iter_max: positive ``int``
697     @param atol: absolute tolerance for the solution     :param atol: absolute tolerance for the solution
698     @type atol: positive C{float}     :type atol: positive ``float``
699     @param rtol: relative tolerance for the solution     :param rtol: relative tolerance for the solution
700     @type rtol: positive C{float}     :type rtol: positive ``float``
701     @param gamma: tolerance safety factor for inner iteration     :param gamma: tolerance safety factor for inner iteration
702     @type gamma: positive C{float}, less than 1     :type gamma: positive ``float``, less than 1
703     @param sub_tol_max: upper bound for inner tolerance     :param subtol_max: upper bound for inner tolerance
704     @type sub_tol_max: positive C{float}, less than 1     :type subtol_max: positive ``float``, less than 1
705     @return: an approximation of the solution with the desired accuracy     :return: an approximation of the solution with the desired accuracy
706     @rtype: same type as the initial guess     :rtype: same type as the initial guess
707     """     """
708     lmaxit=iter_max     lmaxit=iter_max
709     if atol<0: raise ValueError,"atol needs to be non-negative."     if atol<0: raise ValueError,"atol needs to be non-negative."
710     if rtol<0: raise ValueError,"rtol needs to be non-negative."     if rtol<0: raise ValueError,"rtol needs to be non-negative."
711     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."
712     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
713     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
714    
715     F=defect(x)     F=defect(x)
716     fnrm=defect.norm(F)     fnrm=defect.norm(F)
717     stop_tol=atol + rtol*fnrm     stop_tol=atol + rtol*fnrm
718     sub_tol=sub_tol_max     subtol=subtol_max
719     if verbose: print "NewtonGMRES: initial residual = %e."%fnrm     if verbose: print "NewtonGMRES: initial residual = %e."%fnrm
720     if verbose: print "             tolerance = %e."%sub_tol     if verbose: print "             tolerance = %e."%subtol
721     iter=1     iter=1
722     #     #
723     # main iteration loop     # main iteration loop
# Line 695  def NewtonGMRES(defect, x, iter_max=100, Line 725  def NewtonGMRES(defect, x, iter_max=100,
725     while not fnrm<=stop_tol:     while not fnrm<=stop_tol:
726              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
727              #              #
728          #   adjust sub_tol_          #   adjust subtol_
729          #          #
730              if iter > 1:              if iter > 1:
731             rat=fnrm/fnrmo             rat=fnrm/fnrmo
732                 sub_tol_old=sub_tol                 subtol_old=subtol
733             sub_tol=gamma*rat**2             subtol=gamma*rat**2
734             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)
735             sub_tol=max(min(sub_tol,sub_tol_max), .5*stop_tol/fnrm)             subtol=max(min(subtol,subtol_max), .5*stop_tol/fnrm)
736          #          #
737          # calculate newton increment xc          # calculate newton increment xc
738              #     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 740  def NewtonGMRES(defect, x, iter_max=100,
740              #     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
741              #              #
742              #              #
743              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
744              try:              try:
745                 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)
746              except MaxIterReached:              except MaxIterReached:
747                 raise MaxIterReached,"maximum number of %s steps reached."%iter_max                 raise MaxIterReached,"maximum number of %s steps reached."%iter_max
748              if sub_iter<0:              if sub_iter<0:
# Line 731  def NewtonGMRES(defect, x, iter_max=100, Line 761  def NewtonGMRES(defect, x, iter_max=100,
761  def __givapp(c,s,vin):  def __givapp(c,s,vin):
762      """      """
763      Applies a sequence of Givens rotations (c,s) recursively to the vector      Applies a sequence of Givens rotations (c,s) recursively to the vector
764      C{vin}      ``vin``
765    
766      @warning: C{vin} is altered.      :warning: ``vin`` is altered.
767      """      """
768      vrot=vin      vrot=vin
769      if isinstance(c,float):      if isinstance(c,float):
# Line 838  def GMRES(r, Aprod, x, bilinearform, ato Line 868  def GMRES(r, Aprod, x, bilinearform, ato
868     """     """
869     Solver for     Solver for
870    
871     M{Ax=b}     *Ax=b*
872    
873     with a general operator A (more details required!).     with a general operator A (more details required!).
874     It uses the generalized minimum residual method (GMRES).     It uses the generalized minimum residual method (GMRES).
875    
876     The iteration is terminated if     The iteration is terminated if
877    
878     M{|r| <= atol+rtol*|r0|}     *|r| <= atol+rtol*|r0|*
879    
880     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
881    
882     M{|r| = sqrt( bilinearform(r,r))}     *|r| = sqrt( bilinearform(r,r))*
883    
884     @param r: initial residual M{r=b-Ax}. C{r} is altered.     :param r: initial residual *r=b-Ax*. ``r`` is altered.
885     @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)
886     @param x: an initial guess for the solution     :param x: an initial guess for the solution
887     @type x: same like C{r}     :type x: same like ``r``
888     @param Aprod: returns the value Ax     :param Aprod: returns the value Ax
889     @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
890                  argument C{x}. The returned object needs to be of the same                  argument ``x``. The returned object needs to be of the same
891                  type like argument C{r}.                  type like argument ``r``.
892     @param bilinearform: inner product C{<x,r>}     :param bilinearform: inner product ``<x,r>``
893     @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
894                         type like argument C{x} and C{r}. The returned value is                         type like argument ``x`` and ``r``. The returned value is
895                         a C{float}.                         a ``float``.
896     @param atol: absolute tolerance     :param atol: absolute tolerance
897     @type atol: non-negative C{float}     :type atol: non-negative ``float``
898     @param rtol: relative tolerance     :param rtol: relative tolerance
899     @type rtol: non-negative C{float}     :type rtol: non-negative ``float``
900     @param iter_max: maximum number of iteration steps     :param iter_max: maximum number of iteration steps
901     @type iter_max: C{int}     :type iter_max: ``int``
902     @param iter_restart: in order to save memory the orthogonalization process     :param iter_restart: in order to save memory the orthogonalization process
903                          is terminated after C{iter_restart} steps and the                          is terminated after ``iter_restart`` steps and the
904                          iteration is restarted.                          iteration is restarted.
905     @type iter_restart: C{int}     :type iter_restart: ``int``
906     @return: the solution approximation and the corresponding residual     :return: the solution approximation and the corresponding residual
907     @rtype: C{tuple}     :rtype: ``tuple``
908     @warning: C{r} and C{x} are altered.     :warning: ``r`` and ``x`` are altered.
909     """     """
910     m=iter_restart     m=iter_restart
911     restarted=False     restarted=False
# Line 1007  def MINRES(r, Aprod, x, Msolve, bilinear Line 1037  def MINRES(r, Aprod, x, Msolve, bilinear
1037      """      """
1038      Solver for      Solver for
1039    
1040      M{Ax=b}      *Ax=b*
1041    
1042      with a symmetric and positive definite operator A (more details required!).      with a symmetric and positive definite operator A (more details required!).
1043      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 1045  def MINRES(r, Aprod, x, Msolve, bilinear
1045    
1046      The iteration is terminated if      The iteration is terminated if
1047    
1048      M{|r| <= atol+rtol*|r0|}      *|r| <= atol+rtol*|r0|*
1049    
1050      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
1051    
1052      M{|r| = sqrt( bilinearform(Msolve(r),r))}      *|r| = sqrt( bilinearform(Msolve(r),r))*
1053    
1054      For details on the preconditioned conjugate gradient method see the book:      For details on the preconditioned conjugate gradient method see the book:
1055    
# Line 1027  def MINRES(r, Aprod, x, Msolve, bilinear Line 1057  def MINRES(r, Aprod, x, Msolve, bilinear
1057      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,
1058      C. Romine, and H. van der Vorst}.      C. Romine, and H. van der Vorst}.
1059    
1060      @param r: initial residual M{r=b-Ax}. C{r} is altered.      :param r: initial residual *r=b-Ax*. ``r`` is altered.
1061      @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)
1062      @param x: an initial guess for the solution      :param x: an initial guess for the solution
1063      @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)
1064      @param Aprod: returns the value Ax      :param Aprod: returns the value Ax
1065      @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
1066                   argument C{x}. The returned object needs to be of the same                   argument ``x``. The returned object needs to be of the same
1067                   type like argument C{r}.                   type like argument ``r``.
1068      @param Msolve: solves Mx=r      :param Msolve: solves Mx=r
1069      @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
1070                    argument C{r}. The returned object needs to be of the same                    argument ``r``. The returned object needs to be of the same
1071                    type like argument C{x}.                    type like argument ``x``.
1072      @param bilinearform: inner product C{<x,r>}      :param bilinearform: inner product ``<x,r>``
1073      @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
1074                          type like argument C{x} and C{r} is. The returned value                          type like argument ``x`` and ``r`` is. The returned value
1075                          is a C{float}.                          is a ``float``.
1076      @param atol: absolute tolerance      :param atol: absolute tolerance
1077      @type atol: non-negative C{float}      :type atol: non-negative ``float``
1078      @param rtol: relative tolerance      :param rtol: relative tolerance
1079      @type rtol: non-negative C{float}      :type rtol: non-negative ``float``
1080      @param iter_max: maximum number of iteration steps      :param iter_max: maximum number of iteration steps
1081      @type iter_max: C{int}      :type iter_max: ``int``
1082      @return: the solution approximation and the corresponding residual      :return: the solution approximation and the corresponding residual
1083      @rtype: C{tuple}      :rtype: ``tuple``
1084      @warning: C{r} and C{x} are altered.      :warning: ``r`` and ``x`` are altered.
1085      """      """
1086      #------------------------------------------------------------------      #------------------------------------------------------------------
1087      # 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 1215  def TFQMR(r, Aprod, x, bilinearform, ato
1215    """    """
1216    Solver for    Solver for
1217    
1218    M{Ax=b}    *Ax=b*
1219    
1220    with a general operator A (more details required!).    with a general operator A (more details required!).
1221    It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).    It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).
1222    
1223    The iteration is terminated if    The iteration is terminated if
1224    
1225    M{|r| <= atol+rtol*|r0|}    *|r| <= atol+rtol*|r0|*
1226    
1227    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
1228    
1229    M{|r| = sqrt( bilinearform(r,r))}    *|r| = sqrt( bilinearform(r,r))*
1230    
1231    @param r: initial residual M{r=b-Ax}. C{r} is altered.    :param r: initial residual *r=b-Ax*. ``r`` is altered.
1232    @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)
1233    @param x: an initial guess for the solution    :param x: an initial guess for the solution
1234    @type x: same like C{r}    :type x: same like ``r``
1235    @param Aprod: returns the value Ax    :param Aprod: returns the value Ax
1236    @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
1237                 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
1238                 like argument C{r}.                 like argument ``r``.
1239    @param bilinearform: inner product C{<x,r>}    :param bilinearform: inner product ``<x,r>``
1240    @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
1241                        type like argument C{x} and C{r}. The returned value is                        type like argument ``x`` and ``r``. The returned value is
1242                        a C{float}.                        a ``float``.
1243    @param atol: absolute tolerance    :param atol: absolute tolerance
1244    @type atol: non-negative C{float}    :type atol: non-negative ``float``
1245    @param rtol: relative tolerance    :param rtol: relative tolerance
1246    @type rtol: non-negative C{float}    :type rtol: non-negative ``float``
1247    @param iter_max: maximum number of iteration steps    :param iter_max: maximum number of iteration steps
1248    @type iter_max: C{int}    :type iter_max: ``int``
1249    @rtype: C{tuple}    :rtype: ``tuple``
1250    @warning: C{r} and C{x} are altered.    :warning: ``r`` and ``x`` are altered.
1251    """    """
1252    u1=0    u1=0
1253    u2=0    u2=0
# Line 1287  def TFQMR(r, Aprod, x, bilinearform, ato Line 1317  def TFQMR(r, Aprod, x, bilinearform, ato
1317    
1318  class ArithmeticTuple(object):  class ArithmeticTuple(object):
1319     """     """
1320     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
1321     ArithmeticTuple and C{a} is a float.     ArithmeticTuple and ``a`` is a float.
1322    
1323     Example of usage::     Example of usage::
1324    
# Line 1302  class ArithmeticTuple(object): Line 1332  class ArithmeticTuple(object):
1332     """     """
1333     def __init__(self,*args):     def __init__(self,*args):
1334         """         """
1335         Initializes object with elements C{args}.         Initializes object with elements ``args``.
1336    
1337         @param args: tuple of objects that support inplace add (x+=y) and         :param args: tuple of objects that support inplace add (x+=y) and
1338                      scaling (x=a*y)                      scaling (x=a*y)
1339         """         """
1340         self.__items=list(args)         self.__items=list(args)
# Line 1313  class ArithmeticTuple(object): Line 1343  class ArithmeticTuple(object):
1343         """         """
1344         Returns the number of items.         Returns the number of items.
1345    
1346         @return: number of items         :return: number of items
1347         @rtype: C{int}         :rtype: ``int``
1348         """         """
1349         return len(self.__items)         return len(self.__items)
1350    
# Line 1322  class ArithmeticTuple(object): Line 1352  class ArithmeticTuple(object):
1352         """         """
1353         Returns item at specified position.         Returns item at specified position.
1354    
1355         @param index: index of item to be returned         :param index: index of item to be returned
1356         @type index: C{int}         :type index: ``int``
1357         @return: item with index C{index}         :return: item with index ``index``
1358         """         """
1359         return self.__items.__getitem__(index)         return self.__items.__getitem__(index)
1360    
1361     def __mul__(self,other):     def __mul__(self,other):
1362         """         """
1363         Scales by C{other} from the right.         Scales by ``other`` from the right.
1364    
1365         @param other: scaling factor         :param other: scaling factor
1366         @type other: C{float}         :type other: ``float``
1367         @return: itemwise self*other         :return: itemwise self*other
1368         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1369         """         """
1370         out=[]         out=[]
1371         try:         try:
# Line 1349  class ArithmeticTuple(object): Line 1379  class ArithmeticTuple(object):
1379    
1380     def __rmul__(self,other):     def __rmul__(self,other):
1381         """         """
1382         Scales by C{other} from the left.         Scales by ``other`` from the left.
1383    
1384         @param other: scaling factor         :param other: scaling factor
1385         @type other: C{float}         :type other: ``float``
1386         @return: itemwise other*self         :return: itemwise other*self
1387         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1388         """         """
1389         out=[]         out=[]
1390         try:         try:
# Line 1368  class ArithmeticTuple(object): Line 1398  class ArithmeticTuple(object):
1398    
1399     def __div__(self,other):     def __div__(self,other):
1400         """         """
1401         Scales by (1/C{other}) from the right.         Scales by (1/``other``) from the right.
1402    
1403         @param other: scaling factor         :param other: scaling factor
1404         @type other: C{float}         :type other: ``float``
1405         @return: itemwise self/other         :return: itemwise self/other
1406         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1407         """         """
1408         return self*(1/other)         return self*(1/other)
1409    
1410     def __rdiv__(self,other):     def __rdiv__(self,other):
1411         """         """
1412         Scales by (1/C{other}) from the left.         Scales by (1/``other``) from the left.
1413    
1414         @param other: scaling factor         :param other: scaling factor
1415         @type other: C{float}         :type other: ``float``
1416         @return: itemwise other/self         :return: itemwise other/self
1417         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1418         """         """
1419         out=[]         out=[]
1420         try:         try:
# Line 1398  class ArithmeticTuple(object): Line 1428  class ArithmeticTuple(object):
1428    
1429     def __iadd__(self,other):     def __iadd__(self,other):
1430         """         """
1431         Inplace addition of C{other} to self.         Inplace addition of ``other`` to self.
1432    
1433         @param other: increment         :param other: increment
1434         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1435         """         """
1436         if len(self) != len(other):         if len(self) != len(other):
1437             raise ValueError,"tuple lengths must match."             raise ValueError,"tuple lengths must match."
# Line 1411  class ArithmeticTuple(object): Line 1441  class ArithmeticTuple(object):
1441    
1442     def __add__(self,other):     def __add__(self,other):
1443         """         """
1444         Adds C{other} to self.         Adds ``other`` to self.
1445    
1446         @param other: increment         :param other: increment
1447         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1448         """         """
1449         out=[]         out=[]
1450         try:         try:
# Line 1428  class ArithmeticTuple(object): Line 1458  class ArithmeticTuple(object):
1458    
1459     def __sub__(self,other):     def __sub__(self,other):
1460         """         """
1461         Subtracts C{other} from self.         Subtracts ``other`` from self.
1462    
1463         @param other: decrement         :param other: decrement
1464         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1465         """         """
1466         out=[]         out=[]
1467         try:         try:
# Line 1445  class ArithmeticTuple(object): Line 1475  class ArithmeticTuple(object):
1475    
1476     def __isub__(self,other):     def __isub__(self,other):
1477         """         """
1478         Inplace subtraction of C{other} from self.         Inplace subtraction of ``other`` from self.
1479    
1480         @param other: decrement         :param other: decrement
1481         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1482         """         """
1483         if len(self) != len(other):         if len(self) != len(other):
1484             raise ValueError,"tuple length must match."             raise ValueError,"tuple length must match."
# Line 1471  class HomogeneousSaddlePointProblem(obje Line 1501  class HomogeneousSaddlePointProblem(obje
1501        This class provides a framework for solving linear homogeneous saddle        This class provides a framework for solving linear homogeneous saddle
1502        point problems of the form::        point problems of the form::
1503    
1504            M{Av+B^*p=f}            *Av+B^*p=f*
1505            M{Bv     =0}            *Bv     =0*
1506    
1507        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
1508        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*.
1509          *A* may depend weakly on *v* and *p*.
1510        """        """
1511        def __init__(self, adaptSubTolerance=True, **kwargs):        def __init__(self, **kwargs):
1512      """      """
1513      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}  
1514      """      """
1515            self.resetControlParameters()
1516          self.setTolerance()          self.setTolerance()
1517          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
1518      self.__adaptSubTolerance=adaptSubTolerance        def resetControlParameters(self, K_p=1., K_v=1., rtol_max=0.01, rtol_min = 1.e-7, chi_max=0.5, reduction_factor=0.3, theta = 0.1):
1519        #=============================================================           """
1520        def initialize(self):           sets a control parameter
1521          """  
1522          Initializes the problem (overwrite).           :param K_p: initial value for constant to adjust pressure tolerance
1523          """           :type K_p: ``float``
1524          pass           :param K_v: initial value for constant to adjust velocity tolerance
1525             :type K_v: ``float``
1526             :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1527             :type rtol_max: ``float``
1528             :param chi_max: maximum tolerable converegence rate.
1529             :type chi_max: ``float``
1530             :param reduction_factor: reduction factor for adjustment factors.
1531             :type reduction_factor: ``float``
1532             """
1533             self.setControlParameter(K_p, K_v, rtol_max, rtol_min, chi_max, reduction_factor, theta)
1534    
1535          def setControlParameter(self,K_p=None, K_v=None, rtol_max=None, rtol_min=None, chi_max=None, reduction_factor=None, theta=None):
1536             """
1537             sets a control parameter
1538    
1539    
1540             :param K_p: initial value for constant to adjust pressure tolerance
1541             :type K_p: ``float``
1542             :param K_v: initial value for constant to adjust velocity tolerance
1543             :type K_v: ``float``
1544             :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1545             :type rtol_max: ``float``
1546             :param chi_max: maximum tolerable converegence rate.
1547             :type chi_max: ``float``
1548             :type reduction_factor: ``float``
1549             """
1550             if not K_p == None:
1551                if K_p<1:
1552                   raise ValueError,"K_p need to be greater or equal to 1."
1553             else:
1554                K_p=self.__K_p
1555    
1556             if not K_v == None:
1557                if K_v<1:
1558                   raise ValueError,"K_v need to be greater or equal to 1."
1559             else:
1560                K_v=self.__K_v
1561    
1562             if not rtol_max == None:
1563                if rtol_max<=0 or rtol_max>=1:
1564                   raise ValueError,"rtol_max needs to be positive and less than 1."
1565             else:
1566                rtol_max=self.__rtol_max
1567    
1568             if not rtol_min == None:
1569                if rtol_min<=0 or rtol_min>=1:
1570                   raise ValueError,"rtol_min needs to be positive and less than 1."
1571             else:
1572                rtol_min=self.__rtol_min
1573    
1574             if not chi_max == None:
1575                if chi_max<=0 or chi_max>=1:
1576                   raise ValueError,"chi_max needs to be positive and less than 1."
1577             else:
1578                chi_max = self.__chi_max
1579    
1580             if not reduction_factor == None:
1581                if reduction_factor<=0 or reduction_factor>1:
1582                   raise ValueError,"reduction_factor need to be between zero and one."
1583             else:
1584                reduction_factor=self.__reduction_factor
1585    
1586             if not theta == None:
1587                if theta<=0 or theta>1:
1588                   raise ValueError,"theta need to be between zero and one."
1589             else:
1590                theta=self.__theta
1591    
1592             if rtol_min>=rtol_max:
1593                 raise ValueError,"rtol_max = %e needs to be greater than rtol_min = %e"%(rtol_max,rtol_min)
1594             self.__chi_max = chi_max
1595             self.__rtol_max = rtol_max
1596             self.__K_p = K_p
1597             self.__K_v = K_v
1598             self.__reduction_factor = reduction_factor
1599             self.__theta = theta
1600             self.__rtol_min=rtol_min
1601    
1602          #=============================================================
1603        def inner_pBv(self,p,Bv):        def inner_pBv(self,p,Bv):
1604           """           """
1605           Returns inner product of element p and Bv (overwrite).           Returns inner product of element p and Bv (overwrite).
1606    
1607           @param p: a pressure increment           :param p: a pressure increment
1608           @param v: a residual           :param Bv: a residual
1609           @return: inner product of element p and Bv           :return: inner product of element p and Bv
1610           @rtype: C{float}           :rtype: ``float``
1611           @note: used if PCG is applied.           :note: used if PCG is applied.
1612           """           """
1613           raise NotImplementedError,"no inner product for p and Bv implemented."           raise NotImplementedError,"no inner product for p and Bv implemented."
1614    
# Line 1510  class HomogeneousSaddlePointProblem(obje Line 1616  class HomogeneousSaddlePointProblem(obje
1616           """           """
1617           Returns inner product of p0 and p1 (overwrite).           Returns inner product of p0 and p1 (overwrite).
1618    
1619           @param p0: a pressure           :param p0: a pressure
1620           @param p1: a pressure           :param p1: a pressure
1621           @return: inner product of p0 and p1           :return: inner product of p0 and p1
1622           @rtype: C{float}           :rtype: ``float``
1623           """           """
1624           raise NotImplementedError,"no inner product for p implemented."           raise NotImplementedError,"no inner product for p implemented."
1625        
# Line 1521  class HomogeneousSaddlePointProblem(obje Line 1627  class HomogeneousSaddlePointProblem(obje
1627           """           """
1628           Returns the norm of v (overwrite).           Returns the norm of v (overwrite).
1629    
1630           @param v: a velovity           :param v: a velovity
1631           @return: norm of v           :return: norm of v
1632           @rtype: non-negative C{float}           :rtype: non-negative ``float``
1633           """           """
1634           raise NotImplementedError,"no norm of v implemented."           raise NotImplementedError,"no norm of v implemented."
1635        def getV(self, p, v0):        def getDV(self, p, v, tol):
1636           """           """
1637           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)
1638    
1639           @param p: a pressure           :param p: pressure
1640           @param v0: a initial guess for the value v to return.           :param v: pressure
1641           @return: v given as M{v= A^{-1} (f-B^*p)}           :return: dv given as *dv= A^{-1} (f-A v-B^*p)*
1642             :note: Only *A* may depend on *v* and *p*
1643           """           """
1644           raise NotImplementedError,"no v calculation implemented."           raise NotImplementedError,"no dv calculation implemented."
1645    
1646                    
1647        def Bv(self,v):        def Bv(self,v, tol):
1648          """          """
1649          Returns Bv (overwrite).          Returns Bv with accuracy `tol` (overwrite)
1650    
1651          @rtype: equal to the type of p          :rtype: equal to the type of p
1652          @note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1653          """          """
1654          raise NotImplementedError, "no operator B implemented."          raise NotImplementedError, "no operator B implemented."
1655    
# Line 1550  class HomogeneousSaddlePointProblem(obje Line 1657  class HomogeneousSaddlePointProblem(obje
1657          """          """
1658          Returns the norm of Bv (overwrite).          Returns the norm of Bv (overwrite).
1659    
1660          @rtype: equal to the type of p          :rtype: equal to the type of p
1661          @note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1662          """          """
1663          raise NotImplementedError, "no norm of Bv implemented."          raise NotImplementedError, "no norm of Bv implemented."
1664    
1665        def solve_AinvBt(self,p):        def solve_AinvBt(self,dp, tol):
1666           """           """
1667           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}           Solves *A dv=B^*dp* with accuracy `tol`
          (overwrite).  
1668    
1669           @param p: a pressure increment           :param dp: a pressure increment
1670           @return: the solution of M{Av=B^*p}           :return: the solution of *A dv=B^*dp*
1671           @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.
1672           """           """
1673           raise NotImplementedError,"no operator A implemented."           raise NotImplementedError,"no operator A implemented."
1674    
1675        def solve_prec(self,Bv):        def solve_prec(self,Bv, tol):
1676           """           """
1677           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).  
1678    
1679           @rtype: equal to the type of p           :rtype: equal to the type of p
1680           @note: boundary conditions on p should be zero!           :note: boundary conditions on p should be zero!
1681           """           """
1682           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  
1683        #=============================================================        #=============================================================
1684        def __Aprod_PCG(self,p):        def __Aprod_PCG(self,dp):
1685            dv=self.solve_AinvBt(p)            dv=self.solve_AinvBt(dp, self.__subtol)
1686            return ArithmeticTuple(dv,self.Bv(dv))            return ArithmeticTuple(dv,self.Bv(dv, self.__subtol))
1687    
1688        def __inner_PCG(self,p,r):        def __inner_PCG(self,p,r):
1689           return self.inner_pBv(p,r[1])           return self.inner_pBv(p,r[1])
1690    
1691        def __Msolve_PCG(self,r):        def __Msolve_PCG(self,r):
1692            return self.solve_prec(r[1])            return self.solve_prec(r[1], self.__subtol)
1693        #=============================================================        #=============================================================
1694        def __Aprod_GMRES(self,p):        def __Aprod_GMRES(self,p):
1695            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)
1696        def __inner_GMRES(self,p0,p1):        def __inner_GMRES(self,p0,p1):
1697           return self.inner_p(p0,p1)           return self.inner_p(p0,p1)
1698    
1699        #=============================================================        #=============================================================
1700        def norm_p(self,p):        def norm_p(self,p):
1701            """            """
1702            calculates the norm of C{p}            calculates the norm of ``p``
1703                        
1704            @param p: a pressure            :param p: a pressure
1705            @return: the norm of C{p} using the inner product for pressure            :return: the norm of ``p`` using the inner product for pressure
1706            @rtype: C{float}            :rtype: ``float``
1707            """            """
1708            f=self.inner_p(p,p)            f=self.inner_p(p,p)
1709            if f<0: raise ValueError,"negative pressure norm."            if f<0: raise ValueError,"negative pressure norm."
1710            return math.sqrt(f)            return math.sqrt(f)
       def adaptSubTolerance(self):  
       """  
       Returns True if tolerance adaption for subproblem is choosen.  
       """  
           self.__adaptSubTolerance  
1711                
1712        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):
1713           """           """
1714           Solves the saddle point problem using initial guesses v and p.           Solves the saddle point problem using initial guesses v and p.
1715    
1716           @param v: initial guess for velocity           :param v: initial guess for velocity
1717           @param p: initial guess for pressure           :param p: initial guess for pressure
1718           @type v: L{Data}           :type v: `Data`
1719           @type p: L{Data}           :type p: `Data`
1720           @param usePCG: indicates the usage of the PCG rather than GMRES scheme.           :param usePCG: indicates the usage of the PCG rather than GMRES scheme.
1721           @param max_iter: maximum number of iteration steps per correction           :param max_iter: maximum number of iteration steps per correction
1722                            attempt                            attempt
1723           @param verbose: if True, shows information on the progress of the           :param verbose: if True, shows information on the progress of the
1724                           saddlepoint problem solver.                           saddlepoint problem solver.
1725           @param iter_restart: restart the iteration after C{iter_restart} steps           :param iter_restart: restart the iteration after ``iter_restart`` steps
1726                                (only used if useUzaw=False)                                (only used if useUzaw=False)
1727           @type usePCG: C{bool}           :type usePCG: ``bool``
1728           @type max_iter: C{int}           :type max_iter: ``int``
1729           @type verbose: C{bool}           :type verbose: ``bool``
1730           @type iter_restart: C{int}           :type iter_restart: ``int``
1731           @rtype: C{tuple} of L{Data} objects           :rtype: ``tuple`` of `Data` objects
1732             :note: typically this method is overwritten by a subclass. It provides a wrapper for the ``_solve`` method.
1733             """
1734             return self._solve(v=v,p=p,max_iter=max_iter,verbose=verbose, usePCG=usePCG, iter_restart=iter_restart, max_correction_steps=max_correction_steps)
1735    
1736          def _solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1737             """
1738             see `_solve` method.
1739           """           """
1740           self.verbose=verbose           self.verbose=verbose
1741           rtol=self.getTolerance()           rtol=self.getTolerance()
1742           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
1743       if self.adaptSubTolerance(): self.setSubProblemTolerance()  
1744             K_p=self.__K_p
1745             K_v=self.__K_v
1746           correction_step=0           correction_step=0
1747           converged=False           converged=False
1748             chi=None
1749             eps=None
1750    
1751             if self.verbose: print "HomogeneousSaddlePointProblem: start iteration: rtol= %e, atol=%e"%(rtol, atol)
1752           while not converged:           while not converged:
1753                # calculate velocity for current pressure:  
1754                v=self.getV(p,v)               # get tolerance for velecity increment:
1755                Bv=self.Bv(v)               if chi == None:
1756                norm_v=self.norm_v(v)                  rtol_v=self.__rtol_max
1757                norm_Bv=self.norm_Bv(Bv)               else:
1758                ATOL=norm_v*rtol+atol                  rtol_v=min(chi/K_v,self.__rtol_max)
1759                if self.verbose: print "HomogeneousSaddlePointProblem: norm v= %e, norm_Bv= %e, tolerance = %e."%(norm_v, norm_Bv,ATOL)               rtol_v=max(rtol_v, self.__rtol_min)
1760                if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."               if self.verbose: print "HomogeneousSaddlePointProblem: step %s: rtol_v= %e"%(correction_step,rtol_v)
1761                if norm_Bv <= ATOL:               # get velocity increment:
1762                   converged=True               dv1=self.getDV(p,v,rtol_v)
1763                else:               v1=v+dv1
1764                   correction_step+=1               Bv1=self.Bv(v1, rtol_v)
1765                   if correction_step>max_correction_steps:               norm_Bv1=self.norm_Bv(Bv1)
1766                        raise CorrectionFailed,"Given up after %d correction steps."%correction_step               norm_dv1=self.norm_v(dv1)
1767                   dp=self.solve_prec(Bv)               if self.verbose: print "HomogeneousSaddlePointProblem: step %s: norm_Bv1 = %e, norm_dv1 = %e"%(correction_step, norm_Bv1, norm_dv1)
1768                   if usePCG:               if norm_dv1*self.__theta < norm_Bv1:
1769                     norm2=self.inner_pBv(dp,Bv)                  # get tolerance for pressure increment:
1770                     if norm2<0: raise ValueError,"negative PCG norm."                  large_Bv1=True
1771                     norm2=math.sqrt(norm2)                  if chi == None or eps == None:
1772                   else:                     rtol_p=self.__rtol_max
1773                     norm2=self.norm_p(dp)                  else:
1774                   ATOL_ITER=ATOL/norm_Bv*norm2*0.5                     rtol_p=min(chi**2*eps/K_p/norm_Bv1, self.__rtol_max)
1775                   if self.verbose: print "HomogeneousSaddlePointProblem: tolerance for solver: %e"%ATOL_ITER                  self.__subtol=max(rtol_p**2, self.__rtol_min)
1776                   if usePCG:                  if self.verbose: print "HomogeneousSaddlePointProblem: step %s: rtol_p= %e"%(correction_step,rtol_p)
1777                         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)                  # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1
1778                    if usePCG:
1779                        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)
1780                        v2=r[0]
1781                        Bv2=r[1]
1782                    else:
1783                        # don't use!!!!
1784                        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)
1785                        dv2=self.solve_AinvBt(dp, self.__subtol)
1786                        v2=v1-dv2
1787                        Bv2=self.Bv(v2, self.__subtol)
1788                    p2=p+dp
1789                 else:
1790                    large_Bv1=False
1791                    v2=v1
1792                    p2=p
1793                 # update business:
1794                 norm_dv2=self.norm_v(v2-v)
1795                 norm_v2=self.norm_v(v2)
1796                 if self.verbose: print "HomogeneousSaddlePointProblem: step %s: v2 = %e, norm_dv2 = %e"%(correction_step, norm_v2, self.norm_v(v2-v))
1797                 eps, eps_old = max(norm_Bv1, norm_dv2), eps
1798                 if eps_old == None:
1799                      chi, chi_old = None, chi
1800                 else:
1801                      chi, chi_old = min(eps/ eps_old, self.__chi_max), chi
1802                 if eps != None:
1803                     if chi !=None:
1804                        if self.verbose: print "HomogeneousSaddlePointProblem: step %s: convergence rate = %e, correction = %e"%(correction_step,chi, eps)
1805                   else:                   else:
1806                         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)                      if self.verbose: print "HomogeneousSaddlePointProblem: step %s: correction = %e"%(correction_step, eps)
1807           if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached."               if eps <= rtol*norm_v2+atol :
1808                     converged = True
1809                 else:
1810                     if correction_step>=max_correction_steps:
1811                          raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1812                     if chi_old!=None:
1813                        K_p=max(1,self.__reduction_factor*K_p,(chi-chi_old)/chi_old**2*K_p)
1814                        K_v=max(1,self.__reduction_factor*K_v,(chi-chi_old)/chi_old**2*K_p)
1815                        if self.verbose: print "HomogeneousSaddlePointProblem: step %s: new adjustment factor K = %e"%(correction_step,K_p)
1816                 correction_step+=1
1817                 v,p =v2, p2
1818             if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step
1819       return v,p       return v,p
   
1820        #========================================================================        #========================================================================
1821        def setTolerance(self,tolerance=1.e-4):        def setTolerance(self,tolerance=1.e-4):
1822           """           """
1823           Sets the relative tolerance for (v,p).           Sets the relative tolerance for (v,p).
1824    
1825           @param tolerance: tolerance to be used           :param tolerance: tolerance to be used
1826           @type tolerance: non-negative C{float}           :type tolerance: non-negative ``float``
1827           """           """
1828           if tolerance<0:           if tolerance<0:
1829               raise ValueError,"tolerance must be positive."               raise ValueError,"tolerance must be positive."
# Line 1689  class HomogeneousSaddlePointProblem(obje Line 1833  class HomogeneousSaddlePointProblem(obje
1833           """           """
1834           Returns the relative tolerance.           Returns the relative tolerance.
1835    
1836           @return: relative tolerance           :return: relative tolerance
1837           @rtype: C{float}           :rtype: ``float``
1838           """           """
1839           return self.__rtol           return self.__rtol
1840    
# Line 1698  class HomogeneousSaddlePointProblem(obje Line 1842  class HomogeneousSaddlePointProblem(obje
1842           """           """
1843           Sets the absolute tolerance.           Sets the absolute tolerance.
1844    
1845           @param tolerance: tolerance to be used           :param tolerance: tolerance to be used
1846           @type tolerance: non-negative C{float}           :type tolerance: non-negative ``float``
1847           """           """
1848           if tolerance<0:           if tolerance<0:
1849               raise ValueError,"tolerance must be non-negative."               raise ValueError,"tolerance must be non-negative."
# Line 1709  class HomogeneousSaddlePointProblem(obje Line 1853  class HomogeneousSaddlePointProblem(obje
1853           """           """
1854           Returns the absolute tolerance.           Returns the absolute tolerance.
1855    
1856           @return: absolute tolerance           :return: absolute tolerance
1857           @rtype: C{float}           :rtype: ``float``
1858           """           """
1859           return self.__atol           return self.__atol
1860    
       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)  
1861    
1862  def MaskFromBoundaryTag(domain,*tags):  def MaskFromBoundaryTag(domain,*tags):
1863     """     """
# Line 1730  def MaskFromBoundaryTag(domain,*tags): Line 1866  def MaskFromBoundaryTag(domain,*tags):
1866    
1867     Usage: m=MaskFromBoundaryTag(domain, "left", "right")     Usage: m=MaskFromBoundaryTag(domain, "left", "right")
1868    
1869     @param domain: domain to be used     :param domain: domain to be used
1870     @type domain: L{escript.Domain}     :type domain: `escript.Domain`
1871     @param tags: boundary tags     :param tags: boundary tags
1872     @type tags: C{str}     :type tags: ``str``
1873     @return: a mask which marks samples that are touching the boundary tagged     :return: a mask which marks samples that are touching the boundary tagged
1874              by any of the given tags              by any of the given tags
1875     @rtype: L{escript.Data} of rank 0     :rtype: `escript.Data` of rank 0
1876     """     """
1877     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1878     d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))     d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
# Line 1751  def MaskFromTag(domain,*tags): Line 1887  def MaskFromTag(domain,*tags):
1887    
1888     Usage: m=MaskFromTag(domain, "ham")     Usage: m=MaskFromTag(domain, "ham")
1889    
1890     @param domain: domain to be used     :param domain: domain to be used
1891     @type domain: L{escript.Domain}     :type domain: `escript.Domain`
1892     @param tags: boundary tags     :param tags: boundary tags
1893     @type tags: C{str}     :type tags: ``str``
1894     @return: a mask which marks samples that are touching the boundary tagged     :return: a mask which marks samples that are touching the boundary tagged
1895              by any of the given tags              by any of the given tags
1896     @rtype: L{escript.Data} of rank 0     :rtype: `escript.Data` of rank 0
1897     """     """
1898     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1899     d=escript.Scalar(0.,escript.Function(domain))     d=escript.Scalar(0.,escript.Function(domain))

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

  ViewVC Help
Powered by ViewVC 1.1.26