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

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

  ViewVC Help
Powered by ViewVC 1.1.26