/[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 1956 by gross, Mon Nov 3 05:08:42 2008 UTC revision 2647 by gross, Fri Sep 4 05:25:25 2009 UTC
# Line 1  Line 1 
1    
2  ########################################################  ########################################################
3  #  #
4  # Copyright (c) 2003-2008 by University of Queensland  # Copyright (c) 2003-2009 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-2009 by University of Queensland
15  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
16  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
17  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
18  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
19  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
20  __url__="http://www.uq.edu.au/esscc/escript-finley"  __url__="https://launchpad.net/escript-finley"
21    
22  """  """
23  Provides some tools related to PDEs.  Provides some tools related to PDEs.
24    
25  Currently includes:  Currently includes:
26      - Projector - to project a discontinuous      - Projector - to project a discontinuous function onto a continuous function
27      - Locator - to trace values in data objects at a certain location      - Locator - to trace values in data objects at a certain location
28      - TimeIntegrationManager - to handel extraplotion in time      - TimeIntegrationManager - to handle extrapolation in time
29      - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme      - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme
30    
31  @var __author__: name of author  :var __author__: name of author
32  @var __copyright__: copyrights  :var __copyright__: copyrights
33  @var __license__: licence agreement  :var __license__: licence agreement
34  @var __url__: url entry point on documentation  :var __url__: url entry point on documentation
35  @var __version__: version  :var __version__: version
36  @var __date__: date of the version  :var __date__: date of the version
37  """  """
38    
39  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
# Line 41  __author__="Lutz Gross, l.gross@uq.edu.a Line 41  __author__="Lutz Gross, l.gross@uq.edu.a
41    
42  import escript  import escript
43  import linearPDEs  import linearPDEs
44  import numarray  import numpy
45  import util  import util
46  import math  import math
47    
 ##### Added by Artak  
 # from Numeric import zeros,Int,Float64  
 ###################################  
   
   
48  class TimeIntegrationManager:  class TimeIntegrationManager:
49    """    """
50    a simple mechanism to manage time dependend values.    A simple mechanism to manage time dependend values.
51    
52    typical usage is::    Typical usage is::
53    
54       dt=0.1 # time increment       dt=0.1 # time increment
55       tm=TimeIntegrationManager(inital_value,p=1)       tm=TimeIntegrationManager(inital_value,p=1)
# 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 inital_value is the initial value and p is order used for extrapolation       Sets up the value manager where ``inital_values`` are the initial values
67         and p is the order used for extrapolation.
68       """       """
69       if kwargs.has_key("p"):       if kwargs.has_key("p"):
70              self.__p=kwargs["p"]              self.__p=kwargs["p"]
# Line 85  class TimeIntegrationManager: Line 81  class TimeIntegrationManager:
81    
82    def getTime(self):    def getTime(self):
83        return self.__t        return self.__t
84    
85    def getValue(self):    def getValue(self):
86        out=self.__v_mem[0]        out=self.__v_mem[0]
87        if len(out)==1:        if len(out)==1:
# Line 94  class TimeIntegrationManager: Line 91  class TimeIntegrationManager:
91    
92    def checkin(self,dt,*values):    def checkin(self,dt,*values):
93        """        """
94        adds new values to the manager. the p+1 last value get lost        Adds new values to the manager. The p+1 last values are lost.
95        """        """
96        o=min(self.__order+1,self.__p)        o=min(self.__order+1,self.__p)
97        self.__order=min(self.__order+1,self.__p)        self.__order=min(self.__order+1,self.__p)
# Line 111  class TimeIntegrationManager: Line 108  class TimeIntegrationManager:
108    
109    def extrapolate(self,dt):    def extrapolate(self,dt):
110        """        """
111        extrapolates to 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 126  class TimeIntegrationManager: Line 123  class TimeIntegrationManager:
123           return out[0]           return out[0]
124        else:        else:
125           return out           return out
126    
127    
128  class Projector:  class Projector:
129    """    """
130    The Projector is a factory which projects a discontiuous function onto a    The Projector is a factory which projects a discontinuous function onto a
131    continuous function on the a given domain.    continuous function on a given domain.
132    """    """
133    def __init__(self, domain, reduce = True, fast=True):    def __init__(self, domain, reduce=True, fast=True):
134      """      """
135      Create 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 (default is True)      :param reduce: Flag to reduce projection order
139      @param fast: Flag to use a fast method based on matrix lumping (default is true)      :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:
143        self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)          self.__pde.getSolverOptions().setSolverMethod(linearPDEs.SolverOptions.LUMPING)
144      self.__pde.setSymmetryOn()      self.__pde.setSymmetryOn()
145      self.__pde.setReducedOrderTo(reduce)      self.__pde.setReducedOrderTo(reduce)
146      self.__pde.setValue(D = 1.)      self.__pde.setValue(D = 1.)
147      return      return
148      def getSolverOptions(self):
149        """
150        Returns the solver options of the PDE solver.
151        
152        :rtype: `linearPDEs.SolverOptions`
153        """
154    
155    def __call__(self, input_data):    def __call__(self, input_data):
156      """      """
157      Projects input_data onto a continuous function      Projects ``input_data`` onto a continuous function.
158    
159      @param input_data: The input_data to be projected.      :param input_data: the data to be projected
160      """      """
161      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
162      self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())      self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())
# Line 186  class Projector: Line 189  class Projector:
189    
190  class NoPDE:  class NoPDE:
191       """       """
192       solves the following problem for u:       Solves the following problem for u:
193    
194       M{kronecker[i,j]*D[j]*u[j]=Y[i]}       *kronecker[i,j]*D[j]*u[j]=Y[i]*
195    
196       with constraint       with constraint
197    
198       M{u[j]=r[j]}  where M{q[j]>0}       *u[j]=r[j]*  where *q[j]>0*
199    
200       where D, Y, r and q are given functions of rank 1.       where *D*, *Y*, *r* and *q* are given functions of rank 1.
201    
202       In the case of scalars this takes the form       In the case of scalars this takes the form
203    
204       M{D*u=Y}       *D*u=Y*
205    
206       with constraint       with constraint
207    
208       M{u=r}  where M{q>0}       *u=r* where *q>0*
209    
210       where D, Y, r and q are given scalar functions.       where *D*, *Y*, *r* and *q* are given scalar functions.
211    
212       The constraint is overwriting any other condition.       The constraint overwrites any other condition.
213    
214       @note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention       :note: This class is similar to the `linearPDEs.LinearPDE` class with
215              that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole              A=B=C=X=0 but has the intention that all input parameters are given
216              thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.              in `Solution` or `ReducedSolution`.
217       """       """
218         # The whole thing is a bit strange and I blame Rob Woodcock (CSIRO) for
219         # this.
220       def __init__(self,domain,D=None,Y=None,q=None,r=None):       def __init__(self,domain,D=None,Y=None,q=None,r=None):
221           """           """
222           initialize the problem           Initializes the problem.
223    
224           @param domain: domain of the PDE.           :param domain: domain of the PDE
225           @type domain: L{Domain}           :type domain: `Domain`
226           @param D: coefficient of the solution.           :param D: coefficient of the solution
227           @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type D: ``float``, ``int``, ``numpy.ndarray``, `Data`
228           @param Y: right hand side           :param Y: right hand side
229           @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type Y: ``float``, ``int``, ``numpy.ndarray``, `Data`
230           @param q: location of constraints           :param q: location of constraints
231           @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type q: ``float``, ``int``, ``numpy.ndarray``, `Data`
232           @param r: value of solution at locations of constraints           :param r: value of solution at locations of constraints
233           @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type r: ``float``, ``int``, ``numpy.ndarray``, `Data`
234           """           """
235           self.__domain=domain           self.__domain=domain
236           self.__D=D           self.__D=D
# Line 234  class NoPDE: Line 239  class NoPDE:
239           self.__r=r           self.__r=r
240           self.__u=None           self.__u=None
241           self.__function_space=escript.Solution(self.__domain)           self.__function_space=escript.Solution(self.__domain)
242    
243       def setReducedOn(self):       def setReducedOn(self):
244           """           """
245           sets the L{FunctionSpace} of the solution to L{ReducedSolution}           Sets the `FunctionSpace` of the solution to `ReducedSolution`.
246           """           """
247           self.__function_space=escript.ReducedSolution(self.__domain)           self.__function_space=escript.ReducedSolution(self.__domain)
248           self.__u=None           self.__u=None
249    
250       def setReducedOff(self):       def setReducedOff(self):
251           """           """
252           sets the L{FunctionSpace} of the solution to L{Solution}           Sets the `FunctionSpace` of the solution to `Solution`.
253           """           """
254           self.__function_space=escript.Solution(self.__domain)           self.__function_space=escript.Solution(self.__domain)
255           self.__u=None           self.__u=None
256            
257       def setValue(self,D=None,Y=None,q=None,r=None):       def setValue(self,D=None,Y=None,q=None,r=None):
258           """           """
259           assigns values to the parameters.           Assigns values to the parameters.
260    
261           @param D: coefficient of the solution.           :param D: coefficient of the solution
262           @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type D: ``float``, ``int``, ``numpy.ndarray``, `Data`
263           @param Y: right hand side           :param Y: right hand side
264           @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type Y: ``float``, ``int``, ``numpy.ndarray``, `Data`
265           @param q: location of constraints           :param q: location of constraints
266           @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type q: ``float``, ``int``, ``numpy.ndarray``, `Data`
267           @param r: value of solution at locations of constraints           :param r: value of solution at locations of constraints
268           @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type r: ``float``, ``int``, ``numpy.ndarray``, `Data`
269           """           """
270           if not D==None:           if not D==None:
271              self.__D=D              self.__D=D
# Line 276  class NoPDE: Line 282  class NoPDE:
282    
283       def getSolution(self):       def getSolution(self):
284           """           """
285           returns the solution           Returns the solution.
286            
287           @return: the solution of the problem           :return: the solution of the problem
288           @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.           :rtype: `Data` object in the `FunctionSpace` `Solution` or
289                     `ReducedSolution`
290           """           """
291           if self.__u==None:           if self.__u==None:
292              if self.__D==None:              if self.__D==None:
# Line 296  class NoPDE: Line 303  class NoPDE:
303                  self.__u*=(1.-q)                  self.__u*=(1.-q)
304                  if not self.__r==None: self.__u+=q*self.__r                  if not self.__r==None: self.__u+=q*self.__r
305           return self.__u           return self.__u
306                
307  class Locator:  class Locator:
308       """       """
309       Locator provides access to the values of data objects at a given       Locator provides access to the values of data objects at a given spatial
310       spatial coordinate x.         coordinate x.
311        
312       In fact, a Locator object finds the sample in the set of samples of a       In fact, a Locator object finds the sample in the set of samples of a
313       given function space or domain where which is closest to the given       given function space or domain which is closest to the given point x.
      point x.  
