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

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

  ViewVC Help
Powered by ViewVC 1.1.26