/[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 3911 by jfenwick, Thu Jun 14 01:01:03 2012 UTC
# Line 1  Line 1 
1    
2  ########################################################  ########################################################
3  #  #
4  # Copyright (c) 2003-2009 by University of Queensland  # Copyright (c) 2003-2012 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-2012 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"
40    
41    
42  import escript  from . import escript
43  import linearPDEs  from . import linearPDEs
44  import numpy  import numpy
45  import util  from . 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 "p" in kwargs:
70              self.__p=kwargs["p"]              self.__p=kwargs["p"]
71       else:       else:
72              self.__p=1              self.__p=1
73       if kwargs.has_key("time"):       if "time" in kwargs:
74              self.__t=kwargs["time"]              self.__t=kwargs["time"]
75       else:       else:
76              self.__t=0.              self.__t=0.
# 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:
302                 raise ValueError,"coefficient D is undefined"                 raise ValueError("coefficient D is undefined")
303              D=escript.Data(self.__D,self.__function_space)              D=escript.Data(self.__D,self.__function_space)
304              if D.getRank()>1:              if D.getRank()>1:
305                 raise ValueError,"coefficient D must have rank 0 or 1"                 raise ValueError("coefficient D must have rank 0 or 1")
306              if self.__Y==None:              if self.__Y==None:
307                 self.__u=escript.Data(0.,D.getShape(),self.__function_space)                 self.__u=escript.Data(0.,D.getShape(),self.__function_space)
308              else:              else:
# 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
347             except TypeError:             except TypeError:
348               iterative=False               iterative=False
349           xxx=self.__function_space.getX()
350         if iterative:         if iterative:
351             self.__id=[]             self.__id=[]
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 425  class Locator: Line 430  class Locator:
430                  return out                  return out
431          else:          else:
432             return data             return data
433              
434    #     def setValue(self, data, v):
435    #       """
436    #       Sets the value of the ``data`` at the Locator.
437    #       """
438    #       data.expand()   # Need to ensure that this is done globally
439    #       if isinstance(data, escript.Data):
440    #    id=self.getId()
441    #    if isinstance(id, list):
442    #      for i in id:
443    #         data._setTupleForGlobalDataPoint(i[1], i[0], v)
444    #    else:
445    #      data._setTupleForGlobalDataPoint(id[1], id[0], v)
446    #       else:
447    #    raise TypeError, "setValue: Invalid argument type."
448    
449    
450    def getInfLocator(arg):
451        """
452        Return a Locator for a point with the inf value over all arg.
453        """
454        if not isinstance(arg, escript.Data):
455           raise TypeError("getInfLocator: Unknown argument type.")
456        a_inf=util.inf(arg)
457        loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords
458        x=arg.getFunctionSpace().getX()
459        x_min=x.getTupleForGlobalDataPoint(*loc)
460        return Locator(arg.getFunctionSpace(),x_min)
461    
462    def getSupLocator(arg):
463        """
464        Return a Locator for a point with the sup value over all arg.
465        """
466        if not isinstance(arg, escript.Data):
467           raise TypeError("getInfLocator: Unknown argument type.")
468        a_inf=util.sup(arg)
469        loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords
470        x=arg.getFunctionSpace().getX()
471        x_min=x.getTupleForGlobalDataPoint(*loc)
472        return Locator(arg.getFunctionSpace(),x_min)
473        
474    
475  class SolverSchemeException(Exception):  class SolverSchemeException(Exception):
476     """     """
# Line 467  def PCG(r, Aprod, x, Msolve, bilinearfor Line 513  def PCG(r, Aprod, x, Msolve, bilinearfor
513     """     """
514     Solver for     Solver for
515    
516     M{Ax=b}     *Ax=b*
517    
518     with a symmetric and positive definite operator A (more details required!).     with a symmetric and positive definite operator A (more details required!).
519     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 521  def PCG(r, Aprod, x, Msolve, bilinearfor
521    
522     The iteration is terminated if     The iteration is terminated if
523    
524     M{|r| <= atol+rtol*|r0|}     *|r| <= atol+rtol*|r0|*
525    
526     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
527    
528     M{|r| = sqrt( bilinearform(Msolve(r),r))}     *|r| = sqrt( bilinearform(Msolve(r),r))*
529    
530     For details on the preconditioned conjugate gradient method see the book:     For details on the preconditioned conjugate gradient method see the book:
531    
# Line 487  def PCG(r, Aprod, x, Msolve, bilinearfor Line 533  def PCG(r, Aprod, x, Msolve, bilinearfor
533     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,
534     C. Romine, and H. van der Vorst}.     C. Romine, and H. van der Vorst}.
535    
536     @param r: initial residual M{r=b-Ax}. C{r} is altered.     :param r: initial residual *r=b-Ax*. ``r`` is altered.
537     @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)
538     @param x: an initial guess for the solution     :param x: an initial guess for the solution
539     @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)
540     @param Aprod: returns the value Ax     :param Aprod: returns the value Ax
541     @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
542                  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
543                  like argument C{r}.                  like argument ``r``.
544     @param Msolve: solves Mx=r     :param Msolve: solves Mx=r
545     @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
546                   argument C{r}. The returned object needs to be of the same                   argument ``r``. The returned object needs to be of the same
547                   type like argument C{x}.                   type like argument ``x``.
548     @param bilinearform: inner product C{<x,r>}     :param bilinearform: inner product ``<x,r>``
549     @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
550                         type like argument C{x} and C{r} is. The returned value                         type like argument ``x`` and ``r`` is. The returned value
551                         is a C{float}.                         is a ``float``.
552     @param atol: absolute tolerance     :param atol: absolute tolerance
553     @type atol: non-negative C{float}     :type atol: non-negative ``float``
554     @param rtol: relative tolerance     :param rtol: relative tolerance
555     @type rtol: non-negative C{float}     :type rtol: non-negative ``float``
556     @param iter_max: maximum number of iteration steps     :param iter_max: maximum number of iteration steps
557     @type iter_max: C{int}     :type iter_max: ``int``
558     @return: the solution approximation and the corresponding residual     :return: the solution approximation and the corresponding residual
559     @rtype: C{tuple}     :rtype: ``tuple``
560     @warning: C{r} and C{x} are altered.     :warning: ``r`` and ``x`` are altered.
561     """     """
562     iter=0     iter=0
563     rhat=Msolve(r)     rhat=Msolve(r)
564     d = rhat     d = rhat
565     rhat_dot_r = bilinearform(rhat, r)     rhat_dot_r = bilinearform(rhat, r)
566     if rhat_dot_r<0: raise NegativeNorm,"negative norm."     if rhat_dot_r<0: raise NegativeNorm("negative norm.")
567     norm_r0=math.sqrt(rhat_dot_r)     norm_r0=math.sqrt(rhat_dot_r)
568     atol2=atol+rtol*norm_r0     atol2=atol+rtol*norm_r0
569     if atol2<=0:     if atol2<=0:
570        raise ValueError,"Non-positive tolarance."        raise ValueError("Non-positive tolarance.")
571     atol2=max(atol2, 100. * util.EPSILON * norm_r0)     atol2=max(atol2, 100. * util.EPSILON * norm_r0)
572    
573     if verbose: print "PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_r0, atol2)     if verbose: print(("PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_r0, atol2)))
574    
575    
576     while not math.sqrt(rhat_dot_r) <= atol2:     while not math.sqrt(rhat_dot_r) <= atol2:
577         iter+=1         iter+=1
578         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)
579    
580         q=Aprod(d)         q=Aprod(d)
581         alpha = rhat_dot_r / bilinearform(d, q)         alpha = rhat_dot_r / bilinearform(d, q)
582         x += alpha * d         x += alpha * d
583         if isinstance(q,ArithmeticTuple):         if isinstance(q,ArithmeticTuple):
584         r += q * (-alpha)      # Doing it the other way calls the float64.__mul__ not AT.__rmul__            r += q * (-alpha)      # Doing it the other way calls the float64.__mul__ not AT.__rmul__
585         else:         else:
586             r += (-alpha) * q             r += (-alpha) * q
587         rhat=Msolve(r)         rhat=Msolve(r)
# Line 545  def PCG(r, Aprod, x, Msolve, bilinearfor Line 591  def PCG(r, Aprod, x, Msolve, bilinearfor
591         d=rhat         d=rhat
592    
593         rhat_dot_r = rhat_dot_r_new         rhat_dot_r = rhat_dot_r_new
594         if rhat_dot_r<0: raise NegativeNorm,"negative norm."         if rhat_dot_r<0: raise NegativeNorm("negative norm.")
595         if verbose: print "PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))         if verbose: print(("PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))))
596     if verbose: print "PCG: tolerance reached after %s steps."%iter     if verbose: print(("PCG: tolerance reached after %s steps."%iter))
597     return x,r,math.sqrt(rhat_dot_r)     return x,r,math.sqrt(rhat_dot_r)
598    
599  class Defect(object):  class Defect(object):
# Line 564  class Defect(object): Line 610  class Defect(object):
610          """          """
611          Returns the inner product of x0 and x1          Returns the inner product of x0 and x1
612    
613          @param x0: value for x0          :param x0: value for x0
614          @param x1: value for x1          :param x1: value for x1
615          @return: the inner product of x0 and x1          :return: the inner product of x0 and x1
616          @rtype: C{float}          :rtype: ``float``
617          """          """
618          return 0          return 0
619    
620      def norm(self,x):      def norm(self,x):
621          """          """
622          Returns the norm of argument C{x}.          Returns the norm of argument ``x``.
623    
624          @param x: a value          :param x: a value
625          @return: norm of argument x          :return: norm of argument x
626          @rtype: C{float}          :rtype: ``float``
627          @note: by default C{sqrt(self.bilinearform(x,x)} is returned.          :note: by default ``sqrt(self.bilinearform(x,x)`` is returned.
628          """          """
629          s=self.bilinearform(x,x)          s=self.bilinearform(x,x)
630          if s<0: raise NegativeNorm,"negative norm."          if s<0: raise NegativeNorm("negative norm.")
631          return math.sqrt(s)          return math.sqrt(s)
632    
633      def eval(self,x):      def eval(self,x):
634          """          """
635          Returns the value F of a given C{x}.          Returns the value F of a given ``x``.
636    
637          @param x: value for which the defect C{F} is evaluated          :param x: value for which the defect ``F`` is evaluated
638          @return: value of the defect at C{x}          :return: value of the defect at ``x``
639          """          """
640          return 0          return 0
641    
642      def __call__(self,x):      def __call__(self,x):
643          return self.eval(x)          return self.eval(x)
644    
645      def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):      def setDerivativeIncrementLength(self,inc=1000.*math.sqrt(util.EPSILON)):
646          """          """
647          Sets the relative length of the increment used to approximate the          Sets the relative length of the increment used to approximate the
648          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
649          direction of v with x as a starting point.          direction of v with x as a starting point.
650    
651          @param inc: relative increment length          :param inc: relative increment length
652          @type inc: positive C{float}          :type inc: positive ``float``
653          """          """
654          if inc<=0: raise ValueError,"positive increment required."          if inc<=0: raise ValueError("positive increment required.")
655          self.__inc=inc          self.__inc=inc
656    
657      def getDerivativeIncrementLength(self):      def getDerivativeIncrementLength(self):
658          """          """
659          Returns the relative increment length used to approximate the          Returns the relative increment length used to approximate the
660          derivative of the defect.          derivative of the defect.
661          @return: value of the defect at C{x}          :return: value of the defect at ``x``
662          @rtype: positive C{float}          :rtype: positive ``float``
663          """          """
664          return self.__inc          return self.__inc
665    
666      def derivative(self, F0, x0, v, v_is_normalised=True):      def derivative(self, F0, x0, v, v_is_normalised=True):
667          """          """
668          Returns the directional derivative at C{x0} in the direction of C{v}.          Returns the directional derivative at ``x0`` in the direction of ``v``.
669    
670          @param F0: value of this defect at x0          :param F0: value of this defect at x0
671          @param x0: value at which derivative is calculated          :param x0: value at which derivative is calculated
672          @param v: direction          :param v: direction
673          @param v_is_normalised: True to indicate that C{v} is nomalized          :param v_is_normalised: True to indicate that ``v`` is nomalized
674                                  (self.norm(v)=0)                                  (self.norm(v)=0)
675          @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``
676          @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
677                 used but this method maybe overwritten to use exact evaluation.                 used but this method maybe overwritten to use exact evaluation.
678          """          """
679          normx=self.norm(x0)          normx=self.norm(x0)
# Line 645  class Defect(object): Line 691  class Defect(object):
691          return (F1-F0)/epsnew          return (F1-F0)/epsnew
692    
693  ######################################  ######################################
694  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):
695     """     """
696     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
697     criterion:     criterion:
698    
699     M{norm(F(x) <= atol + rtol * norm(F(x0)}     *norm(F(x) <= atol + rtol * norm(F(x0)*
700    
701     where M{x0} is the initial guess.     where *x0* is the initial guess.
702    
703     @param defect: object defining the function M{F}. C{defect.norm} defines the     :param defect: object defining the function *F*. ``defect.norm`` defines the
704                    M{norm} used in the stopping criterion.                    *norm* used in the stopping criterion.
705     @type defect: L{Defect}     :type defect: `Defect`
706     @param x: initial guess for the solution, C{x} is altered.     :param x: initial guess for the solution, ``x`` is altered.
707     @type x: any object type allowing basic operations such as     :type x: any object type allowing basic operations such as
708              C{numpy.ndarray}, L{Data}              ``numpy.ndarray``, `Data`
709     @param iter_max: maximum number of iteration steps     :param iter_max: maximum number of iteration steps
710     @type iter_max: positive C{int}     :type iter_max: positive ``int``
711     @param sub_iter_max: maximum number of inner iteration steps     :param sub_iter_max: maximum number of inner iteration steps
712     @type sub_iter_max: positive C{int}     :type sub_iter_max: positive ``int``
713     @param atol: absolute tolerance for the solution     :param atol: absolute tolerance for the solution
714     @type atol: positive C{float}     :type atol: positive ``float``
715     @param rtol: relative tolerance for the solution     :param rtol: relative tolerance for the solution
716     @type rtol: positive C{float}     :type rtol: positive ``float``
717     @param gamma: tolerance safety factor for inner iteration     :param gamma: tolerance safety factor for inner iteration
718     @type gamma: positive C{float}, less than 1     :type gamma: positive ``float``, less than 1
719     @param sub_tol_max: upper bound for inner tolerance     :param subtol_max: upper bound for inner tolerance
720     @type sub_tol_max: positive C{float}, less than 1     :type subtol_max: positive ``float``, less than 1
721     @return: an approximation of the solution with the desired accuracy     :return: an approximation of the solution with the desired accuracy
722     @rtype: same type as the initial guess     :rtype: same type as the initial guess
723     """     """
724     lmaxit=iter_max     lmaxit=iter_max
725     if atol<0: raise ValueError,"atol needs to be non-negative."     if atol<0: raise ValueError("atol needs to be non-negative.")
726     if rtol<0: raise ValueError,"rtol needs to be non-negative."     if rtol<0: raise ValueError("rtol needs to be non-negative.")
727     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.")
728     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)
729     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)
730    
731     F=defect(x)     F=defect(x)
732     fnrm=defect.norm(F)     fnrm=defect.norm(F)
733     stop_tol=atol + rtol*fnrm     stop_tol=atol + rtol*fnrm
734     sub_tol=sub_tol_max     subtol=subtol_max
735     if verbose: print "NewtonGMRES: initial residual = %e."%fnrm     if verbose: print(("NewtonGMRES: initial residual = %e."%fnrm))
736     if verbose: print "             tolerance = %e."%sub_tol     if verbose: print(("             tolerance = %e."%subtol))
737     iter=1     iter=1
738     #     #
739     # main iteration loop     # main iteration loop
740     #     #
741     while not fnrm<=stop_tol:     while not fnrm<=stop_tol:
742              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)
743              #              #
744          #   adjust sub_tol_          #   adjust subtol_
745          #          #
746              if iter > 1:              if iter > 1:
747             rat=fnrm/fnrmo                 rat=fnrm/fnrmo
748                 sub_tol_old=sub_tol                 subtol_old=subtol
749             sub_tol=gamma*rat**2                 subtol=gamma*rat**2
750             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)
751             sub_tol=max(min(sub_tol,sub_tol_max), .5*stop_tol/fnrm)                 subtol=max(min(subtol,subtol_max), .5*stop_tol/fnrm)
752          #          #
753          # calculate newton increment xc          # calculate newton increment xc
754              #     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 756  def NewtonGMRES(defect, x, iter_max=100,
756              #     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
757              #              #
758              #              #
759              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))
760              try:              try:
761                 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)
762              except MaxIterReached:              except MaxIterReached:
763                 raise MaxIterReached,"maximum number of %s steps reached."%iter_max                 raise MaxIterReached("maximum number of %s steps reached."%iter_max)
764              if sub_iter<0:              if sub_iter<0:
765                 iter+=sub_iter_max                 iter+=sub_iter_max
766              else:              else:
767                 iter+=sub_iter                 iter+=sub_iter
768              # ====              # ====
769          x+=xc              x+=xc
770              F=defect(x)              F=defect(x)
771          iter+=1              iter+=1
772              fnrmo, fnrm=fnrm, defect.norm(F)              fnrmo, fnrm=fnrm, defect.norm(F)
773              if verbose: print "             step %s: residual %e."%(iter,fnrm)              if verbose: print(("             step %s: residual %e."%(iter,fnrm)))
774     if verbose: print "NewtonGMRES: completed after %s steps."%iter     if verbose: print(("NewtonGMRES: completed after %s steps."%iter))
775     return x     return x
776    
777  def __givapp(c,s,vin):  def __givapp(c,s,vin):
778      """      """
779      Applies a sequence of Givens rotations (c,s) recursively to the vector      Applies a sequence of Givens rotations (c,s) recursively to the vector
780      C{vin}      ``vin``
781    
782      @warning: C{vin} is altered.      :warning: ``vin`` is altered.
783      """      """
784      vrot=vin      vrot=vin
785      if isinstance(c,float):      if isinstance(c,float):
# Line 741  def __givapp(c,s,vin): Line 787  def __givapp(c,s,vin):
787      else:      else:
788          for i in range(len(c)):          for i in range(len(c)):
789              w1=c[i]*vrot[i]-s[i]*vrot[i+1]              w1=c[i]*vrot[i]-s[i]*vrot[i+1]
790          w2=s[i]*vrot[i]+c[i]*vrot[i+1]              w2=s[i]*vrot[i]+c[i]*vrot[i+1]
791              vrot[i]=w1              vrot[i]=w1
792              vrot[i+1]=w2              vrot[i+1]=w2
793      return vrot      return vrot
# Line 761  def __FDGMRES(F0, defect, x0, atol, iter Line 807  def __FDGMRES(F0, defect, x0, atol, iter
807     iter=0     iter=0
808     while rho > atol and iter<iter_restart-1:     while rho > atol and iter<iter_restart-1:
809          if iter  >= iter_max:          if iter  >= iter_max:
810              raise MaxIterReached,"maximum number of %s steps reached."%iter_max              raise MaxIterReached("maximum number of %s steps reached."%iter_max)
811    
812          p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)          p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
813          v.append(p)          v.append(p)
# Line 823  def __FDGMRES(F0, defect, x0, atol, iter Line 869  def __FDGMRES(F0, defect, x0, atol, iter
869            i=i-1            i=i-1
870       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
871       for i in range(iter-1):       for i in range(iter-1):
872      xhat += v[i]*y[i]         xhat += v[i]*y[i]
873     else :     else :
874        xhat=v[0] * 0        xhat=v[0] * 0
875    
# Line 838  def GMRES(r, Aprod, x, bilinearform, ato Line 884  def GMRES(r, Aprod, x, bilinearform, ato
884     """     """
885     Solver for     Solver for
886    
887     M{Ax=b}     *Ax=b*
888    
889     with a general operator A (more details required!).     with a general operator A (more details required!).
890     It uses the generalized minimum residual method (GMRES).     It uses the generalized minimum residual method (GMRES).
891    
892     The iteration is terminated if     The iteration is terminated if
893    
894     M{|r| <= atol+rtol*|r0|}     *|r| <= atol+rtol*|r0|*
895    
896     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
897    
898     M{|r| = sqrt( bilinearform(r,r))}     *|r| = sqrt( bilinearform(r,r))*
899    
900     @param r: initial residual M{r=b-Ax}. C{r} is altered.     :param r: initial residual *r=b-Ax*. ``r`` is altered.
901     @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)
902     @param x: an initial guess for the solution     :param x: an initial guess for the solution
903     @type x: same like C{r}     :type x: same like ``r``
904     @param Aprod: returns the value Ax     :param Aprod: returns the value Ax
905     @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
906                  argument C{x}. The returned object needs to be of the same                  argument ``x``. The returned object needs to be of the same
907                  type like argument C{r}.                  type like argument ``r``.
908     @param bilinearform: inner product C{<x,r>}     :param bilinearform: inner product ``<x,r>``
909     @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
910                         type like argument C{x} and C{r}. The returned value is                         type like argument ``x`` and ``r``. The returned value is
911                         a C{float}.                         a ``float``.
912     @param atol: absolute tolerance     :param atol: absolute tolerance
913     @type atol: non-negative C{float}     :type atol: non-negative ``float``
914     @param rtol: relative tolerance     :param rtol: relative tolerance
915     @type rtol: non-negative C{float}     :type rtol: non-negative ``float``
916     @param iter_max: maximum number of iteration steps     :param iter_max: maximum number of iteration steps
917     @type iter_max: C{int}     :type iter_max: ``int``
918     @param iter_restart: in order to save memory the orthogonalization process     :param iter_restart: in order to save memory the orthogonalization process
919                          is terminated after C{iter_restart} steps and the                          is terminated after ``iter_restart`` steps and the
920                          iteration is restarted.                          iteration is restarted.
921     @type iter_restart: C{int}     :type iter_restart: ``int``
922     @return: the solution approximation and the corresponding residual     :return: the solution approximation and the corresponding residual
923     @rtype: C{tuple}     :rtype: ``tuple``
924     @warning: C{r} and C{x} are altered.     :warning: ``r`` and ``x`` are altered.
925     """     """
926     m=iter_restart     m=iter_restart
927     restarted=False     restarted=False
928     iter=0     iter=0
929     if rtol>0:     if rtol>0:
930        r_dot_r = bilinearform(r, r)        r_dot_r = bilinearform(r, r)
931        if r_dot_r<0: raise NegativeNorm,"negative norm."        if r_dot_r<0: raise NegativeNorm("negative norm.")
932        atol2=atol+rtol*math.sqrt(r_dot_r)        atol2=atol+rtol*math.sqrt(r_dot_r)
933        if verbose: print "GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)        if verbose: print(("GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)))
934     else:     else:
935        atol2=atol        atol2=atol
936        if verbose: print "GMRES: absolute tolerance = %e"%atol2        if verbose: print(("GMRES: absolute tolerance = %e"%atol2))
937     if atol2<=0:     if atol2<=0:
938        raise ValueError,"Non-positive tolarance."        raise ValueError("Non-positive tolarance.")
939    
940     while True:     while True:
941        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)
942        if restarted:        if restarted:
943           r2 = r-Aprod(x-x2)           r2 = r-Aprod(x-x2)
944        else:        else:
# Line 901  def GMRES(r, Aprod, x, bilinearform, ato Line 947  def GMRES(r, Aprod, x, bilinearform, ato
947        x,stopped=_GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose,P_R=P_R)        x,stopped=_GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose,P_R=P_R)
948        iter+=iter_restart        iter+=iter_restart
949        if stopped: break        if stopped: break
950        if verbose: print "GMRES: restart."        if verbose: print("GMRES: restart.")
951        restarted=True        restarted=True
952     if verbose: print "GMRES: tolerance has been reached."     if verbose: print("GMRES: tolerance has been reached.")
953     return x     return x
954    
955  def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False, P_R=None):  def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False, P_R=None):
# Line 916  def _GMRESm(r, Aprod, x, bilinearform, a Line 962  def _GMRESm(r, Aprod, x, bilinearform, a
962     v=[]     v=[]
963    
964     r_dot_r = bilinearform(r, r)     r_dot_r = bilinearform(r, r)
965     if r_dot_r<0: raise NegativeNorm,"negative norm."     if r_dot_r<0: raise NegativeNorm("negative norm.")
966     rho=math.sqrt(r_dot_r)     rho=math.sqrt(r_dot_r)
967    
968     v.append(r/rho)     v.append(r/rho)
969     g[0]=rho     g[0]=rho
970    
971     if verbose: print "GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)     if verbose: print(("GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)))
972     while not (rho<=atol or iter==iter_restart):     while not (rho<=atol or iter==iter_restart):
973    
974      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)
975    
976          if P_R!=None:          if P_R!=None:
977              p=Aprod(P_R(v[iter]))              p=Aprod(P_R(v[iter]))
978          else:          else:
979          p=Aprod(v[iter])              p=Aprod(v[iter])
980      v.append(p)          v.append(p)
981    
982      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))          v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
983    
984  # Modified Gram-Schmidt  # Modified Gram-Schmidt
985      for j in range(iter+1):          for j in range(iter+1):
986        h[j,iter]=bilinearform(v[j],v[iter+1])            h[j,iter]=bilinearform(v[j],v[iter+1])
987        v[iter+1]-=h[j,iter]*v[j]            v[iter+1]-=h[j,iter]*v[j]
988    
989      h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))          h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
990      v_norm2=h[iter+1,iter]          v_norm2=h[iter+1,iter]
991    
992  # Reorthogonalize if needed  # Reorthogonalize if needed
993      if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)          if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
994       for j in range(iter+1):           for j in range(iter+1):
995          hr=bilinearform(v[j],v[iter+1])              hr=bilinearform(v[j],v[iter+1])
996              h[j,iter]=h[j,iter]+hr              h[j,iter]=h[j,iter]+hr
997              v[iter+1] -= hr*v[j]              v[iter+1] -= hr*v[j]
998    
999       v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))           v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
1000       h[iter+1,iter]=v_norm2           h[iter+1,iter]=v_norm2
1001    
1002  #   watch out for happy breakdown  #   watch out for happy breakdown
1003          if not v_norm2 == 0:          if not v_norm2 == 0:
1004           v[iter+1]=v[iter+1]/h[iter+1,iter]           v[iter+1]=v[iter+1]/h[iter+1,iter]
1005    
1006  #   Form and store the information for the new Givens rotation  #   Form and store the information for the new Givens rotation
1007      if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])          if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])
1008      mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])          mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
1009    
1010      if mu!=0 :          if mu!=0 :
1011          c[iter]=h[iter,iter]/mu                  c[iter]=h[iter,iter]/mu
1012          s[iter]=-h[iter+1,iter]/mu                  s[iter]=-h[iter+1,iter]/mu
1013          h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]                  h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
1014          h[iter+1,iter]=0.0                  h[iter+1,iter]=0.0
1015                  gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])                  gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
1016                  g[iter]=gg[0]                  g[iter]=gg[0]
1017                  g[iter+1]=gg[1]                  g[iter+1]=gg[1]
1018  # Update the residual norm  # Update the residual norm
1019    
1020          rho=abs(g[iter+1])          rho=abs(g[iter+1])
1021          if verbose: print "GMRES: iteration step %s: residual %e"%(iter,rho)          if verbose: print(("GMRES: iteration step %s: residual %e"%(iter,rho)))
1022      iter+=1          iter+=1
1023    
1024  # At this point either iter > iter_max or rho < tol.  # At this point either iter > iter_max or rho < tol.
1025  # It's time to compute x and leave.  # It's time to compute x and leave.
1026    
1027     if verbose: print "GMRES: iteration stopped after %s step."%iter     if verbose: print(("GMRES: iteration stopped after %s step."%iter))
1028     if iter > 0 :     if iter > 0 :
1029       y=numpy.zeros(iter,numpy.float64)       y=numpy.zeros(iter,numpy.float64)
1030       y[iter-1] = g[iter-1] / h[iter-1,iter-1]       y[iter-1] = g[iter-1] / h[iter-1,iter-1]
# Line 989  def _GMRESm(r, Aprod, x, bilinearform, a Line 1035  def _GMRESm(r, Aprod, x, bilinearform, a
1035            i=i-1            i=i-1
1036       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
1037       for i in range(iter-1):       for i in range(iter-1):
1038      xhat += v[i]*y[i]         xhat += v[i]*y[i]
1039     else:     else:
1040       xhat=v[0] * 0       xhat=v[0] * 0
1041     if P_R!=None:     if P_R!=None:
# Line 1007  def MINRES(r, Aprod, x, Msolve, bilinear Line 1053  def MINRES(r, Aprod, x, Msolve, bilinear
1053      """      """
1054      Solver for      Solver for
1055    
1056      M{Ax=b}      *Ax=b*
1057    
1058      with a symmetric and positive definite operator A (more details required!).      with a symmetric and positive definite operator A (more details required!).
1059      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 1061  def MINRES(r, Aprod, x, Msolve, bilinear
1061    
1062      The iteration is terminated if      The iteration is terminated if
1063    
1064      M{|r| <= atol+rtol*|r0|}      *|r| <= atol+rtol*|r0|*
1065    
1066      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
1067    
1068      M{|r| = sqrt( bilinearform(Msolve(r),r))}      *|r| = sqrt( bilinearform(Msolve(r),r))*
1069    
1070      For details on the preconditioned conjugate gradient method see the book:      For details on the preconditioned conjugate gradient method see the book:
1071    
# Line 1027  def MINRES(r, Aprod, x, Msolve, bilinear Line 1073  def MINRES(r, Aprod, x, Msolve, bilinear
1073      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,
1074      C. Romine, and H. van der Vorst}.      C. Romine, and H. van der Vorst}.
1075    
1076      @param r: initial residual M{r=b-Ax}. C{r} is altered.      :param r: initial residual *r=b-Ax*. ``r`` is altered.
1077      @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)
1078      @param x: an initial guess for the solution      :param x: an initial guess for the solution
1079      @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)
1080      @param Aprod: returns the value Ax      :param Aprod: returns the value Ax
1081      @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
1082                   argument C{x}. The returned object needs to be of the same                   argument ``x``. The returned object needs to be of the same
1083                   type like argument C{r}.                   type like argument ``r``.
1084      @param Msolve: solves Mx=r      :param Msolve: solves Mx=r
1085      @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
1086                    argument C{r}. The returned object needs to be of the same                    argument ``r``. The returned object needs to be of the same
1087                    type like argument C{x}.                    type like argument ``x``.
1088      @param bilinearform: inner product C{<x,r>}      :param bilinearform: inner product ``<x,r>``
1089      @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
1090                          type like argument C{x} and C{r} is. The returned value                          type like argument ``x`` and ``r`` is. The returned value
1091                          is a C{float}.                          is a ``float``.
1092      @param atol: absolute tolerance      :param atol: absolute tolerance
1093      @type atol: non-negative C{float}      :type atol: non-negative ``float``
1094      @param rtol: relative tolerance      :param rtol: relative tolerance
1095      @type rtol: non-negative C{float}      :type rtol: non-negative ``float``
1096      @param iter_max: maximum number of iteration steps      :param iter_max: maximum number of iteration steps
1097      @type iter_max: C{int}      :type iter_max: ``int``
1098      @return: the solution approximation and the corresponding residual      :return: the solution approximation and the corresponding residual
1099      @rtype: C{tuple}      :rtype: ``tuple``
1100      @warning: C{r} and C{x} are altered.      :warning: ``r`` and ``x`` are altered.
1101      """      """
1102      #------------------------------------------------------------------      #------------------------------------------------------------------
1103      # Set up y and v for the first Lanczos vector v1.      # Set up y and v for the first Lanczos vector v1.
# Line 1062  def MINRES(r, Aprod, x, Msolve, bilinear Line 1108  def MINRES(r, Aprod, x, Msolve, bilinear
1108      y = Msolve(r)      y = Msolve(r)
1109      beta1 = bilinearform(y,r)      beta1 = bilinearform(y,r)
1110    
1111      if beta1< 0: raise NegativeNorm,"negative norm."      if beta1< 0: raise NegativeNorm("negative norm.")
1112    
1113      #  If r = 0 exactly, stop with x      #  If r = 0 exactly, stop with x
1114      if beta1==0: return x      if beta1==0: return x
# Line 1097  def MINRES(r, Aprod, x, Msolve, bilinear Line 1143  def MINRES(r, Aprod, x, Msolve, bilinear
1143      # --------------------------------------------------------------------      # --------------------------------------------------------------------
1144      while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL      while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL
1145    
1146      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)
1147          iter    = iter  +  1          iter    = iter  +  1
1148    
1149          #-----------------------------------------------------------------          #-----------------------------------------------------------------
# Line 1127  def MINRES(r, Aprod, x, Msolve, bilinear Line 1173  def MINRES(r, Aprod, x, Msolve, bilinear
1173          y = Msolve(r2)          y = Msolve(r2)
1174          oldb   = beta                           # oldb = betak          oldb   = beta                           # oldb = betak
1175          beta   = bilinearform(y,r2)             # beta = betak+1^2          beta   = bilinearform(y,r2)             # beta = betak+1^2
1176          if beta < 0: raise NegativeNorm,"negative norm."          if beta < 0: raise NegativeNorm("negative norm.")
1177    
1178          beta   = math.sqrt( beta )          beta   = math.sqrt( beta )
1179          tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta          tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
# Line 1185  def TFQMR(r, Aprod, x, bilinearform, ato Line 1231  def TFQMR(r, Aprod, x, bilinearform, ato
1231    """    """
1232    Solver for    Solver for
1233    
1234    M{Ax=b}    *Ax=b*
1235    
1236    with a general operator A (more details required!).    with a general operator A (more details required!).
1237    It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).    It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).
1238    
1239    The iteration is terminated if    The iteration is terminated if
1240    
1241    M{|r| <= atol+rtol*|r0|}    *|r| <= atol+rtol*|r0|*
1242    
1243    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
1244    
1245    M{|r| = sqrt( bilinearform(r,r))}    *|r| = sqrt( bilinearform(r,r))*
1246    
1247    @param r: initial residual M{r=b-Ax}. C{r} is altered.    :param r: initial residual *r=b-Ax*. ``r`` is altered.
1248    @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)
1249    @param x: an initial guess for the solution    :param x: an initial guess for the solution
1250    @type x: same like C{r}    :type x: same like ``r``
1251    @param Aprod: returns the value Ax    :param Aprod: returns the value Ax
1252    @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
1253                 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
1254                 like argument C{r}.                 like argument ``r``.
1255    @param bilinearform: inner product C{<x,r>}    :param bilinearform: inner product ``<x,r>``
1256    @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
1257                        type like argument C{x} and C{r}. The returned value is                        type like argument ``x`` and ``r``. The returned value is
1258                        a C{float}.                        a ``float``.
1259    @param atol: absolute tolerance    :param atol: absolute tolerance
1260    @type atol: non-negative C{float}    :type atol: non-negative ``float``
1261    @param rtol: relative tolerance    :param rtol: relative tolerance
1262    @type rtol: non-negative C{float}    :type rtol: non-negative ``float``
1263    @param iter_max: maximum number of iteration steps    :param iter_max: maximum number of iteration steps
1264    @type iter_max: C{int}    :type iter_max: ``int``
1265    @rtype: C{tuple}    :rtype: ``tuple``
1266    @warning: C{r} and C{x} are altered.    :warning: ``r`` and ``x`` are altered.
1267    """    """
1268    u1=0    u1=0
1269    u2=0    u2=0
# Line 1234  def TFQMR(r, Aprod, x, bilinearform, ato Line 1280  def TFQMR(r, Aprod, x, bilinearform, ato
1280    theta = 0.0;    theta = 0.0;
1281    eta = 0.0;    eta = 0.0;
1282    rho=bilinearform(r,r)    rho=bilinearform(r,r)
1283    if rho < 0: raise NegativeNorm,"negative norm."    if rho < 0: raise NegativeNorm("negative norm.")
1284    tau = math.sqrt(rho)    tau = math.sqrt(rho)
1285    norm_r0=tau    norm_r0=tau
1286    while tau>atol+rtol*norm_r0:    while tau>atol+rtol*norm_r0:
1287      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)
1288    
1289      sigma = bilinearform(r,v)      sigma = bilinearform(r,v)
1290      if sigma == 0.0: raise IterationBreakDown,'TFQMR breakdown, sigma=0'      if sigma == 0.0: raise IterationBreakDown('TFQMR breakdown, sigma=0')
1291    
1292      alpha = rho / sigma      alpha = rho / sigma
1293    
# Line 1269  def TFQMR(r, Aprod, x, bilinearform, ato Line 1315  def TFQMR(r, Aprod, x, bilinearform, ato
1315  #  #
1316  #  Try to terminate the iteration at each pass through the loop  #  Try to terminate the iteration at each pass through the loop
1317  #  #
1318      if rho == 0.0: raise IterationBreakDown,'TFQMR breakdown, rho=0'      if rho == 0.0: raise IterationBreakDown('TFQMR breakdown, rho=0')
1319    
1320      rhon = bilinearform(r,w)      rhon = bilinearform(r,w)
1321      beta = rhon / rho;      beta = rhon / rho;
# Line 1287  def TFQMR(r, Aprod, x, bilinearform, ato Line 1333  def TFQMR(r, Aprod, x, bilinearform, ato
1333    
1334  class ArithmeticTuple(object):  class ArithmeticTuple(object):
1335     """     """
1336     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
1337     ArithmeticTuple and C{a} is a float.     ArithmeticTuple and ``a`` is a float.
1338    
1339     Example of usage::     Example of usage::
1340    
# Line 1302  class ArithmeticTuple(object): Line 1348  class ArithmeticTuple(object):
1348     """     """
1349     def __init__(self,*args):     def __init__(self,*args):
1350         """         """
1351         Initializes object with elements C{args}.         Initializes object with elements ``args``.
1352    
1353         @param args: tuple of objects that support inplace add (x+=y) and         :param args: tuple of objects that support inplace add (x+=y) and
1354                      scaling (x=a*y)                      scaling (x=a*y)
1355         """         """
1356         self.__items=list(args)         self.__items=list(args)
# Line 1313  class ArithmeticTuple(object): Line 1359  class ArithmeticTuple(object):
1359         """         """
1360         Returns the number of items.         Returns the number of items.
1361    
1362         @return: number of items         :return: number of items
1363         @rtype: C{int}         :rtype: ``int``
1364         """         """
1365         return len(self.__items)         return len(self.__items)
1366    
# Line 1322  class ArithmeticTuple(object): Line 1368  class ArithmeticTuple(object):
1368         """         """
1369         Returns item at specified position.         Returns item at specified position.
1370    
1371         @param index: index of item to be returned         :param index: index of item to be returned
1372         @type index: C{int}         :type index: ``int``
1373         @return: item with index C{index}         :return: item with index ``index``
1374         """         """
1375         return self.__items.__getitem__(index)         return self.__items.__getitem__(index)
1376    
1377     def __mul__(self,other):     def __mul__(self,other):
1378         """         """
1379         Scales by C{other} from the right.         Scales by ``other`` from the right.
1380    
1381         @param other: scaling factor         :param other: scaling factor
1382         @type other: C{float}         :type other: ``float``
1383         @return: itemwise self*other         :return: itemwise self*other
1384         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1385         """         """
1386         out=[]         out=[]
1387         try:         try:
1388             l=len(other)             l=len(other)
1389             if l!=len(self):             if l!=len(self):
1390                 raise ValueError,"length of arguments don't match."                 raise ValueError("length of arguments don't match.")
1391             for i in range(l): out.append(self[i]*other[i])             for i in range(l): out.append(self[i]*other[i])
1392         except TypeError:         except TypeError:
1393             for i in range(len(self)): out.append(self[i]*other)             for i in range(len(self)): out.append(self[i]*other)
# Line 1349  class ArithmeticTuple(object): Line 1395  class ArithmeticTuple(object):
1395    
1396     def __rmul__(self,other):     def __rmul__(self,other):
1397         """         """
1398         Scales by C{other} from the left.         Scales by ``other`` from the left.
1399    
1400         @param other: scaling factor         :param other: scaling factor
1401         @type other: C{float}         :type other: ``float``
1402         @return: itemwise other*self         :return: itemwise other*self
1403         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1404         """         """
1405         out=[]         out=[]
1406         try:         try:
1407             l=len(other)             l=len(other)
1408             if l!=len(self):             if l!=len(self):
1409                 raise ValueError,"length of arguments don't match."                 raise ValueError("length of arguments don't match.")
1410             for i in range(l): out.append(other[i]*self[i])             for i in range(l): out.append(other[i]*self[i])
1411         except TypeError:         except TypeError:
1412             for i in range(len(self)): out.append(other*self[i])             for i in range(len(self)): out.append(other*self[i])
# Line 1368  class ArithmeticTuple(object): Line 1414  class ArithmeticTuple(object):
1414    
1415     def __div__(self,other):     def __div__(self,other):
1416         """         """
1417         Scales by (1/C{other}) from the right.         Scales by (1/``other``) from the right.
1418    
1419         @param other: scaling factor         :param other: scaling factor
1420         @type other: C{float}         :type other: ``float``
1421         @return: itemwise self/other         :return: itemwise self/other
1422         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1423         """         """
1424         return self*(1/other)         return self*(1/other)
1425    
1426     def __rdiv__(self,other):     def __rdiv__(self,other):
1427         """         """
1428         Scales by (1/C{other}) from the left.         Scales by (1/``other``) from the left.
1429    
1430         @param other: scaling factor         :param other: scaling factor
1431         @type other: C{float}         :type other: ``float``
1432         @return: itemwise other/self         :return: itemwise other/self
1433         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1434         """         """
1435         out=[]         out=[]
1436         try:         try:
1437             l=len(other)             l=len(other)
1438             if l!=len(self):             if l!=len(self):
1439                 raise ValueError,"length of arguments don't match."                 raise ValueError("length of arguments don't match.")
1440             for i in range(l): out.append(other[i]/self[i])             for i in range(l): out.append(other[i]/self[i])
1441         except TypeError:         except TypeError:
1442             for i in range(len(self)): out.append(other/self[i])             for i in range(len(self)): out.append(other/self[i])
# Line 1398  class ArithmeticTuple(object): Line 1444  class ArithmeticTuple(object):
1444    
1445     def __iadd__(self,other):     def __iadd__(self,other):
1446         """         """
1447         Inplace addition of C{other} to self.         Inplace addition of ``other`` to self.
1448    
1449         @param other: increment         :param other: increment
1450         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1451         """         """
1452         if len(self) != len(other):         if len(self) != len(other):
1453             raise ValueError,"tuple lengths must match."             raise ValueError("tuple lengths must match.")
1454         for i in range(len(self)):         for i in range(len(self)):
1455             self.__items[i]+=other[i]             self.__items[i]+=other[i]
1456         return self         return self
1457    
1458     def __add__(self,other):     def __add__(self,other):
1459         """         """
1460         Adds C{other} to self.         Adds ``other`` to self.
1461    
1462         @param other: increment         :param other: increment
1463         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1464         """         """
1465         out=[]         out=[]
1466         try:         try:
1467             l=len(other)             l=len(other)
1468             if l!=len(self):             if l!=len(self):
1469                 raise ValueError,"length of arguments don't match."                 raise ValueError("length of arguments don't match.")
1470             for i in range(l): out.append(self[i]+other[i])             for i in range(l): out.append(self[i]+other[i])
1471         except TypeError:         except TypeError:
1472             for i in range(len(self)): out.append(self[i]+other)             for i in range(len(self)): out.append(self[i]+other)
# Line 1428  class ArithmeticTuple(object): Line 1474  class ArithmeticTuple(object):
1474    
1475     def __sub__(self,other):     def __sub__(self,other):
1476         """         """
1477         Subtracts C{other} from self.         Subtracts ``other`` from self.
1478    
1479         @param other: decrement         :param other: decrement
1480         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1481         """         """
1482         out=[]         out=[]
1483         try:         try:
1484             l=len(other)             l=len(other)
1485             if l!=len(self):             if l!=len(self):
1486                 raise ValueError,"length of arguments don't match."                 raise ValueError("length of arguments don't match.")
1487             for i in range(l): out.append(self[i]-other[i])             for i in range(l): out.append(self[i]-other[i])
1488         except TypeError:         except TypeError:
1489             for i in range(len(self)): out.append(self[i]-other)             for i in range(len(self)): out.append(self[i]-other)
# Line 1445  class ArithmeticTuple(object): Line 1491  class ArithmeticTuple(object):
1491    
1492     def __isub__(self,other):     def __isub__(self,other):
1493         """         """
1494         Inplace subtraction of C{other} from self.         Inplace subtraction of ``other`` from self.
1495    
1496         @param other: decrement         :param other: decrement
1497         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1498         """         """
1499         if len(self) != len(other):         if len(self) != len(other):
1500             raise ValueError,"tuple length must match."             raise ValueError("tuple length must match.")
1501         for i in range(len(self)):         for i in range(len(self)):
1502             self.__items[i]-=other[i]             self.__items[i]-=other[i]
1503         return self         return self
# Line 1471  class HomogeneousSaddlePointProblem(obje Line 1517  class HomogeneousSaddlePointProblem(obje
1517        This class provides a framework for solving linear homogeneous saddle        This class provides a framework for solving linear homogeneous saddle
1518        point problems of the form::        point problems of the form::
1519    
1520            M{Av+B^*p=f}            *Av+B^*p=f*
1521            M{Bv     =0}            *Bv     =0*
1522    
1523        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
1524        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*.
1525          *A* may depend weakly on *v* and *p*.
1526        """        """
1527        def __init__(self, adaptSubTolerance=True, **kwargs):        def __init__(self, **kwargs):
     """  
     initializes the saddle point problem  
       
     @param adaptSubTolerance: If True the tolerance for subproblem is set automatically.  
     @type adaptSubTolerance: C{bool}  
     """  
         self.setTolerance()  
         self.setAbsoluteTolerance()  
     self.__adaptSubTolerance=adaptSubTolerance  
       #=============================================================  
       def initialize(self):  