314       """       """
315    
316       def __init__(self,where,x=numarray.zeros((3,))):       def __init__(self,where,x=numpy.zeros((3,))):
317         """         """
318         Initializes a Locator to access values in Data objects on the Doamin         Initializes a Locator to access values in Data objects on the Doamin
319         or FunctionSpace where for the sample point which         or FunctionSpace for the sample point which is closest to the given
320         closest to the given point x.         point x.
321    
322         @param where: function space         :param where: function space
323         @type where: L{escript.FunctionSpace}         :type where: `escript.FunctionSpace`
324         @param x: coefficient of the solution.         :param x: location(s) of the Locator
325         @type x: L{numarray.NumArray} or C{list} of L{numarray.NumArray}         :type x: ``numpy.ndarray`` or ``list`` of ``numpy.ndarray``
326         """         """
327         if isinstance(where,escript.FunctionSpace):         if isinstance(where,escript.FunctionSpace):
328            self.__function_space=where            self.__function_space=where
329         else:         else:
330            self.__function_space=escript.ContinuousFunction(where)            self.__function_space=escript.ContinuousFunction(where)
331           iterative=False
332         if isinstance(x, list):         if isinstance(x, list):
333               if len(x)==0:
334                  raise "ValueError", "At least one point must be given."
335               try:
336                 iter(x[0])
337                 iterative=True
338               except TypeError:
339                 iterative=False
340           if iterative:
341             self.__id=[]             self.__id=[]
342             for p in x:             for p in x:
343                self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())                self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())
# Line 350  class Locator: Line 365  class Locator:
365    
366       def getX(self):       def getX(self):
367          """          """
368      Returns the exact coordinates of the Locator.          Returns the exact coordinates of the Locator.
369      """          """
370          return self(self.getFunctionSpace().getX())          return self(self.getFunctionSpace().getX())
371    
372       def getFunctionSpace(self):       def getFunctionSpace(self):
373          """          """
374      Returns the function space of the Locator.          Returns the function space of the Locator.
375      """          """
376          return self.__function_space          return self.__function_space
377    
378       def getId(self,item=None):       def getId(self,item=None):
379          """          """
380      Returns the identifier of the location.          Returns the identifier of the location.
381      """          """
382          if item == None:          if item == None:
383             return self.__id             return self.__id
384          else:          else:
# Line 375  class Locator: Line 390  class Locator:
390    
391       def __call__(self,data):       def __call__(self,data):
392          """          """
393      Returns the value of data at the Locator of a Data object otherwise          Returns the value of data at the Locator of a Data object.
394      the object is returned.          """
     """  
395          return self.getValue(data)          return self.getValue(data)
396    
397       def getValue(self,data):       def getValue(self,data):
398          """          """
399      Returns the value of data at the Locator if data is a Data object          Returns the value of ``data`` at the Locator if ``data`` is a `Data`
400      otherwise the object is returned.          object otherwise the object is returned.
401      """          """
402          if isinstance(data,escript.Data):          if isinstance(data,escript.Data):
403             if data.getFunctionSpace()==self.getFunctionSpace():             dat=util.interpolate(data,self.getFunctionSpace())
              dat=data  
            else:  
              dat=data.interpolate(self.getFunctionSpace())  
404             id=self.getId()             id=self.getId()
405             r=data.getRank()             r=data.getRank()
406             if isinstance(id,list):             if isinstance(id,list):
407                 out=[]                 out=[]
408                 for i in id:                 for i in id:
409                    o=data.getValueOfGlobalDataPoint(*i)                    o=numpy.array(dat.getTupleForGlobalDataPoint(*i))
410                    if data.getRank()==0:                    if data.getRank()==0:
411                       out.append(o[0])                       out.append(o[0])
412                    else:                    else:
413                       out.append(o)                       out.append(o)
414                 return out                 return out
415             else:             else:
416               out=data.getValueOfGlobalDataPoint(*id)               out=numpy.array(dat.getTupleForGlobalDataPoint(*id))
417               if data.getRank()==0:               if data.getRank()==0:
418                  return out[0]                  return out[0]
419               else:               else:
# Line 412  class Locator: Line 423  class Locator:
423    
424  class SolverSchemeException(Exception):  class SolverSchemeException(Exception):
425     """     """
426     exceptions thrown by solvers     This is a generic exception thrown by solvers.
427     """     """
428     pass     pass
429    
430  class IndefinitePreconditioner(SolverSchemeException):  class IndefinitePreconditioner(SolverSchemeException):
431     """     """
432     the preconditioner is not positive definite.     Exception thrown if the preconditioner is not positive definite.
433     """     """
434     pass     pass
435    
436  class MaxIterReached(SolverSchemeException):  class MaxIterReached(SolverSchemeException):
437     """     """
438     maxium number of iteration steps is reached.     Exception thrown if the maximum number of iteration steps is reached.
439     """     """
440     pass     pass
441  class IterationBreakDown(SolverSchemeException):  
442    class CorrectionFailed(SolverSchemeException):
443     """     """
444     iteration scheme econouters an incurable breakdown.     Exception thrown if no convergence has been achieved in the solution
445       correction scheme.
446     """     """
447     pass     pass
448  class NegativeNorm(SolverSchemeException):  
449    class IterationBreakDown(SolverSchemeException):
450     """     """
451     a norm calculation returns a negative norm.     Exception thrown if the iteration scheme encountered an incurable breakdown.
452     """     """
453     pass     pass
454    
455  class IterationHistory(object):  class NegativeNorm(SolverSchemeException):
456     """     """
457     The IterationHistory class is used to define a stopping criterium. It keeps track of the     Exception thrown if a norm calculation returns a negative norm.
    residual norms. The stoppingcriterium indicates termination if the residual norm has been reduced by  
    a given tolerance.  
458     """     """
459     def __init__(self,tolerance=math.sqrt(util.EPSILON),verbose=False):     pass
       """  
       Initialization  
   
       @param tolerance: tolerance  
       @type tolerance: positive C{float}  
       @param verbose: switches on the printing out some information  
       @type verbose: C{bool}  
       """  
       if not tolerance>0.:  
           raise ValueError,"tolerance needs to be positive."  
       self.tolerance=tolerance  
       self.verbose=verbose  
       self.history=[]  
    def stoppingcriterium(self,norm_r,r,x):  
        """  
        returns True if the C{norm_r} is C{tolerance}*C{norm_r[0]} where C{norm_r[0]}  is the residual norm at the first call.  
   
         
        @param norm_r: current residual norm  
        @type norm_r: non-negative C{float}  
        @param r: current residual (not used)  
        @param x: current solution approximation (not used)  
        @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.  
        @rtype: C{bool}  
   
        """  
        self.history.append(norm_r)  
        if self.verbose: print "iter: %s:  inner(rhat,r) = %e"#(len(self.history)-1, self.history[-1])  
        return self.history[-1]<=self.tolerance * self.history[0]  
460    
461     def stoppingcriterium2(self,norm_r,norm_b,solver="GMRES",TOL=None):  def PCG(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100, initial_guess=True, verbose=False):
462         """     """
463         returns True if the C{norm_r} is C{tolerance}*C{norm_b}     Solver for
464    
465             *Ax=b*
        @param norm_r: current residual norm  
        @type norm_r: non-negative C{float}  
        @param norm_b: norm of right hand side  
        @type norm_b: non-negative C{float}  
        @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.  
        @rtype: C{bool}  
466    
467         """     with a symmetric and positive definite operator A (more details required!).
468         if TOL==None:     It uses the conjugate gradient method with preconditioner M providing an
469            TOL=self.tolerance     approximation of A.
        self.history.append(norm_r)  
        if self.verbose: print "iter: %s:  norm(r) = %e"#(len(self.history)-1, self.history[-1])  
        return self.history[-1]<=TOL * norm_b  
470    
471  def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):     The iteration is terminated if
    """  
    Solver for  
472    
473     M{Ax=b}     *|r| <= atol+rtol*|r0|*
474    
475     with a symmetric and positive definite operator A (more details required!).     where *r0* is the initial residual and *|.|* is the energy norm. In fact
    It uses the conjugate gradient method with preconditioner M providing an approximation of A.  
476    
477     The iteration is terminated if the C{stoppingcriterium} function return C{True}.     *|r| = sqrt( bilinearform(Msolve(r),r))*
478    
479     For details on the preconditioned conjugate gradient method see the book:     For details on the preconditioned conjugate gradient method see the book:
480    
481     Templates for the Solution of Linear Systems by R. Barrett, M. Berry,     I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
482     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,
483     C. Romine, and H. van der Vorst.     C. Romine, and H. van der Vorst}.
484    
485     @param b: the right hand side of the liner system. C{b} is altered.     :param r: initial residual *r=b-Ax*. ``r`` is altered.
486     @type b: 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)
487     @param Aprod: returns the value Ax     :param x: an initial guess for the solution
488     @type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}. The returned object needs to be of the same type like argument C{b}.     :type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
489     @param Msolve: solves Mx=r     :param Aprod: returns the value Ax
490     @type Msolve: function C{Msolve(r)} where C{r} is of the same type like argument C{b}. The returned object needs to be of the same     :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
491  type like argument C{x}.                  argument ``x``. The returned object needs to be of the same type
492     @param bilinearform: inner product C{<x,r>}                  like argument ``r``.
493     @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is . The returned value is a C{float}.     :param Msolve: solves Mx=r
494     @param stoppingcriterium: function which returns True if a stopping criterium is meet. C{stoppingcriterium} has the arguments C{norm_r}, C{r} and C{x} giving the current norm of the residual (=C{sqrt(bilinearform(Msolve(r),r)}), the current residual and the current solution approximation. C{stoppingcriterium} is called in each iteration step.     :type Msolve: function ``Msolve(r)`` where ``r`` is of the same type like
495     @type stoppingcriterium: function that returns C{True} or C{False}                   argument ``r``. The returned object needs to be of the same
496     @param x: an initial guess for the solution. If no C{x} is given 0*b is used.                   type like argument ``x``.
497     @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)     :param bilinearform: inner product ``<x,r>``
498     @param iter_max: maximum number of iteration steps.     :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
499     @type iter_max: C{int}                         type like argument ``x`` and ``r`` is. The returned value
500     @return: the solution approximation and the corresponding residual                         is a ``float``.
501     @rtype: C{tuple}     :param atol: absolute tolerance
502     @warning: C{b} and C{x} are altered.     :type atol: non-negative ``float``
503       :param rtol: relative tolerance
504       :type rtol: non-negative ``float``
505       :param iter_max: maximum number of iteration steps
506       :type iter_max: ``int``
507       :return: the solution approximation and the corresponding residual
508       :rtype: ``tuple``
509       :warning: ``r`` and ``x`` are altered.
510     """     """
511     iter=0     iter=0
    if x==None:  
       x=0*b  
    else:  
       b += (-1)*Aprod(x)  
    r=b  
