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

  ViewVC Help
Powered by ViewVC 1.1.26