1528          """          """
1529          Initializes the problem (overwrite).          initializes the saddle point problem
1530          """          """
1531          pass          self.resetControlParameters()
1532            self.setTolerance()
1533            self.setAbsoluteTolerance()
1534          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):
1535             """
1536             sets a control parameter
1537    
1538             :param K_p: initial value for constant to adjust pressure tolerance
1539             :type K_p: ``float``
1540             :param K_v: initial value for constant to adjust velocity tolerance
1541             :type K_v: ``float``
1542             :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1543             :type rtol_max: ``float``
1544             :param chi_max: maximum tolerable converegence rate.
1545             :type chi_max: ``float``
1546             :param reduction_factor: reduction factor for adjustment factors.
1547             :type reduction_factor: ``float``
1548             """
1549             self.setControlParameter(K_p, K_v, rtol_max, rtol_min, chi_max, reduction_factor, theta)
1550    
1551          def setControlParameter(self,K_p=None, K_v=None, rtol_max=None, rtol_min=None, chi_max=None, reduction_factor=None, theta=None):
1552             """
1553             sets a control parameter
1554    
1555    
1556             :param K_p: initial value for constant to adjust pressure tolerance
1557             :type K_p: ``float``
1558             :param K_v: initial value for constant to adjust velocity tolerance
1559             :type K_v: ``float``
1560             :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1561             :type rtol_max: ``float``
1562             :param chi_max: maximum tolerable converegence rate.
1563             :type chi_max: ``float``
1564             :type reduction_factor: ``float``
1565             """
1566             if not K_p == None:
1567                if K_p<1:
1568                   raise ValueError("K_p need to be greater or equal to 1.")
1569             else:
1570                K_p=self.__K_p
1571    
1572             if not K_v == None:
1573                if K_v<1:
1574                   raise ValueError("K_v need to be greater or equal to 1.")
1575             else:
1576                K_v=self.__K_v
1577    
1578             if not rtol_max == None:
1579                if rtol_max<=0 or rtol_max>=1:
1580                   raise ValueError("rtol_max needs to be positive and less than 1.")
1581             else:
1582                rtol_max=self.__rtol_max
1583    
1584             if not rtol_min == None:
1585                if rtol_min<=0 or rtol_min>=1:
1586                   raise ValueError("rtol_min needs to be positive and less than 1.")
1587             else:
1588                rtol_min=self.__rtol_min
1589    
1590             if not chi_max == None:
1591                if chi_max<=0 or chi_max>=1:
1592                   raise ValueError("chi_max needs to be positive and less than 1.")
1593             else:
1594                chi_max = self.__chi_max
1595    
1596             if not reduction_factor == None:
1597                if reduction_factor<=0 or reduction_factor>1:
1598                   raise ValueError("reduction_factor need to be between zero and one.")
1599             else:
1600                reduction_factor=self.__reduction_factor
1601    
1602             if not theta == None:
1603                if theta<=0 or theta>1:
1604                   raise ValueError("theta need to be between zero and one.")
1605             else:
1606                theta=self.__theta
1607    
1608             if rtol_min>=rtol_max:
1609                 raise ValueError("rtol_max = %e needs to be greater than rtol_min = %e"%(rtol_max,rtol_min))
1610             self.__chi_max = chi_max
1611             self.__rtol_max = rtol_max
1612             self.__K_p = K_p
1613             self.__K_v = K_v
1614             self.__reduction_factor = reduction_factor
1615             self.__theta = theta
1616             self.__rtol_min=rtol_min
1617    
1618          #=============================================================
1619        def inner_pBv(self,p,Bv):        def inner_pBv(self,p,Bv):
1620           """           """
1621           Returns inner product of element p and Bv (overwrite).           Returns inner product of element p and Bv (overwrite).
1622    
1623           @param p: a pressure increment           :param p: a pressure increment
1624           @param v: a residual           :param Bv: a residual
1625           @return: inner product of element p and Bv           :return: inner product of element p and Bv
1626           @rtype: C{float}           :rtype: ``float``
1627           @note: used if PCG is applied.           :note: used if PCG is applied.
1628           """           """
1629           raise NotImplementedError,"no inner product for p and Bv implemented."           raise NotImplementedError("no inner product for p and Bv implemented.")
1630    
1631        def inner_p(self,p0,p1):        def inner_p(self,p0,p1):
1632           """           """
1633           Returns inner product of p0 and p1 (overwrite).           Returns inner product of p0 and p1 (overwrite).
1634    
1635           @param p0: a pressure           :param p0: a pressure
1636           @param p1: a pressure           :param p1: a pressure
1637           @return: inner product of p0 and p1           :return: inner product of p0 and p1
1638           @rtype: C{float}           :rtype: ``float``
1639           """           """
1640           raise NotImplementedError,"no inner product for p implemented."           raise NotImplementedError("no inner product for p implemented.")
1641        
1642        def norm_v(self,v):        def norm_v(self,v):
1643           """           """
1644           Returns the norm of v (overwrite).           Returns the norm of v (overwrite).
1645    
1646           @param v: a velovity           :param v: a velovity
1647           @return: norm of v           :return: norm of v
1648           @rtype: non-negative C{float}           :rtype: non-negative ``float``
1649           """           """
1650           raise NotImplementedError,"no norm of v implemented."           raise NotImplementedError("no norm of v implemented.")
1651        def getV(self, p, v0):        def getDV(self, p, v, tol):
1652           """           """
1653           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)
1654    
1655           @param p: a pressure           :param p: pressure
1656           @param v0: a initial guess for the value v to return.           :param v: pressure
1657           @return: v given as M{v= A^{-1} (f-B^*p)}           :return: dv given as *dv= A^{-1} (f-A v-B^*p)*
1658             :note: Only *A* may depend on *v* and *p*
1659           """           """
1660           raise NotImplementedError,"no v calculation implemented."           raise NotImplementedError("no dv calculation implemented.")
1661    
1662                    
1663        def Bv(self,v):        def Bv(self,v, tol):
1664          """          """
1665          Returns Bv (overwrite).          Returns Bv with accuracy `tol` (overwrite)
1666    
1667          @rtype: equal to the type of p          :rtype: equal to the type of p
1668          @note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1669          """          """
1670          raise NotImplementedError, "no operator B implemented."          raise NotImplementedError("no operator B implemented.")
1671    
1672        def norm_Bv(self,Bv):        def norm_Bv(self,Bv):
1673          """          """
1674          Returns the norm of Bv (overwrite).          Returns the norm of Bv (overwrite).
1675    
1676          @rtype: equal to the type of p          :rtype: equal to the type of p
1677          @note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1678          """          """
1679          raise NotImplementedError, "no norm of Bv implemented."          raise NotImplementedError("no norm of Bv implemented.")
1680    
1681        def solve_AinvBt(self,p):        def solve_AinvBt(self,dp, tol):
1682           """           """
1683           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}           Solves *A dv=B^*dp* with accuracy `tol`
          (overwrite).  