512     rhat=Msolve(r)     rhat=Msolve(r)
513     d = rhat     d = rhat
514     rhat_dot_r = bilinearform(rhat, r)     rhat_dot_r = bilinearform(rhat, r)
515     if rhat_dot_r<0: raise NegativeNorm,"negative norm."     if rhat_dot_r<0: raise NegativeNorm,"negative norm."
516       norm_r0=math.sqrt(rhat_dot_r)
517       atol2=atol+rtol*norm_r0
518       if atol2<=0:
519          raise ValueError,"Non-positive tolarance."
520       atol2=max(atol2, 100. * util.EPSILON * norm_r0)
521    
522     while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):     if verbose: print "PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_r0, atol2)
523    
524    
525       while not math.sqrt(rhat_dot_r) <= atol2:
526         iter+=1         iter+=1
527         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
528    
529         q=Aprod(d)         q=Aprod(d)
530         alpha = rhat_dot_r / bilinearform(d, q)         alpha = rhat_dot_r / bilinearform(d, q)
531         x += alpha * d         x += alpha * d
532         r += (-alpha) * q         if isinstance(q,ArithmeticTuple):
533           r += q * (-alpha)      # Doing it the other way calls the float64.__mul__ not AT.__rmul__
534           else:
535               r += (-alpha) * q
536         rhat=Msolve(r)         rhat=Msolve(r)
537         rhat_dot_r_new = bilinearform(rhat, r)         rhat_dot_r_new = bilinearform(rhat, r)
538         beta = rhat_dot_r_new / rhat_dot_r         beta = rhat_dot_r_new / rhat_dot_r
# Line 557  type like argument C{x}. Line 541  type like argument C{x}.
541    
542         rhat_dot_r = rhat_dot_r_new         rhat_dot_r = rhat_dot_r_new
543         if rhat_dot_r<0: raise NegativeNorm,"negative norm."         if rhat_dot_r<0: raise NegativeNorm,"negative norm."
544           if verbose: print "PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))
545     return x,r     if verbose: print "PCG: tolerance reached after %s steps."%iter
546       return x,r,math.sqrt(rhat_dot_r)
547    
548  class Defect(object):  class Defect(object):
549      """      """
550      defines a non-linear defect F(x) of a variable x      Defines a non-linear defect F(x) of a variable x.
551      """      """
552      def __init__(self):      def __init__(self):
553          """          """
554          initialize defect          Initializes defect.
555          """          """
556          self.setDerivativeIncrementLength()          self.setDerivativeIncrementLength()
557    
558      def bilinearform(self, x0, x1):      def bilinearform(self, x0, x1):
559          """          """
560          returns the inner product of x0 and x1          Returns the inner product of x0 and x1
561          @param x0: a value for x  
562          @param x1: a value for x          :param x0: value for x0
563          @return: the inner product of x0 and x1          :param x1: value for x1
564          @rtype: C{float}          :return: the inner product of x0 and x1
565            :rtype: ``float``
566          """          """
567          return 0          return 0
568          
569      def norm(self,x):      def norm(self,x):
570          """          """
571          the norm of argument C{x}          Returns the norm of argument ``x``.
572    
573          @param x: a value for x          :param x: a value
574          @return: norm of argument x          :return: norm of argument x
575          @rtype: C{float}          :rtype: ``float``
576          @note: by default C{sqrt(self.bilinearform(x,x)} is retrurned.          :note: by default ``sqrt(self.bilinearform(x,x)`` is returned.
577          """          """
578          s=self.bilinearform(x,x)          s=self.bilinearform(x,x)
579          if s<0: raise NegativeNorm,"negative norm."          if s<0: raise NegativeNorm,"negative norm."
580          return math.sqrt(s)          return math.sqrt(s)
581    
   
582      def eval(self,x):      def eval(self,x):
583          """          """
584          returns the value F of a given x          Returns the value F of a given ``x``.
585    
586          @param x: value for which the defect C{F} is evalulated.          :param x: value for which the defect ``F`` is evaluated
587          @return: value of the defect at C{x}          :return: value of the defect at ``x``
588          """          """
589          return 0          return 0
590    
# Line 608  class Defect(object): Line 593  class Defect(object):
593    
594      def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):      def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):
595          """          """
596          sets the relative length of the increment used to approximate the derivative of the defect          Sets the relative length of the increment used to approximate the
597          the increment is inc*norm(x)/norm(v)*v in the direction of v with x as a staring point.          derivative of the defect. The increment is inc*norm(x)/norm(v)*v in the
598            direction of v with x as a starting point.
599    
600          @param inc: relative increment length          :param inc: relative increment length
601          @type inc: positive C{float}          :type inc: positive ``float``
602          """          """
603          if inc<=0: raise ValueError,"positive increment required."          if inc<=0: raise ValueError,"positive increment required."
604          self.__inc=inc          self.__inc=inc
605    
606      def getDerivativeIncrementLength(self):      def getDerivativeIncrementLength(self):
607          """          """
608          returns the relative increment length used to approximate the derivative of the defect          Returns the relative increment length used to approximate the
609          @return: value of the defect at C{x}          derivative of the defect.
610          @rtype: positive C{float}          :return: value of the defect at ``x``
611            :rtype: positive ``float``
612          """          """
613          return self.__inc          return self.__inc
614    
615      def derivative(self, F0, x0, v, v_is_normalised=True):      def derivative(self, F0, x0, v, v_is_normalised=True):
616          """          """
617          returns the directional derivative at x0 in the direction of v          Returns the directional derivative at ``x0`` in the direction of ``v``.
618    
619          @param F0: value of this defect at x0          :param F0: value of this defect at x0
620          @param x0: value at which derivative is calculated.          :param x0: value at which derivative is calculated
621          @param v: direction          :param v: direction
622          @param v_is_normalised: is true to indicate that C{v} is nomalized (self.norm(v)=0)          :param v_is_normalised: True to indicate that ``v`` is nomalized
623          @return: derivative of this defect at x0 in the direction of C{v}                                  (self.norm(v)=0)
624          @note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is used but this method          :return: derivative of this defect at x0 in the direction of ``v``
625          maybe oepsnew verwritten to use exact evalution.          :note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is
626                   used but this method maybe overwritten to use exact evaluation.
627          """          """
628          normx=self.norm(x0)          normx=self.norm(x0)
629          if normx>0:          if normx>0:
# Line 651  class Defect(object): Line 639  class Defect(object):
639          F1=self.eval(x0 + epsnew * v)          F1=self.eval(x0 + epsnew * v)
640          return (F1-F0)/epsnew          return (F1-F0)/epsnew
641    
642  ######################################      ######################################
643  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, sub_tol_max=0.5, gamma=0.9, verbose=False):
644     """     """
645     solves a non-linear problem M{F(x)=0} for unknown M{x} using the stopping criterion:     Solves a non-linear problem *F(x)=0* for unknown *x* using the stopping
646       criterion:
647    
648     M{norm(F(x) <= atol + rtol * norm(F(x0)}     *norm(F(x) <= atol + rtol * norm(F(x0)*
649      
650     where M{x0} is the initial guess.     where *x0* is the initial guess.
651    
652     @param defect: object defining the the function M{F}, C{defect.norm} defines the M{norm} used in the stopping criterion.     :param defect: object defining the function *F*. ``defect.norm`` defines the
653     @type defect: L{Defect}                    *norm* used in the stopping criterion.
654     @param x: initial guess for the solution, C{x} is altered.     :type defect: `Defect`
655     @type x: any object type allowing basic operations such as  L{numarray.NumArray}, L{Data}     :param x: initial guess for the solution, ``x`` is altered.
656     @param iter_max: maximum number of iteration steps     :type x: any object type allowing basic operations such as
657     @type iter_max: positive C{int}              ``numpy.ndarray``, `Data`
658     @param sub_iter_max:     :param iter_max: maximum number of iteration steps
659     @type sub_iter_max:     :type iter_max: positive ``int``
660     @param atol: absolute tolerance for the solution     :param sub_iter_max: maximum number of inner iteration steps
661     @type atol: positive C{float}     :type sub_iter_max: positive ``int``
662     @param rtol: relative tolerance for the solution     :param atol: absolute tolerance for the solution
663     @type rtol: positive C{float}     :type atol: positive ``float``
664     @param gamma: tolerance safety factor for inner iteration     :param rtol: relative tolerance for the solution
665     @type gamma: positive C{float}, less than 1     :type rtol: positive ``float``
666     @param sub_tol_max: upper bound for inner tolerance.     :param gamma: tolerance safety factor for inner iteration
667     @type sub_tol_max: positive C{float}, less than 1     :type gamma: positive ``float``, less than 1
668     @return: an approximation of the solution with the desired accuracy     :param sub_tol_max: upper bound for inner tolerance
669     @rtype: same type as the initial guess.     :type sub_tol_max: positive ``float``, less than 1
670       :return: an approximation of the solution with the desired accuracy
671       :rtype: same type as the initial guess
672     """     """
673     lmaxit=iter_max     lmaxit=iter_max
674     if atol<0: raise ValueError,"atol needs to be non-negative."     if atol<0: raise ValueError,"atol needs to be non-negative."
# Line 697  def NewtonGMRES(defect, x, iter_max=100, Line 688  def NewtonGMRES(defect, x, iter_max=100,
688     # main iteration loop     # main iteration loop
689     #     #
690     while not fnrm<=stop_tol:     while not fnrm<=stop_tol:
691              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
692              #              #
693          #   adjust sub_tol_          #   adjust sub_tol_
694          #          #
695              if iter > 1:              if iter > 1:
696             rat=fnrm/fnrmo             rat=fnrm/fnrmo
# Line 712  def NewtonGMRES(defect, x, iter_max=100, Line 703  def NewtonGMRES(defect, x, iter_max=100,
703              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown
704              #     if iter_restart -1 is returned as sub_iter              #     if iter_restart -1 is returned as sub_iter
705              #     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
706              #              #
707              #                #
708              if verbose: print "             subiteration (GMRES) is called with relative tolerance %e."%sub_tol              if verbose: print "             subiteration (GMRES) is called with relative tolerance %e."%sub_tol
709              try:              try:
710                 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, sub_tol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)
# Line 734  def NewtonGMRES(defect, x, iter_max=100, Line 725  def NewtonGMRES(defect, x, iter_max=100,
725    
726  def __givapp(c,s,vin):  def __givapp(c,s,vin):
727      """      """
728      apply a sequence of Givens rotations (c,s) to the recuirsively to the vector vin      Applies a sequence of Givens rotations (c,s) recursively to the vector
729      @warning: C{vin} is altered.      ``vin``
730    
731        :warning: ``vin`` is altered.
732      """      """
733      vrot=vin      vrot=vin
734      if isinstance(c,float):      if isinstance(c,float):
735          vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]          vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
736      else:      else:
737          for i in range(len(c)):          for i in range(len(c)):
738              w1=c[i]*vrot[i]-s[i]*vrot[i+1]              w1=c[i]*vrot[i]-s[i]*vrot[i+1]
739          w2=s[i]*vrot[i]+c[i]*vrot[i+1]          w2=s[i]*vrot[i]+c[i]*vrot[i+1]
740              vrot[i:i+2]=w1,w2              vrot[i]=w1
741                vrot[i+1]=w2
742      return vrot      return vrot
743    
744  def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):  def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
745     h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)     h=numpy.zeros((iter_restart,iter_restart),numpy.float64)
746     c=numarray.zeros(iter_restart,numarray.Float64)     c=numpy.zeros(iter_restart,numpy.float64)
747     s=numarray.zeros(iter_restart,numarray.Float64)     s=numpy.zeros(iter_restart,numpy.float64)
748     g=numarray.zeros(iter_restart,numarray.Float64)     g=numpy.zeros(iter_restart,numpy.float64)
749     v=[]     v=[]
750    
751     rho=defect.norm(F0)     rho=defect.norm(F0)
752     if rho<=0.: return x0*0     if rho<=0.: return x0*0
753      
754     v.append(-F0/rho)     v.append(-F0/rho)
755     g[0]=rho     g[0]=rho
756     iter=0     iter=0
757     while rho > atol and iter<iter_restart-1:     while rho > atol and iter<iter_restart-1:
758            if iter  >= iter_max:
759      if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max              raise MaxIterReached,"maximum number of %s steps reached."%iter_max
760    
761          p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)          p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
762      v.append(p)          v.append(p)
763    
764      v_norm1=defect.norm(v[iter+1])          v_norm1=defect.norm(v[iter+1])
765    
766          # Modified Gram-Schmidt          # Modified Gram-Schmidt
767      for j in range(iter+1):          for j in range(iter+1):
768           h[j][iter]=defect.bilinearform(v[j],v[iter+1])                h[j,iter]=defect.bilinearform(v[j],v[iter+1])
769           v[iter+1]-=h[j][iter]*v[j]              v[iter+1]-=h[j,iter]*v[j]
770          
771      h[iter+1][iter]=defect.norm(v[iter+1])          h[iter+1,iter]=defect.norm(v[iter+1])
772      v_norm2=h[iter+1][iter]          v_norm2=h[iter+1,iter]
773    
774          # Reorthogonalize if needed          # Reorthogonalize if needed
775      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)
776          for j in range(iter+1):                for j in range(iter+1):
777             hr=defect.bilinearform(v[j],v[iter+1])                  hr=defect.bilinearform(v[j],v[iter+1])
778                 h[j][iter]=h[j][iter]+hr                  h[j,iter]=h[j,iter]+hr
779                 v[iter+1] -= hr*v[j]                  v[iter+1] -= hr*v[j]
780    
781          v_norm2=defect.norm(v[iter+1])              v_norm2=defect.norm(v[iter+1])
782          h[iter+1][iter]=v_norm2              h[iter+1,iter]=v_norm2
783          #   watch out for happy breakdown          #   watch out for happy breakdown
784          if not v_norm2 == 0:          if not v_norm2 == 0:
785                  v[iter+1]=v[iter+1]/h[iter+1][iter]              v[iter+1]=v[iter+1]/h[iter+1,iter]
786    
787          #   Form and store the information for the new Givens rotation          #   Form and store the information for the new Givens rotation
788      if iter > 0 :          if iter > 0 :
789          hhat=numarray.zeros(iter+1,numarray.Float64)              hhat=numpy.zeros(iter+1,numpy.float64)
790          for i in range(iter+1) : hhat[i]=h[i][iter]              for i in range(iter+1) : hhat[i]=h[i,iter]
791          hhat=__givapp(c[0:iter],s[0:iter],hhat);              hhat=__givapp(c[0:iter],s[0:iter],hhat);
792              for i in range(iter+1) : h[i][iter]=hhat[i]              for i in range(iter+1) : h[i,iter]=hhat[i]
793    
794      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])
795    
796      if mu!=0 :          if mu!=0 :
797          c[iter]=h[iter][iter]/mu              c[iter]=h[iter,iter]/mu
798          s[iter]=-h[iter+1][iter]/mu              s[iter]=-h[iter+1,iter]/mu
799          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]
800          h[iter+1][iter]=0.0              h[iter+1,iter]=0.0
801          g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])              gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
802                g[iter]=gg[0]
803                g[iter+1]=gg[1]
804    
805          # Update the residual norm          # Update the residual norm
806          rho=abs(g[iter+1])          rho=abs(g[iter+1])
807      iter+=1          iter+=1
808    
809     # At this point either iter > iter_max or rho < tol.     # At this point either iter > iter_max or rho < tol.
810     # It's time to compute x and leave.             # It's time to compute x and leave.
811     if iter > 0 :     if iter > 0 :
812       y=numarray.zeros(iter,numarray.Float64)           y=numpy.zeros(iter,numpy.float64)
813       y[iter-1] = g[iter-1] / h[iter-1][iter-1]       y[iter-1] = g[iter-1] / h[iter-1,iter-1]
814       if iter > 1 :         if iter > 1 :
815          i=iter-2            i=iter-2
816          while i>=0 :          while i>=0 :
817            y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]            y[i] = ( g[i] - numpy.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
818            i=i-1            i=i-1
819       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
820       for i in range(iter-1):       for i in range(iter-1):
821      xhat += v[i]*y[i]      xhat += v[i]*y[i]
822     else :     else :
823        xhat=v[0] * 0        xhat=v[0] * 0
824    
825     if iter<iter_restart-1:     if iter<iter_restart-1:
826        stopped=iter        stopped=iter
827     else:     else:
828        stopped=-1        stopped=-1
829    
830     return xhat,stopped     return xhat,stopped
831    
832  ##############################################  def GMRES(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100, iter_restart=20, verbose=False,P_R=None):
833  def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):     """
834  ################################################     Solver for
835    
836       *Ax=b*
837    
838       with a general operator A (more details required!).
839       It uses the generalized minimum residual method (GMRES).
840    
841       The iteration is terminated if
842    
843       *|r| <= atol+rtol*|r0|*
844    
845       where *r0* is the initial residual and *|.|* is the energy norm. In fact
846    
847       *|r| = sqrt( bilinearform(r,r))*
848    
849       :param r: initial residual *r=b-Ax*. ``r`` is altered.
850       :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
851       :param x: an initial guess for the solution
852       :type x: same like ``r``
853       :param Aprod: returns the value Ax
854       :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
855                    argument ``x``. The returned object needs to be of the same
856                    type like argument ``r``.
857       :param bilinearform: inner product ``<x,r>``
858       :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
859                           type like argument ``x`` and ``r``. The returned value is
860                           a ``float``.
861       :param atol: absolute tolerance
862       :type atol: non-negative ``float``
863       :param rtol: relative tolerance
864       :type rtol: non-negative ``float``
865       :param iter_max: maximum number of iteration steps
866       :type iter_max: ``int``
867       :param iter_restart: in order to save memory the orthogonalization process
868                            is terminated after ``iter_restart`` steps and the
869                            iteration is restarted.
870       :type iter_restart: ``int``
871       :return: the solution approximation and the corresponding residual
872       :rtype: ``tuple``
873       :warning: ``r`` and ``x`` are altered.
874       """
875     m=iter_restart     m=iter_restart
876       restarted=False
877     iter=0     iter=0
878     xc=x     if rtol>0:
879          r_dot_r = bilinearform(r, r)
880          if r_dot_r<0: raise NegativeNorm,"negative norm."
881          atol2=atol+rtol*math.sqrt(r_dot_r)
882          if verbose: print "GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)
883       else:
884          atol2=atol
885          if verbose: print "GMRES: absolute tolerance = %e"%atol2
886       if atol2<=0:
887          raise ValueError,"Non-positive tolarance."
888    
889     while True:     while True:
890        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
891        xc,stopped=__GMRESm(b*1, Aprod, Msolve, bilinearform, stoppingcriterium, x=xc*1, iter_max=iter_max-iter, iter_restart=m)        if restarted:
892             r2 = r-Aprod(x-x2)
893          else:
894             r2=1*r
895          x2=x*1.
896          x,stopped=_GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose,P_R=P_R)
897          iter+=iter_restart
898        if stopped: break        if stopped: break
899        iter+=iter_restart            if verbose: print "GMRES: restart."
900     return xc        restarted=True
901       if verbose: print "GMRES: tolerance has been reached."
902       return x
903    
904  def __GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):  def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False, P_R=None):
905     iter=0     iter=0
    r=Msolve(b)  
    r_dot_r = bilinearform(r, r)  
    if r_dot_r<0: raise NegativeNorm,"negative norm."  
    norm_b=math.sqrt(r_dot_r)  
906    
907     if x==None:     h=numpy.zeros((iter_restart+1,iter_restart),numpy.float64)
908        x=0*b     c=numpy.zeros(iter_restart,numpy.float64)
909     else:     s=numpy.zeros(iter_restart,numpy.float64)
910        r=Msolve(b-Aprod(x))     g=numpy.zeros(iter_restart+1,numpy.float64)
       r_dot_r = bilinearform(r, r)  
       if r_dot_r<0: raise NegativeNorm,"negative norm."  
     
    h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)  
    c=numarray.zeros(iter_restart,numarray.Float64)  
    s=numarray.zeros(iter_restart,numarray.Float64)  
    g=numarray.zeros(iter_restart,numarray.Float64)  
911     v=[]     v=[]
912    
913       r_dot_r = bilinearform(r, r)
914       if r_dot_r<0: raise NegativeNorm,"negative norm."
915     rho=math.sqrt(r_dot_r)     rho=math.sqrt(r_dot_r)
916      
917     v.append(r/rho)     v.append(r/rho)
918     g[0]=rho     g[0]=rho
919    
920     while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):     if verbose: print "GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)
921       while not (rho<=atol or iter==iter_restart):
922    
923      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
924    
925      p=Msolve(Aprod(v[iter]))          if P_R!=None:
926                p=Aprod(P_R(v[iter]))
927            else:
928            p=Aprod(v[iter])
929      v.append(p)      v.append(p)
930    
931      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
932    
933  # Modified Gram-Schmidt  # Modified Gram-Schmidt
934      for j in range(iter+1):      for j in range(iter+1):
935        h[j][iter]=bilinearform(v[j],v[iter+1])          h[j,iter]=bilinearform(v[j],v[iter+1])
936        v[iter+1]-=h[j][iter]*v[j]        v[iter+1]-=h[j,iter]*v[j]
937          
938      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]))
939      v_norm2=h[iter+1][iter]      v_norm2=h[iter+1,iter]
940    
941  # Reorthogonalize if needed  # Reorthogonalize if needed
942      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)
943       for j in range(iter+1):         for j in range(iter+1):
944          hr=bilinearform(v[j],v[iter+1])          hr=bilinearform(v[j],v[iter+1])
945              h[j][iter]=h[j][iter]+hr              h[j,iter]=h[j,iter]+hr
946              v[iter+1] -= hr*v[j]              v[iter+1] -= hr*v[j]
947    
948       v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
949       h[iter+1][iter]=v_norm2       h[iter+1,iter]=v_norm2
950    
951  #   watch out for happy breakdown  #   watch out for happy breakdown
952          if not v_norm2 == 0:          if not v_norm2 == 0:
953           v[iter+1]=v[iter+1]/h[iter+1][iter]           v[iter+1]=v[iter+1]/h[iter+1,iter]
954    
955  #   Form and store the information for the new Givens rotation  #   Form and store the information for the new Givens rotation
956      if iter > 0 :      if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])
957          hhat=numarray.zeros(iter+1,numarray.Float64)      mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
         for i in range(iter+1) : hhat[i]=h[i][iter]  
         hhat=__givapp(c[0:iter],s[0:iter],hhat);  
             for i in range(iter+1) : h[i][iter]=hhat[i]  
   
     mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])  