1684    
1685           @param p: a pressure increment           :param dp: a pressure increment
1686           @return: the solution of M{Av=B^*p}           :return: the solution of *A dv=B^*dp*
1687           @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.
1688           """           """
1689           raise NotImplementedError,"no operator A implemented."           raise NotImplementedError("no operator A implemented.")
1690    
1691        def solve_prec(self,Bv):        def solve_prec(self,Bv, tol):
1692           """           """
1693           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).  
1694    
1695           @rtype: equal to the type of p           :rtype: equal to the type of p
1696           @note: boundary conditions on p should be zero!           :note: boundary conditions on p should be zero!
          """  
          raise NotImplementedError,"no preconditioner for Schur complement implemented."  
       def setSubProblemTolerance(self):  
          """  
      Updates the tolerance for subproblems  
      @note: method is typically the method is overwritten.  
1697           """           """
1698           pass           raise NotImplementedError("no preconditioner for Schur complement implemented.")
1699        #=============================================================        #=============================================================
1700        def __Aprod_PCG(self,p):        def __Aprod_PCG(self,dp):
1701            dv=self.solve_AinvBt(p)            dv=self.solve_AinvBt(dp, self.__subtol)
1702            return ArithmeticTuple(dv,self.Bv(dv))            return ArithmeticTuple(dv,self.Bv(dv, self.__subtol))
1703    
1704        def __inner_PCG(self,p,r):        def __inner_PCG(self,p,r):
1705           return self.inner_pBv(p,r[1])           return self.inner_pBv(p,r[1])
1706    
1707        def __Msolve_PCG(self,r):        def __Msolve_PCG(self,r):
1708            return self.solve_prec(r[1])            return self.solve_prec(r[1], self.__subtol)
1709        #=============================================================        #=============================================================
1710        def __Aprod_GMRES(self,p):        def __Aprod_GMRES(self,p):
1711            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)
1712        def __inner_GMRES(self,p0,p1):        def __inner_GMRES(self,p0,p1):
1713           return self.inner_p(p0,p1)           return self.inner_p(p0,p1)
1714    
1715        #=============================================================        #=============================================================
1716        def norm_p(self,p):        def norm_p(self,p):
1717            """            """
1718            calculates the norm of C{p}            calculates the norm of ``p``
1719                        
1720            @param p: a pressure            :param p: a pressure
1721            @return: the norm of C{p} using the inner product for pressure            :return: the norm of ``p`` using the inner product for pressure
1722            @rtype: C{float}            :rtype: ``float``
1723            """            """
1724            f=self.inner_p(p,p)            f=self.inner_p(p,p)
1725            if f<0: raise ValueError,"negative pressure norm."            if f<0: raise ValueError("negative pressure norm.")
1726            return math.sqrt(f)            return math.sqrt(f)
       def adaptSubTolerance(self):  
       """  
       Returns True if tolerance adaption for subproblem is choosen.  
       """  
           self.__adaptSubTolerance  
1727                
1728        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):
1729           """           """
1730           Solves the saddle point problem using initial guesses v and p.           Solves the saddle point problem using initial guesses v and p.
1731    
1732           @param v: initial guess for velocity           :param v: initial guess for velocity
1733           @param p: initial guess for pressure           :param p: initial guess for pressure
1734           @type v: L{Data}           :type v: `Data`
1735           @type p: L{Data}           :type p: `Data`
1736           @param usePCG: indicates the usage of the PCG rather than GMRES scheme.           :param usePCG: indicates the usage of the PCG rather than GMRES scheme.
1737           @param max_iter: maximum number of iteration steps per correction           :param max_iter: maximum number of iteration steps per correction
1738                            attempt                            attempt
1739           @param verbose: if True, shows information on the progress of the           :param verbose: if True, shows information on the progress of the
1740                           saddlepoint problem solver.                           saddlepoint problem solver.
1741           @param iter_restart: restart the iteration after C{iter_restart} steps           :param iter_restart: restart the iteration after ``iter_restart`` steps
1742                                (only used if useUzaw=False)                                (only used if useUzaw=False)
1743           @type usePCG: C{bool}           :type usePCG: ``bool``
1744           @type max_iter: C{int}           :type max_iter: ``int``
1745           @type verbose: C{bool}           :type verbose: ``bool``
1746           @type iter_restart: C{int}           :type iter_restart: ``int``
1747           @rtype: C{tuple} of L{Data} objects           :rtype: ``tuple`` of `Data` objects
1748             :note: typically this method is overwritten by a subclass. It provides a wrapper for the ``_solve`` method.
1749             """
1750             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)
1751    
1752          def _solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1753             """
1754             see `_solve` method.
1755           """           """
1756           self.verbose=verbose           self.verbose=verbose
1757           rtol=self.getTolerance()           rtol=self.getTolerance()
1758           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
1759       if self.adaptSubTolerance(): self.setSubProblemTolerance()  
1760             K_p=self.__K_p
1761             K_v=self.__K_v
1762           correction_step=0           correction_step=0
1763           converged=False           converged=False
1764             chi=None
1765             eps=None
1766    
1767             if self.verbose: print(("HomogeneousSaddlePointProblem: start iteration: rtol= %e, atol=%e"%(rtol, atol)))
1768           while not converged:           while not converged:
               # calculate velocity for current pressure:  
               v=self.getV(p,v)  
               Bv=self.Bv(v)  
               norm_v=self.norm_v(v)  
               norm_Bv=self.norm_Bv(Bv)  
               ATOL=norm_v*rtol+atol  
               if self.verbose: print "HomogeneousSaddlePointProblem: norm v= %e, norm_Bv= %e, tolerance = %e."%(norm_v, norm_Bv,ATOL)  
               if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."  
               if norm_Bv <= ATOL:  
                  converged=True  
               else:  
                  correction_step+=1  
                  if correction_step>max_correction_steps:  
                       raise CorrectionFailed,"Given up after %d correction steps."%correction_step  
                  dp=self.solve_prec(Bv)  
                  if usePCG:  
                    norm2=self.inner_pBv(dp,Bv)  
                    if norm2<0: raise ValueError,"negative PCG norm."  
                    norm2=math.sqrt(norm2)  
                  else:  
                    norm2=self.norm_p(dp)  
                  ATOL_ITER=ATOL/norm_Bv*norm2*0.5  
                  if self.verbose: print "HomogeneousSaddlePointProblem: tolerance for solver: %e"%ATOL_ITER  
                  if usePCG:  
                        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)  
                  else:  
                        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: tolerance reached."  
      return v,p  