958    
959      if mu!=0 :      if mu!=0 :
960          c[iter]=h[iter][iter]/mu          c[iter]=h[iter,iter]/mu
961          s[iter]=-h[iter+1][iter]/mu          s[iter]=-h[iter+1,iter]/mu
962          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]
963          h[iter+1][iter]=0.0          h[iter+1,iter]=0.0
964          g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])                  gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
965                    g[iter]=gg[0]
966                    g[iter+1]=gg[1]
967  # Update the residual norm  # Update the residual norm
968                  
969          rho=abs(g[iter+1])          rho=abs(g[iter+1])
970            if verbose: print "GMRES: iteration step %s: residual %e"%(iter,rho)
971      iter+=1      iter+=1
972    
973  # At this point either iter > iter_max or rho < tol.  # At this point either iter > iter_max or rho < tol.
974  # It's time to compute x and leave.          # It's time to compute x and leave.
975    
976     if iter > 0 :     if verbose: print "GMRES: iteration stopped after %s step."%iter
977       y=numarray.zeros(iter,numarray.Float64)         if iter > 0 :
978       y[iter-1] = g[iter-1] / h[iter-1][iter-1]       y=numpy.zeros(iter,numpy.float64)
979       if iter > 1 :         y[iter-1] = g[iter-1] / h[iter-1,iter-1]
980          i=iter-2         if iter > 1 :
981            i=iter-2
982          while i>=0 :          while i>=0 :
983            y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]            y[i] = ( g[i] - numpy.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
984            i=i-1            i=i-1
985       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
986       for i in range(iter-1):       for i in range(iter-1):
987      xhat += v[i]*y[i]      xhat += v[i]*y[i]
988     else : xhat=v[0]     else:
989         xhat=v[0] * 0
990     x += xhat     if P_R!=None:
991     if iter<iter_restart-1:        x += P_R(xhat)
992        stopped=True     else:
993     else:        x += xhat
994       if iter<iter_restart-1:
995          stopped=True
996       else:
997        stopped=False        stopped=False
998    
999     return x,stopped     return x,stopped
1000    
1001  #################################################  def MINRES(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1002  def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):      """
1003  #################################################      Solver for
1004      #  
1005      #  minres solves the system of linear equations Ax = b      *Ax=b*
1006      #  where A is a symmetric matrix (possibly indefinite or singular)  
1007      #  and b is a given vector.      with a symmetric and positive definite operator A (more details required!).
1008      #        It uses the minimum residual method (MINRES) with preconditioner M
1009      #  "A" may be a dense or sparse matrix (preferably sparse!)      providing an approximation of A.
1010      #  or the name of a function such that  
1011      #               y = A(x)      The iteration is terminated if
1012      #  returns the product y = Ax for any given vector x.  
1013      #      *|r| <= atol+rtol*|r0|*
1014      #  "M" defines a positive-definite preconditioner M = C C'.  
1015      #  "M" may be a dense or sparse matrix (preferably sparse!)      where *r0* is the initial residual and *|.|* is the energy norm. In fact
1016      #  or the name of a function such that  
1017      #  solves the system My = x for any given vector x.      *|r| = sqrt( bilinearform(Msolve(r),r))*
1018      #  
1019      #      For details on the preconditioned conjugate gradient method see the book:
1020        
1021        I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
1022        T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
1023        C. Romine, and H. van der Vorst}.
1024    
1025        :param r: initial residual *r=b-Ax*. ``r`` is altered.
1026        :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1027        :param x: an initial guess for the solution
1028        :type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1029        :param Aprod: returns the value Ax
1030        :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
1031                     argument ``x``. The returned object needs to be of the same
1032                     type like argument ``r``.
1033        :param Msolve: solves Mx=r
1034        :type Msolve: function ``Msolve(r)`` where ``r`` is of the same type like
1035                      argument ``r``. The returned object needs to be of the same
1036                      type like argument ``x``.
1037        :param bilinearform: inner product ``<x,r>``
1038        :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
1039                            type like argument ``x`` and ``r`` is. The returned value
1040                            is a ``float``.
1041        :param atol: absolute tolerance
1042        :type atol: non-negative ``float``
1043        :param rtol: relative tolerance
1044        :type rtol: non-negative ``float``
1045        :param iter_max: maximum number of iteration steps
1046        :type iter_max: ``int``
1047        :return: the solution approximation and the corresponding residual
1048        :rtype: ``tuple``
1049        :warning: ``r`` and ``x`` are altered.
1050        """
1051      #------------------------------------------------------------------      #------------------------------------------------------------------
1052      # Set up y and v for the first Lanczos vector v1.      # Set up y and v for the first Lanczos vector v1.
1053      # y  =  beta1 P' v1,  where  P = C**(-1).      # y  =  beta1 P' v1,  where  P = C**(-1).
1054      # v is really P' v1.      # v is really P' v1.
1055      #------------------------------------------------------------------      #------------------------------------------------------------------
1056      if x==None:      r1    = r
1057        x=0*b      y = Msolve(r)
1058      else:      beta1 = bilinearform(y,r)
       b += (-1)*Aprod(x)  
1059    
     r1    = b  
     y = Msolve(b)  
     beta1 = bilinearform(y,b)  
   
1060      if beta1< 0: raise NegativeNorm,"negative norm."      if beta1< 0: raise NegativeNorm,"negative norm."
1061    
1062      #  If b = 0 exactly, stop with x = 0.      #  If r = 0 exactly, stop with x
1063      if beta1==0: return x*0.      if beta1==0: return x
1064    
1065      if beta1> 0:      if beta1> 0: beta1  = math.sqrt(beta1)
       beta1  = math.sqrt(beta1)        
1066    
1067      #------------------------------------------------------------------      #------------------------------------------------------------------
1068      # Initialize quantities.      # Initialize quantities.
# Line 1008  def MINRES(b, Aprod, Msolve, bilinearfor Line 1082  def MINRES(b, Aprod, Msolve, bilinearfor
1082      ynorm2 = 0      ynorm2 = 0
1083      cs     = -1      cs     = -1
1084      sn     = 0      sn     = 0
1085      w      = b*0.      w      = r*0.
1086      w2     = b*0.      w2     = r*0.
1087      r2     = r1      r2     = r1
1088      eps    = 0.0001      eps    = 0.0001
1089    
1090      #---------------------------------------------------------------------      #---------------------------------------------------------------------
1091      # Main iteration loop.      # Main iteration loop.
1092      # --------------------------------------------------------------------      # --------------------------------------------------------------------
1093      while not stoppingcriterium(rnorm,Anorm*ynorm,'MINRES'):    #  checks ||r|| < (||A|| ||x||) * TOL      while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL
1094    
1095      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
1096          iter    = iter  +  1          iter    = iter  +  1
# Line 1035  def MINRES(b, Aprod, Msolve, bilinearfor Line 1109  def MINRES(b, Aprod, Msolve, bilinearfor
1109          #-----------------------------------------------------------------          #-----------------------------------------------------------------
1110          s = 1/beta                 # Normalize previous vector (in y).          s = 1/beta                 # Normalize previous vector (in y).
1111          v = s*y                    # v = vk if P = I          v = s*y                    # v = vk if P = I
1112        
1113          y      = Aprod(v)          y      = Aprod(v)
1114        
1115          if iter >= 2:          if iter >= 2:
1116            y = y - (beta/oldb)*r1            y = y - (beta/oldb)*r1
1117    
1118          alfa   = bilinearform(v,y)              # alphak          alfa   = bilinearform(v,y)              # alphak
1119          y      += (- alfa/beta)*r2          y      += (- alfa/beta)*r2
1120          r1     = r2          r1     = r2
1121          r2     = y          r2     = y
1122          y = Msolve(r2)          y = Msolve(r2)
# Line 1052  def MINRES(b, Aprod, Msolve, bilinearfor Line 1126  def MINRES(b, Aprod, Msolve, bilinearfor
1126    
1127          beta   = math.sqrt( beta )          beta   = math.sqrt( beta )
1128          tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta          tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
1129            
1130          if iter==1:                 # Initialize a few things.          if iter==1:                 # Initialize a few things.
1131            gmax   = abs( alfa )      # alpha1            gmax   = abs( alfa )      # alpha1
1132            gmin   = gmax             # alpha1            gmin   = gmax             # alpha1
# Line 1060  def MINRES(b, Aprod, Msolve, bilinearfor Line 1134  def MINRES(b, Aprod, Msolve, bilinearfor
1134          # Apply previous rotation Qk-1 to get          # Apply previous rotation Qk-1 to get
1135          #   [deltak epslnk+1] = [cs  sn][dbark    0   ]          #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
1136          #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].          #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
1137        
1138          oldeps = epsln          oldeps = epsln
1139          delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak          delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
1140          gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k          gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
# Line 1070  def MINRES(b, Aprod, Msolve, bilinearfor Line 1144  def MINRES(b, Aprod, Msolve, bilinearfor
1144          # Compute the next plane rotation Qk          # Compute the next plane rotation Qk
1145    
1146          gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak          gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
1147          gamma  = max(gamma,eps)          gamma  = max(gamma,eps)
1148          cs     = gbar / gamma             # ck          cs     = gbar / gamma             # ck
1149          sn     = beta / gamma             # sk          sn     = beta / gamma             # sk
1150          phi    = cs * phibar              # phik          phi    = cs * phibar              # phik
# Line 1078  def MINRES(b, Aprod, Msolve, bilinearfor Line 1152  def MINRES(b, Aprod, Msolve, bilinearfor
1152    
1153          # Update  x.          # Update  x.
1154    
1155          denom = 1/gamma          denom = 1/gamma
1156          w1    = w2          w1    = w2
1157          w2    = w          w2    = w
1158          w     = (v - oldeps*w1 - delta*w2) * denom          w     = (v - oldeps*w1 - delta*w2) * denom
1159          x     +=  phi*w          x     +=  phi*w
1160    
# Line 1095  def MINRES(b, Aprod, Msolve, bilinearfor Line 1169  def MINRES(b, Aprod, Msolve, bilinearfor
1169    
1170          # Estimate various norms and test for convergence.          # Estimate various norms and test for convergence.
1171    
1172          Anorm  = math.sqrt( tnorm2 )          Anorm  = math.sqrt( tnorm2 )
1173          ynorm  = math.sqrt( ynorm2 )          ynorm  = math.sqrt( ynorm2 )
1174    
1175          rnorm  = phibar          rnorm  = phibar
1176    
1177      return x      return x
1178    
1179  def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  def TFQMR(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1180      """
1181      Solver for
1182    
1183  # TFQMR solver for linear systems    *Ax=b*
 #  
 #  
 # initialization  
 #  
   errtol = math.sqrt(bilinearform(b,b))  
   norm_b=errtol  
   kmax  = iter_max  
   error = []  
   
   if math.sqrt(bilinearform(x,x)) != 0.0:  
     r = b - Aprod(x)  
   else:  
     r = b  