1769    
1770                 # get tolerance for velecity increment:
1771                 if chi == None:
1772                    rtol_v=self.__rtol_max
1773                 else:
1774                    rtol_v=min(chi/K_v,self.__rtol_max)
1775                 rtol_v=max(rtol_v, self.__rtol_min)
1776                 if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: rtol_v= %e"%(correction_step,rtol_v)))
1777                 # get velocity increment:
1778                 dv1=self.getDV(p,v,rtol_v)
1779                 v1=v+dv1
1780                 Bv1=self.Bv(v1, rtol_v)
1781                 norm_Bv1=self.norm_Bv(Bv1)
1782                 norm_dv1=self.norm_v(dv1)
1783                 if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: norm_Bv1 = %e, norm_dv1 = %e"%(correction_step, norm_Bv1, norm_dv1)))
1784                 if norm_dv1*self.__theta < norm_Bv1:
1785                    # get tolerance for pressure increment:
1786                    large_Bv1=True
1787                    if chi == None or eps == None:
1788                       rtol_p=self.__rtol_max
1789                    else:
1790                       rtol_p=min(chi**2*eps/K_p/norm_Bv1, self.__rtol_max)
1791                    self.__subtol=max(rtol_p**2, self.__rtol_min)
1792                    if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: rtol_p= %e"%(correction_step,rtol_p)))
1793                    # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1
1794                    if usePCG:
1795                        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)
1796                        v2=r[0]
1797                        Bv2=r[1]
1798                    else:
1799                        # don't use!!!!
1800                        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)
1801                        dv2=self.solve_AinvBt(dp, self.__subtol)
1802                        v2=v1-dv2
1803                        Bv2=self.Bv(v2, self.__subtol)
1804                    p2=p+dp
1805                 else:
1806                    large_Bv1=False
1807                    v2=v1
1808                    p2=p
1809                 # update business:
1810                 norm_dv2=self.norm_v(v2-v)
1811                 norm_v2=self.norm_v(v2)
1812                 if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: v2 = %e, norm_dv2 = %e"%(correction_step, norm_v2, self.norm_v(v2-v))))
1813                 eps, eps_old = max(norm_Bv1, norm_dv2), eps
1814                 if eps_old == None:
1815                      chi, chi_old = None, chi
1816                 else:
1817                      chi, chi_old = min(eps/ eps_old, self.__chi_max), chi
1818                 if eps != None:
1819                     if chi !=None:
1820                        if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: convergence rate = %e, correction = %e"%(correction_step,chi, eps)))
1821                     else:
1822                        if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: correction = %e"%(correction_step, eps)))
1823                 if eps <= rtol*norm_v2+atol :
1824                     converged = True
1825                 else:
1826                     if correction_step>=max_correction_steps:
1827                          raise CorrectionFailed("Given up after %d correction steps."%correction_step)
1828                     if chi_old!=None:
1829                        K_p=max(1,self.__reduction_factor*K_p,(chi-chi_old)/chi_old**2*K_p)
1830                        K_v=max(1,self.__reduction_factor*K_v,(chi-chi_old)/chi_old**2*K_p)
1831                        if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: new adjustment factor K = %e"%(correction_step,K_p)))
1832                 correction_step+=1
1833                 v,p =v2, p2
1834             if self.verbose: print(("HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step))
1835             return v,p
1836        #========================================================================        #========================================================================
1837        def setTolerance(self,tolerance=1.e-4):        def setTolerance(self,tolerance=1.e-4):
1838           """           """
1839           Sets the relative tolerance for (v,p).           Sets the relative tolerance for (v,p).
1840    
1841           @param tolerance: tolerance to be used           :param tolerance: tolerance to be used
1842           @type tolerance: non-negative C{float}           :type tolerance: non-negative ``float``
1843           """           """
1844           if tolerance<0:           if tolerance<0:
1845               raise ValueError,"tolerance must be positive."               raise ValueError("tolerance must be positive.")
1846           self.__rtol=tolerance           self.__rtol=tolerance
1847    
1848        def getTolerance(self):        def getTolerance(self):
1849           """           """
1850           Returns the relative tolerance.           Returns the relative tolerance.
1851    
1852           @return: relative tolerance           :return: relative tolerance
1853           @rtype: C{float}           :rtype: ``float``
1854           """           """
1855           return self.__rtol           return self.__rtol
1856    
# Line 1698  class HomogeneousSaddlePointProblem(obje Line 1858  class HomogeneousSaddlePointProblem(obje
1858           """           """
1859           Sets the absolute tolerance.           Sets the absolute tolerance.
1860    
1861           @param tolerance: tolerance to be used           :param tolerance: tolerance to be used
1862           @type tolerance: non-negative C{float}           :type tolerance: non-negative ``float``
1863           """           """
1864           if tolerance<0:           if tolerance<0:
1865               raise ValueError,"tolerance must be non-negative."               raise ValueError("tolerance must be non-negative.")
1866           self.__atol=tolerance           self.__atol=tolerance
1867    
1868        def getAbsoluteTolerance(self):        def getAbsoluteTolerance(self):
1869           """           """
1870           Returns the absolute tolerance.           Returns the absolute tolerance.
1871    
1872           @return: absolute tolerance           :return: absolute tolerance
1873           @rtype: C{float}           :rtype: ``float``
1874           """           """
1875           return self.__atol           return self.__atol
1876    
       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)  
1877    
1878  def MaskFromBoundaryTag(domain,*tags):  def MaskFromBoundaryTag(domain,*tags):
1879     """     """
# Line 1730  def MaskFromBoundaryTag(domain,*tags): Line 1882  def MaskFromBoundaryTag(domain,*tags):
1882    
1883     Usage: m=MaskFromBoundaryTag(domain, "left", "right")     Usage: m=MaskFromBoundaryTag(domain, "left", "right")
1884    
1885     @param domain: domain to be used     :param domain: domain to be used
1886     @type domain: L{escript.Domain}     :type domain: `escript.Domain`
1887     @param tags: boundary tags     :param tags: boundary tags
1888     @type tags: C{str}     :type tags: ``str``
1889     @return: a mask which marks samples that are touching the boundary tagged     :return: a mask which marks samples that are touching the boundary tagged
1890              by any of the given tags              by any of the given tags
1891     @rtype: L{escript.Data} of rank 0     :rtype: `escript.Data` of rank 0
1892     """     """
1893     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1894     d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))     d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
# Line 1751  def MaskFromTag(domain,*tags): Line 1903  def MaskFromTag(domain,*tags):
1903    
1904     Usage: m=MaskFromTag(domain, "ham")     Usage: m=MaskFromTag(domain, "ham")
1905    
1906     @param domain: domain to be used     :param domain: domain to be used
1907     @type domain: L{escript.Domain}     :type domain: `escript.Domain`
1908     @param tags: boundary tags     :param tags: boundary tags
1909     @type tags: C{str}     :type tags: ``str``
1910     @return: a mask which marks samples that are touching the boundary tagged     :return: a mask which marks samples that are touching the boundary tagged
1911              by any of the given tags              by any of the given tags
1912     @rtype: L{escript.Data} of rank 0     :rtype: `escript.Data` of rank 0
1913     """     """
1914     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1915     d=escript.Scalar(0.,escript.Function(domain))     d=escript.Scalar(0.,escript.Function(domain))

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

  ViewVC Help
Powered by ViewVC 1.1.26