1184    
1185    r=Msolve(r)    with a general operator A (more details required!).
1186      It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).
1187    
1188      The iteration is terminated if
1189    
1190      *|r| <= atol+rtol*|r0|*
1191    
1192      where *r0* is the initial residual and *|.|* is the energy norm. In fact
1193    
1194      *|r| = sqrt( bilinearform(r,r))*
1195    
1196      :param r: initial residual *r=b-Ax*. ``r`` is altered.
1197      :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1198      :param x: an initial guess for the solution
1199      :type x: same like ``r``
1200      :param Aprod: returns the value Ax
1201      :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
1202                   argument ``x``. The returned object needs to be of the same type
1203                   like argument ``r``.
1204      :param bilinearform: inner product ``<x,r>``
1205      :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
1206                          type like argument ``x`` and ``r``. The returned value is
1207                          a ``float``.
1208      :param atol: absolute tolerance
1209      :type atol: non-negative ``float``
1210      :param rtol: relative tolerance
1211      :type rtol: non-negative ``float``
1212      :param iter_max: maximum number of iteration steps
1213      :type iter_max: ``int``
1214      :rtype: ``tuple``
1215      :warning: ``r`` and ``x`` are altered.
1216      """
1217    u1=0    u1=0
1218    u2=0    u2=0
1219    y1=0    y1=0
1220    y2=0    y2=0
1221    
1222    w = r    w = r
1223    y1 = r    y1 = r
1224    iter = 0    iter = 0
1225    d = 0    d = 0
1226        v = Aprod(y1)
   v = Msolve(Aprod(y1))  
1227    u1 = v    u1 = v
1228      
1229    theta = 0.0;    theta = 0.0;
1230    eta = 0.0;    eta = 0.0;
1231    tau = math.sqrt(bilinearform(r,r))    rho=bilinearform(r,r)
1232    error = [ error, tau ]    if rho < 0: raise NegativeNorm,"negative norm."
1233    rho = tau * tau    tau = math.sqrt(rho)
1234    m=1    norm_r0=tau
1235  #    while tau>atol+rtol*norm_r0:
 #  TFQMR iteration  
 #  
 #  while ( iter < kmax-1 ):  
     
   while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):  
1236      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
1237    
1238      sigma = bilinearform(r,v)      sigma = bilinearform(r,v)
1239        if sigma == 0.0: raise IterationBreakDown,'TFQMR breakdown, sigma=0'
     if ( sigma == 0.0 ):  
       raise 'TFQMR breakdown, sigma=0'  
       
1240    
1241      alpha = rho / sigma      alpha = rho / sigma
1242    
# Line 1162  def TFQMR(b, Aprod, Msolve, bilinearform Line 1246  def TFQMR(b, Aprod, Msolve, bilinearform
1246  #  #
1247        if ( j == 1 ):        if ( j == 1 ):
1248          y2 = y1 - alpha * v          y2 = y1 - alpha * v
1249          u2 = Msolve(Aprod(y2))          u2 = Aprod(y2)
1250    
1251        m = 2 * (iter+1) - 2 + (j+1)        m = 2 * (iter+1) - 2 + (j+1)
1252        if j==0:        if j==0:
1253           w = w - alpha * u1           w = w - alpha * u1
1254           d = y1 + ( theta * theta * eta / alpha ) * d           d = y1 + ( theta * theta * eta / alpha ) * d
1255        if j==1:        if j==1:
# Line 1180  def TFQMR(b, Aprod, Msolve, bilinearform Line 1264  def TFQMR(b, Aprod, Msolve, bilinearform
1264  #  #
1265  #  Try to terminate the iteration at each pass through the loop  #  Try to terminate the iteration at each pass through the loop
1266  #  #
1267       # if ( tau * math.sqrt ( m + 1 ) <= errtol ):      if rho == 0.0: raise IterationBreakDown,'TFQMR breakdown, rho=0'
      #   error = [ error, tau ]  
      #   total_iters = iter  
      #   break  
         
   
     if ( rho == 0.0 ):  
       raise 'TFQMR breakdown, rho=0'  
       
1268    
1269      rhon = bilinearform(r,w)      rhon = bilinearform(r,w)
1270      beta = rhon / rho;      beta = rhon / rho;
1271      rho = rhon;      rho = rhon;
1272      y1 = w + beta * y2;      y1 = w + beta * y2;
1273      u1 = Msolve(Aprod(y1))      u1 = Aprod(y1)
1274      v = u1 + beta * ( u2 + beta * v )      v = u1 + beta * ( u2 + beta * v )
1275      error = [ error, tau ]  
1276      total_iters = iter      iter += 1
       
     iter = iter + 1  
1277    
1278    return x    return x
1279    
# Line 1208  def TFQMR(b, Aprod, Msolve, bilinearform Line 1282  def TFQMR(b, Aprod, Msolve, bilinearform
1282    
1283  class ArithmeticTuple(object):  class ArithmeticTuple(object):
1284     """     """
1285     tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.     Tuple supporting inplace update x+=y and scaling x=a*y where ``x,y`` is an
1286       ArithmeticTuple and ``a`` is a float.
1287    
1288     example of usage:     Example of usage::
1289    
1290     from esys.escript import Data         from esys.escript import Data
1291     from numarray import array         from numpy import array
1292     a=Data(...)         a=Data(...)
1293     b=array([1.,4.])         b=array([1.,4.])
1294     x=ArithmeticTuple(a,b)         x=ArithmeticTuple(a,b)
1295     y=5.*x         y=5.*x
1296    
1297     """     """
1298     def __init__(self,*args):     def __init__(self,*args):
1299         """         """
1300         initialize object with elements args.         Initializes object with elements ``args``.
1301    
1302         @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)         :param args: tuple of objects that support inplace add (x+=y) and
1303                        scaling (x=a*y)
1304         """         """
1305         self.__items=list(args)         self.__items=list(args)
1306    
1307     def __len__(self):     def __len__(self):
1308         """         """
1309         number of items         Returns the number of items.
1310    
1311         @return: number of items         :return: number of items
1312         @rtype: C{int}         :rtype: ``int``
1313         """         """
1314         return len(self.__items)         return len(self.__items)
1315    
1316     def __getitem__(self,index):     def __getitem__(self,index):
1317         """         """
1318         get an item         Returns item at specified position.
1319    
1320         @param index: item to be returned         :param index: index of item to be returned
1321         @type index: C{int}         :type index: ``int``
1322         @return: item with index C{index}         :return: item with index ``index``
1323         """         """
1324         return self.__items.__getitem__(index)         return self.__items.__getitem__(index)
1325    
1326     def __mul__(self,other):     def __mul__(self,other):
1327         """         """
1328         scaling from the right         Scales by ``other`` from the right.
1329    
1330         @param other: scaling factor         :param other: scaling factor
1331         @type other: C{float}         :type other: ``float``
1332         @return: itemwise self*other         :return: itemwise self*other
1333         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1334         """         """
1335         out=[]         out=[]
1336         try:           try:
1337             l=len(other)             l=len(other)
1338             if l!=len(self):             if l!=len(self):
1339                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1340             for i in range(l): out.append(self[i]*other[i])             for i in range(l): out.append(self[i]*other[i])
1341         except TypeError:         except TypeError:
1342         for i in range(len(self)): out.append(self[i]*other)             for i in range(len(self)): out.append(self[i]*other)
1343         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1344    
1345     def __rmul__(self,other):     def __rmul__(self,other):
1346         """         """
1347         scaling from the left         Scales by ``other`` from the left.
1348    
1349         @param other: scaling factor         :param other: scaling factor
1350         @type other: C{float}         :type other: ``float``
1351         @return: itemwise other*self         :return: itemwise other*self
1352         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1353         """         """
1354         out=[]         out=[]
1355         try:           try:
1356             l=len(other)             l=len(other)
1357             if l!=len(self):             if l!=len(self):
1358                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1359             for i in range(l): out.append(other[i]*self[i])             for i in range(l): out.append(other[i]*self[i])
1360         except TypeError:         except TypeError:
1361         for i in range(len(self)): out.append(other*self[i])             for i in range(len(self)): out.append(other*self[i])
1362         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1363    
1364     def __div__(self,other):     def __div__(self,other):
1365         """         """
1366         dividing from the right         Scales by (1/``other``) from the right.
1367    
1368         @param other: scaling factor         :param other: scaling factor
1369         @type other: C{float}         :type other: ``float``
1370         @return: itemwise self/other         :return: itemwise self/other
1371         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1372         """         """
1373         return self*(1/other)         return self*(1/other)
1374    
1375     def __rdiv__(self,other):     def __rdiv__(self,other):
1376         """         """
1377         dividing from the left         Scales by (1/``other``) from the left.
1378    
1379         @param other: scaling factor         :param other: scaling factor
1380         @type other: C{float}         :type other: ``float``
1381         @return: itemwise other/self         :return: itemwise other/self
1382         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1383         """         """
1384         out=[]         out=[]
1385         try:           try:
1386             l=len(other)             l=len(other)
1387             if l!=len(self):             if l!=len(self):
1388                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1389             for i in range(l): out.append(other[i]/self[i])             for i in range(l): out.append(other[i]/self[i])
1390         except TypeError:         except TypeError:
1391         for i in range(len(self)): out.append(other/self[i])             for i in range(len(self)): out.append(other/self[i])
1392         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1393      
1394     def __iadd__(self,other):     def __iadd__(self,other):
1395         """         """
1396         in-place add of other to self         Inplace addition of ``other`` to self.
1397    
1398         @param other: increment         :param other: increment
1399         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1400         """         """
1401         if len(self) != len(other):         if len(self) != len(other):
1402             raise ValueError,"tuple length must match."             raise ValueError,"tuple lengths must match."
1403         for i in range(len(self)):         for i in range(len(self)):
1404             self.__items[i]+=other[i]             self.__items[i]+=other[i]
1405         return self         return self
1406    
1407     def __add__(self,other):     def __add__(self,other):
1408         """         """
1409         add other to self         Adds ``other`` to self.
1410    
1411         @param other: increment         :param other: increment
1412         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1413         """         """
1414         out=[]         out=[]
1415         try:           try:
1416             l=len(other)             l=len(other)
1417             if l!=len(self):             if l!=len(self):
1418                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1419             for i in range(l): out.append(self[i]+other[i])             for i in range(l): out.append(self[i]+other[i])
1420         except TypeError:         except TypeError:
1421         for i in range(len(self)): out.append(self[i]+other)             for i in range(len(self)): out.append(self[i]+other)
1422         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1423    
1424     def __sub__(self,other):     def __sub__(self,other):
1425         """         """
1426         subtract other from self         Subtracts ``other`` from self.
1427    
1428         @param other: increment         :param other: decrement
1429         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1430         """         """
1431         out=[]         out=[]
1432         try:           try:
1433             l=len(other)             l=len(other)
1434             if l!=len(self):             if l!=len(self):
1435                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1436             for i in range(l): out.append(self[i]-other[i])             for i in range(l): out.append(self[i]-other[i])
1437         except TypeError:         except TypeError:
1438         for i in range(len(self)): out.append(self[i]-other)             for i in range(len(self)): out.append(self[i]-other)
1439         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1440      
1441     def __isub__(self,other):     def __isub__(self,other):
1442         """         """
1443         subtract other from self         Inplace subtraction of ``other`` from self.
1444    
1445         @param other: increment         :param other: decrement
1446         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1447         """         """
1448         if len(self) != len(other):         if len(self) != len(other):
1449             raise ValueError,"tuple length must match."             raise ValueError,"tuple length must match."
# Line 1377  class ArithmeticTuple(object): Line 1453  class ArithmeticTuple(object):
1453    
1454     def __neg__(self):     def __neg__(self):
1455         """         """
1456         negate         Negates values.
   
1457         """         """
1458         out=[]         out=[]
1459         for i in range(len(self)):         for i in range(len(self)):
# Line 1388  class ArithmeticTuple(object): Line 1463  class ArithmeticTuple(object):
1463    
1464  class HomogeneousSaddlePointProblem(object):  class HomogeneousSaddlePointProblem(object):
1465        """        """
1466        This provides a framwork for solving linear homogeneous saddle point problem of the form        This class provides a framework for solving linear homogeneous saddle
1467          point problems of the form::
              Av+B^*p=f  
              Bv    =0  
1468    
1469        for the unknowns v and p and given operators A and B and given right hand side f.            *Av+B^*p=f*
1470        B^* is the adjoint operator of B is the given inner product.            *Bv     =0*
1471    
1472          for the unknowns *v* and *p* and given operators *A* and *B* and
1473          given right hand side *f*. *B^** is the adjoint operator of *B*.
1474        """        """
1475        def __init__(self,**kwargs):        def __init__(self, adaptSubTolerance=True, **kwargs):
1476        """
1477        initializes the saddle point problem
1478        
1479        :param adaptSubTolerance: If True the tolerance for subproblem is set automatically.
1480        :type adaptSubTolerance: ``bool``
1481        """
1482          self.setTolerance()          self.setTolerance()
1483          self.setToleranceReductionFactor()          self.setAbsoluteTolerance()
1484        self.__adaptSubTolerance=adaptSubTolerance
1485          #=============================================================
1486        def initialize(self):        def initialize(self):
1487          """          """
1488          initialize the problem (overwrite)          Initializes the problem (overwrite).
1489          """          """
1490          pass          pass
1491        def B(self,v):  
1492          def inner_pBv(self,p,Bv):
1493           """           """
1494           returns Bv (overwrite)           Returns inner product of element p and Bv (overwrite).
          @rtype: equal to the type of p  
1495    
1496           @note: boundary conditions on p should be zero!           :param p: a pressure increment
1497             :param Bv: a residual
1498             :return: inner product of element p and Bv
1499             :rtype: ``float``
1500             :note: used if PCG is applied.
1501           """           """
1502           pass           raise NotImplementedError,"no inner product for p and Bv implemented."
1503    
1504        def inner(self,p0,p1):        def inner_p(self,p0,p1):
1505           """           """
1506           returns inner product of two element p0 and p1  (overwrite)           Returns inner product of p0 and p1 (overwrite).
           
          @type p0: equal to the type of p  
          @type p1: equal to the type of p  
          @rtype: C{float}  
1507    
1508           @rtype: equal to the type of p           :param p0: a pressure
1509             :param p1: a pressure
1510             :return: inner product of p0 and p1
1511             :rtype: ``float``
1512           """           """
1513           pass           raise NotImplementedError,"no inner product for p implemented."
1514      
1515          def norm_v(self,v):
1516             """
1517             Returns the norm of v (overwrite).
1518    
1519        def solve_A(self,u,p):           :param v: a velovity
1520             :return: norm of v
1521             :rtype: non-negative ``float``
1522           """           """
1523           solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)           raise NotImplementedError,"no norm of v implemented."
1524          def getV(self, p, v0):
1525             """
1526             return the value for v for a given p (overwrite)
1527    
1528           @rtype: equal to the type of v           :param p: a pressure
1529           @note: boundary conditions on v should be zero!           :param v0: a initial guess for the value v to return.
1530             :return: v given as *v= A^{-1} (f-B^*p)*
1531           """           """
1532           pass           raise NotImplementedError,"no v calculation implemented."
1533    
1534            
1535          def Bv(self,v):
1536            """
1537            Returns Bv (overwrite).
1538    
1539            :rtype: equal to the type of p
1540            :note: boundary conditions on p should be zero!
1541            """
1542            raise NotImplementedError, "no operator B implemented."
1543    
1544          def norm_Bv(self,Bv):
1545            """
1546            Returns the norm of Bv (overwrite).
1547    
1548        def solve_prec(self,p):          :rtype: equal to the type of p
1549            :note: boundary conditions on p should be zero!
1550            """
1551            raise NotImplementedError, "no norm of Bv implemented."
1552    
1553          def solve_AinvBt(self,p):
1554           """           """
1555           provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)           Solves *Av=B^*p* with accuracy `self.getSubProblemTolerance()`
1556             (overwrite).
1557    
1558           @rtype: equal to the type of p           :param p: a pressure increment
1559             :return: the solution of *Av=B^*p*
1560             :note: boundary conditions on v should be zero!
1561           """           """
1562           pass           raise NotImplementedError,"no operator A implemented."
1563    
1564        def stoppingcriterium(self,Bv,v,p):        def solve_prec(self,Bv):
1565           """           """
1566           returns a True if iteration is terminated. (overwrite)           Provides a preconditioner for *(BA^{-1}B^ * )* applied to Bv with accuracy
1567             `self.getSubProblemTolerance()` (overwrite).
1568    
1569           @rtype: C{bool}           :rtype: equal to the type of p
1570             :note: boundary conditions on p should be zero!
1571             """
1572             raise NotImplementedError,"no preconditioner for Schur complement implemented."
1573          def setSubProblemTolerance(self):
1574             """
1575         Updates the tolerance for subproblems
1576         :note: method is typically the method is overwritten.
1577           """           """
1578           pass           pass
1579                      #=============================================================
1580        def __inner(self,p,r):        def __Aprod_PCG(self,p):
1581           return self.inner(p,r[1])            dv=self.solve_AinvBt(p)
1582              return ArithmeticTuple(dv,self.Bv(dv))
       def __inner_p(self,p1,p2):  
          return self.inner(p1,p2)  
         
       def __inner_a(self,a1,a2):  
          return self.inner_a(a1,a2)  
   
       def __inner_a1(self,a1,a2):  
          return self.inner(a1[1],a2[1])  
1583    
1584        def __stoppingcriterium(self,norm_r,r,p):        def __inner_PCG(self,p,r):
1585            return self.stoppingcriterium(r[1],r[0],p)           return self.inner_pBv(p,r[1])
1586    
1587        def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):        def __Msolve_PCG(self,r):
1588            return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)            return self.solve_prec(r[1])
1589          #=============================================================
1590          def __Aprod_GMRES(self,p):
1591              return self.solve_prec(self.Bv(self.solve_AinvBt(p)))
1592          def __inner_GMRES(self,p0,p1):
1593             return self.inner_p(p0,p1)
1594    
1595          #=============================================================
1596          def norm_p(self,p):
1597              """
1598              calculates the norm of ``p``
1599              
1600              :param p: a pressure
1601              :return: the norm of ``p`` using the inner product for pressure
1602              :rtype: ``float``
1603              """
1604              f=self.inner_p(p,p)
1605              if f<0: raise ValueError,"negative pressure norm."
1606              return math.sqrt(f)
1607          def adaptSubTolerance(self):
1608          """
1609          Returns True if tolerance adaption for subproblem is choosen.
1610          """
1611              return self.__adaptSubTolerance
1612          
1613          def solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1614             """
1615             Solves the saddle point problem using initial guesses v and p.
1616    
1617        def setTolerance(self,tolerance=1.e-8):           :param v: initial guess for velocity
1618                self.__tol=tolerance           :param p: initial guess for pressure
1619        def getTolerance(self):           :type v: `Data`
1620                return self.__tol           :type p: `Data`
1621        def setToleranceReductionFactor(self,reduction=0.01):           :param usePCG: indicates the usage of the PCG rather than GMRES scheme.
1622                self.__reduction=reduction           :param max_iter: maximum number of iteration steps per correction
1623        def getSubProblemTolerance(self):                            attempt
1624                return self.__reduction*self.getTolerance()           :param verbose: if True, shows information on the progress of the
1625                             saddlepoint problem solver.
1626             :param iter_restart: restart the iteration after ``iter_restart`` steps
1627                                  (only used if useUzaw=False)
1628             :type usePCG: ``bool``
1629             :type max_iter: ``int``
1630             :type verbose: ``bool``
1631             :type iter_restart: ``int``
1632             :rtype: ``tuple`` of `Data` objects
1633             """
1634             self.verbose=verbose
1635             rtol=self.getTolerance()
1636             atol=self.getAbsoluteTolerance()
1637         if self.adaptSubTolerance(): self.setSubProblemTolerance()
1638             correction_step=0
1639             converged=False
1640             while not converged:
1641                  # calculate velocity for current pressure:
1642                  v=self.getV(p,v)
1643                  Bv=self.Bv(v)
1644                  norm_v=self.norm_v(v)
1645                  norm_Bv=self.norm_Bv(Bv)
1646                  ATOL=norm_v*rtol+atol
1647                  if self.verbose: print "HomogeneousSaddlePointProblem: norm v= %e, norm_Bv= %e, tolerance = %e."%(norm_v, norm_Bv,ATOL)
1648                  if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."
1649                  if norm_Bv <= ATOL:
1650                     converged=True
1651                  else:
1652                     correction_step+=1
1653                     if correction_step>max_correction_steps:
1654                          raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1655                     dp=self.solve_prec(Bv)
1656                     if usePCG:
1657                       norm2=self.inner_pBv(dp,Bv)
1658                       if norm2<0: raise ValueError,"negative PCG norm."
1659                       norm2=math.sqrt(norm2)
1660                     else:
1661                       norm2=self.norm_p(dp)
1662                     ATOL_ITER=ATOL/norm_Bv*norm2*0.5
1663                     if self.verbose: print "HomogeneousSaddlePointProblem: tolerance for solver: %e"%ATOL_ITER
1664                     if usePCG:
1665                           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)
1666                     else:
1667                           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)
1668             if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached."
1669         return v,p
1670    
1671        def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):        #========================================================================
1672                """        def setTolerance(self,tolerance=1.e-4):
1673                solves the saddle point problem using initial guesses v and p.           """
1674             Sets the relative tolerance for (v,p).
               @param max_iter: maximum number of iteration steps.  
               """  
               self.verbose=verbose  
               self.show_details=show_details and self.verbose  
   
               # assume p is known: then v=A^-1(f-B^*p)  
               # which leads to BA^-1B^*p = BA^-1f    
   
           # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)        
           self.__z=v+self.solve_A(v,p*0)  
               Bz=self.B(self.__z)  
               #  
           #   solve BA^-1B^*p = Bz  
               #  
               #  
               #  
               self.iter=0  
           if solver=='GMRES':        
                 if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter  
                 p=GMRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)  
                 # solve Au=f-B^*p  
                 #       A(u-v)=f-B^*p-Av  
                 #       u=v+(u-v)  
         u=v+self.solve_A(v,p)  
   
           if solver=='TFQMR':        
                 if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter  
                 p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)  
                 # solve Au=f-B^*p  
                 #       A(u-v)=f-B^*p-Av  
                 #       u=v+(u-v)  
         u=v+self.solve_A(v,p)  
   
           if solver=='MINRES':        
                 if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter  
                 p=MINRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)  
                 # solve Au=f-B^*p  
                 #       A(u-v)=f-B^*p-Av  
                 #       u=v+(u-v)  
         u=v+self.solve_A(v,p)  
                 
           if solver=='GMRESC':        
                 if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter  
                 p0=self.solve_prec1(Bz)  
             #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))  
                 #p-=alfa  
                 x=GMRES(ArithmeticTuple(self.__z*1.,p0*1),self.__Anext,self.__Mempty,self.__inner_a,self.__stoppingcriterium2,iter_max=max_iter, x=ArithmeticTuple(v*1,p*1),iter_restart=20)  
                 #x=NewtonGMRES(ArithmeticTuple(self.__z*1.,p0*1),self.__Aprod_Newton2,self.__Mempty,self.__inner_a,self.__stoppingcriterium2,iter_max=max_iter, x=ArithmeticTuple(v*1,p*1),atol=0,rtol=self.getTolerance())  
   
                 # solve Au=f-B^*p  
                 #       A(u-v)=f-B^*p-Av  
                 #       u=v+(u-v)  
             p=x[1]  
         u=v+self.solve_A(v,p)        
         #p=x[1]  
         #u=x[0]  
   
               if solver=='PCG':  
                 #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv  
                 #  
                 #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)  
                 #                           A(v-z)= f -Az - B^*p (v-z=0 on fixed_u_mask)  
                 if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter  
                 p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)  
             u=r[0]    
                 # print "DIFF=",util.integrate(p)  
1675    
1676                # print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)           :param tolerance: tolerance to be used
1677             :type tolerance: non-negative ``float``
1678             """
1679             if tolerance<0:
1680                 raise ValueError,"tolerance must be positive."
1681             self.__rtol=tolerance
1682    
1683            return u,p        def getTolerance(self):
1684             """
1685             Returns the relative tolerance.
1686    
1687        def __Msolve(self,r):           :return: relative tolerance
1688            return self.solve_prec(r[1])           :rtype: ``float``
1689             """
1690             return self.__rtol
1691    
1692        def __Msolve2(self,r):        def setAbsoluteTolerance(self,tolerance=0.):
1693            return self.solve_prec(r*1)           """
1694             Sets the absolute tolerance.
1695    
1696        def __Mempty(self,r):           :param tolerance: tolerance to be used
1697            return r           :type tolerance: non-negative ``float``
1698             """
1699             if tolerance<0:
1700                 raise ValueError,"tolerance must be non-negative."
1701             self.__atol=tolerance
1702    
1703          def getAbsoluteTolerance(self):
1704             """
1705             Returns the absolute tolerance.
1706    
1707        def __Aprod(self,p):           :return: absolute tolerance
1708            # return BA^-1B*p           :rtype: ``float``
1709            #solve Av =B^*p as Av =f-Az-B^*(-p)           """
1710            v=self.solve_A(self.__z,-p)           return self.__atol
           return ArithmeticTuple(v, self.B(v))  
   
       def __Anext(self,x):  
           # return next v,p  
           #solve Av =-B^*p as Av =f-Az-B^*p  
   
       pc=x[1]  
           v=self.solve_A(self.__z,-pc)  
       p=self.solve_prec1(self.B(v))  
   
           return ArithmeticTuple(v,p)  
   
   
       def __Aprod2(self,p):  
           # return BA^-1B*p  
           #solve Av =B^*p as Av =f-Az-B^*(-p)  
       v=self.solve_A(self.__z,-p)  
           return self.B(v)  
   
       def __Aprod_Newton(self,p):  
           # return BA^-1B*p - Bz  
           #solve Av =-B^*p as Av =f-Az-B^*p  
       v=self.solve_A(self.__z,-p)  
           return self.B(v-self.__z)  
   
       def __Aprod_Newton2(self,x):  
           # return BA^-1B*p - Bz  
           #solve Av =-B^*p as Av =f-Az-B^*p  
           pc=x[1]  
       v=self.solve_A(self.__z,-pc)  
           p=self.solve_prec1(self.B(v-self.__z))  
           return ArithmeticTuple(v,p)  
1711    
1712          def getSubProblemTolerance(self):
1713             """
1714             Sets the relative tolerance to solve the subproblem(s).
1715             """
1716             return max(200.*util.EPSILON,self.getTolerance()**2)
1717    
1718  def MaskFromBoundaryTag(domain,*tags):  def MaskFromBoundaryTag(domain,*tags):
1719     """     """
1720     creates a mask on the Solution(domain) function space which one for samples     Creates a mask on the Solution(domain) function space where the value is
1721     that touch the boundary tagged by tags.     one for samples that touch the boundary tagged by tags.
1722    
1723     usage: m=MaskFromBoundaryTag(domain,"left", "right")     Usage: m=MaskFromBoundaryTag(domain, "left", "right")
1724    
1725     @param domain: a given domain     :param domain: domain to be used
1726     @type domain: L{escript.Domain}     :type domain: `escript.Domain`
1727     @param tags: boundray tags     :param tags: boundary tags
1728     @type tags: C{str}     :type tags: ``str``
1729     @return: a mask which marks samples that are touching the boundary tagged by any of the given tags.     :return: a mask which marks samples that are touching the boundary tagged
1730     @rtype: L{escript.Data} of rank 0              by any of the given tags
1731       :rtype: `escript.Data` of rank 0
1732     """     """
1733     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1734     d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))     d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
1735     for t in tags: d.setTaggedValue(t,1.)     for t in tags: d.setTaggedValue(t,1.)
1736     pde.setValue(y=d)     pde.setValue(y=d)
1737     return util.whereNonZero(pde.getRightHandSide())     return util.whereNonZero(pde.getRightHandSide())
 #============================================================================================================================================  
 class SaddlePointProblem(object):  
    """  
    This implements a solver for a saddlepoint problem  
   
    M{f(u,p)=0}  
    M{g(u)=0}  
   
    for u and p. The problem is solved with an inexact Uszawa scheme for p:  
   
    M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}  
    M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}  
1738    
1739     where Q_f is an approximation of the Jacobiean A_f of f with respect to u  and Q_f is an approximation of  def MaskFromTag(domain,*tags):
    A_g A_f^{-1} A_g with A_g is the jacobiean of g with respect to p. As a the construction of a 'proper'  
    Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays  
    in fact the role of a preconditioner.  
1740     """     """
1741     def __init__(self,verbose=False,*args):     Creates a mask on the Solution(domain) function space where the value is
1742         """     one for samples that touch regions tagged by tags.
        initializes the problem  
   
        @param verbose: switches on the printing out some information  
        @type verbose: C{bool}  
        @note: this method may be overwritten by a particular saddle point problem  
        """  
        print "SaddlePointProblem should not be used anymore!"  
        if not isinstance(verbose,bool):  
             raise TypeError("verbose needs to be of type bool.")  
        self.__verbose=verbose  
        self.relaxation=1.  
        DeprecationWarning("SaddlePointProblem should not be used anymore.")  
   
    def trace(self,text):  
        """  
        prints text if verbose has been set  
   
        @param text: a text message  
        @type text: C{str}  
        """  
        if self.__verbose: print "%s: %s"%(str(self),text)  
   
    def solve_f(self,u,p,tol=1.e-8):  
        """  
        solves  
1743    
1744         A_f du = f(u,p)     Usage: m=MaskFromTag(domain, "ham")
1745    
1746         with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.     :param domain: domain to be used
1747       :type domain: `escript.Domain`
1748         @param u: current approximation of u     :param tags: boundary tags
1749         @type u: L{escript.Data}     :type tags: ``str``
1750         @param p: current approximation of p     :return: a mask which marks samples that are touching the boundary tagged
1751         @type p: L{escript.Data}              by any of the given tags
1752         @param tol: tolerance expected for du     :rtype: `escript.Data` of rank 0
1753         @type tol: C{float}     """
1754         @return: increment du     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1755         @rtype: L{escript.Data}     d=escript.Scalar(0.,escript.Function(domain))
1756         @note: this method has to be overwritten by a particular saddle point problem     for t in tags: d.setTaggedValue(t,1.)
1757         """     pde.setValue(Y=d)
1758         pass     return util.whereNonZero(pde.getRightHandSide())
   
    def solve_g(self,u,tol=1.e-8):  
        """  
        solves  
   
        Q_g dp = g(u)  
   
        with Q_g is a preconditioner for A_g A_f^{-1} A_g with  A_g is the jacobiean of g with respect to p.  
   
        @param u: current approximation of u  
        @type u: L{escript.Data}  
        @param tol: tolerance expected for dp  
        @type tol: C{float}  
        @return: increment dp  
        @rtype: L{escript.Data}  
        @note: this method has to be overwritten by a particular saddle point problem  
        """  
        pass  
   
    def inner(self,p0,p1):  
        """  
        inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)  
        @return: inner product of p0 and p1  
        @rtype: C{float}  
        """  
        pass  
1759    
    subiter_max=3  
    def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):  
         """  
         runs the solver  
1760    
         @param u0: initial guess for C{u}  
         @type u0: L{esys.escript.Data}  
         @param p0: initial guess for C{p}  
         @type p0: L{esys.escript.Data}  
         @param tolerance: tolerance for relative error in C{u} and C{p}  
         @type tolerance: positive C{float}  
         @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}  
         @type tolerance_u: positive C{float}  
         @param iter_max: maximum number of iteration steps.  
         @type iter_max: C{int}  
         @param accepted_reduction: if the norm  g cannot be reduced by C{accepted_reduction} backtracking to adapt the  
                                    relaxation factor. If C{accepted_reduction=None} no backtracking is used.  
         @type accepted_reduction: positive C{float} or C{None}  
         @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.  
         @type relaxation: C{float} or C{None}  
         """  
         tol=1.e-2  
         if tolerance_u==None: tolerance_u=tolerance  
         if not relaxation==None: self.relaxation=relaxation  
         if accepted_reduction ==None:  
               angle_limit=0.  
         elif accepted_reduction>=1.:  
               angle_limit=0.  
         else:  
               angle_limit=util.sqrt(1-accepted_reduction**2)  
         self.iter=0  
         u=u0  
         p=p0  
         #  
         #   initialize things:  
         #  
         converged=False  
         #  
         #  start loop:  
         #  
         #  initial search direction is g  
         #  
         while not converged :  
             if self.iter>iter_max:  
                 raise ArithmeticError("no convergence after %s steps."%self.iter)  
             f_new=self.solve_f(u,p,tol)  
             norm_f_new = util.Lsup(f_new)  
             u_new=u-f_new  
             g_new=self.solve_g(u_new,tol)  
             self.iter+=1  
             norm_g_new = util.sqrt(self.inner(g_new,g_new))  
             if norm_f_new==0. and norm_g_new==0.: return u, p  
             if self.iter>1 and not accepted_reduction==None:  
                #  
                #   did we manage to reduce the norm of G? I  
                #   if not we start a backtracking procedure  
                #  
                # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g  
                if norm_g_new > accepted_reduction * norm_g:  
                   sub_iter=0  
                   s=self.relaxation  
                   d=g  
                   g_last=g  
                   self.trace("    start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))  
                   while sub_iter < self.subiter_max and  norm_g_new > accepted_reduction * norm_g:  
                      dg= g_new-g_last  
                      norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)  
                      rad=self.inner(g_new,dg)/self.relaxation  
                      # print "   ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit  
                      # print "   ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g  
                      if abs(rad) < norm_dg*norm_g_new * angle_limit:  
                          if sub_iter>0: self.trace("    no further improvements expected from backtracking.")  
                          break  
                      r=self.relaxation  
                      self.relaxation= - rad/norm_dg**2  
                      s+=self.relaxation  
                      #####  
                      # a=g_new+self.relaxation*dg/r  
                      # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation  
                      #####  
                      g_last=g_new  
                      p+=self.relaxation*d  
                      f_new=self.solve_f(u,p,tol)  
                      u_new=u-f_new  
                      g_new=self.solve_g(u_new,tol)  
                      self.iter+=1  
                      norm_f_new = util.Lsup(f_new)  
                      norm_g_new = util.sqrt(self.inner(g_new,g_new))  
                      # print "   ",sub_iter," new g norm",norm_g_new  
                      self.trace("    %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))  
                      #  
                      #   can we expect reduction of g?  
                      #  
                      # u_last=u_new  
                      sub_iter+=1  
                   self.relaxation=s  
             #  
             #  check for convergence:  
             #  
             norm_u_new = util.Lsup(u_new)  
             p_new=p+self.relaxation*g_new  
             norm_p_new = util.sqrt(self.inner(p_new,p_new))  
             self.trace("%s th step: f = %s, f/u = %s, g = %s, g/p = %s, relaxation = %s."%(self.iter, norm_f_new ,norm_f_new/norm_u_new, norm_g_new, norm_g_new/norm_p_new, self.relaxation))  
   
             if self.iter>1:  
                dg2=g_new-g  
                df2=f_new-f  
                norm_dg2=util.sqrt(self.inner(dg2,dg2))  
                norm_df2=util.Lsup(df2)  
                # print norm_g_new, norm_g, norm_dg, norm_p, tolerance  
                tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new  
                tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new  
                if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:  
                    converged=True  
             f, norm_f, u, norm_u, g, norm_g, p, norm_p = f_new, norm_f_new, u_new, norm_u_new, g_new, norm_g_new, p_new, norm_p_new  
         self.trace("convergence after %s steps."%self.iter)  
         return u,p  

Legend:
Removed from v.1956  
changed lines
  Added in v.2647

  ViewVC Help
Powered by ViewVC 1.1.26