/[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 1552 by gross, Thu May 8 08:52:41 2008 UTC revision 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC
# Line 1  Line 1 
1    
2    ########################################################
3  #  #
4  # $Id$  # Copyright (c) 2003-2009 by University of Queensland
5  #  # Earth Systems Science Computational Center (ESSCC)
6  #######################################################  # http://www.uq.edu.au/esscc
7  #  #
8  #           Copyright 2003-2007 by ACceSS MNRF  # Primary Business: Queensland, Australia
9  #       Copyright 2007 by University of Queensland  # Licensed under the Open Software License version 3.0
10  #  # http://www.opensource.org/licenses/osl-3.0.php
11  #                http://esscc.uq.edu.au  #
12  #        Primary Business: Queensland, Australia  ########################################################
13  #  Licensed under the Open Software License version 3.0  
14  #     http://www.opensource.org/licenses/osl-3.0.php  __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15  #  Earth Systems Science Computational Center (ESSCC)
16  #######################################################  http://www.uq.edu.au/esscc
17  #  Primary Business: Queensland, Australia"""
18    __license__="""Licensed under the Open Software License version 3.0
19    http://www.opensource.org/licenses/osl-3.0.php"""
20    __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 32  Currently includes: Line 37  Currently includes:
37  """  """
38    
39  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
 __copyright__="""  Copyright (c) 2006 by ACcESS MNRF  
                     http://www.access.edu.au  
                 Primary Business: Queensland, Australia"""  
 __license__="""Licensed under the Open Software License version 3.0  
              http://www.opensource.org/licenses/osl-3.0.php"""  
 __url__="http://www.iservo.edu.au/esys"  
 __version__="$Revision$"  
 __date__="$Date$"  
40    
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 55  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 71  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 88  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 97  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 114  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 129  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 189  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 237  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 279  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 299  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 353  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 378  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 415  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       if verbose: print "PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_r0, atol2)
528    
529    
530     while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):     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 560  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       if verbose: print "PCG: tolerance reached after %s steps."%iter
551       return x,r,math.sqrt(rhat_dot_r)
552    
553     return x,r  class Defect(object):
554        """
555        Defines a non-linear defect F(x) of a variable x.
556        """
557        def __init__(self):
558            """
559            Initializes defect.
560            """
561            self.setDerivativeIncrementLength()
562    
563        def bilinearform(self, x0, x1):
564            """
565            Returns the inner product of x0 and x1
566    
567  ############################          @param x0: value for x0
568  # Added by Artak          @param x1: value for x1
569  #################################3          @return: the inner product of x0 and x1
570            @rtype: C{float}
571            """
572            return 0
573    
574  #Apply a sequence of k Givens rotations, used within gmres codes      def norm(self,x):
575  # vrot=givapp(c, s, vin, k)          """
576  def givapp(c,s,vin):          Returns the norm of argument C{x}.
     vrot=vin # warning: vin is altered!!!!  
     if isinstance(c,float):  
         vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]  
     else:  
         for i in range(len(c)):  
             w1=c[i]*vrot[i]-s[i]*vrot[i+1]  
         w2=s[i]*vrot[i]+c[i]*vrot[i+1]  
             vrot[i:i+2]=w1,w2  
     return vrot  
577    
578  ##############################################          @param x: a value
579  def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):          @return: norm of argument x
580  ################################################          @rtype: C{float}
581     m=iter_restart          @note: by default C{sqrt(self.bilinearform(x,x)} is returned.
582     iter=0          """
583     xc=x          s=self.bilinearform(x,x)
584     while True:          if s<0: raise NegativeNorm,"negative norm."
585        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max          return math.sqrt(s)
       xc,stopped=GMRESm(b*1, Aprod, Msolve, bilinearform, stoppingcriterium, x=xc*1, iter_max=iter_max-iter, iter_restart=m)  
       if stopped: break  
       iter+=iter_restart      
    return xc  
586    
587  def GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):      def eval(self,x):
588     iter=0          """
589     r=Msolve(b)          Returns the value F of a given C{x}.
    r_dot_r = bilinearform(r, r)  
    if r_dot_r<0: raise NegativeNorm,"negative norm."  
    norm_b=math.sqrt(r_dot_r)  
590    
591     if x==None:          @param x: value for which the defect C{F} is evaluated
592        x=0*b          @return: value of the defect at C{x}
593     else:          """
594        r=Msolve(b-Aprod(x))          return 0
       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)  
    v=[]  
595    
596     rho=math.sqrt(r_dot_r)      def __call__(self,x):
597              return self.eval(x)
    v.append(r/rho)  
    g[0]=rho  
    while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):  
598    
599      if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max      def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):
600            """
601            Sets the relative length of the increment used to approximate the
602            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
606      p=Msolve(Aprod(v[iter]))          @type inc: positive C{float}
607            """
608            if inc<=0: raise ValueError,"positive increment required."
609            self.__inc=inc
610    
611      v.append(p)      def getDerivativeIncrementLength(self):
612            """
613            Returns the relative increment length used to approximate the
614            derivative of the defect.
615            @return: value of the defect at C{x}
616            @rtype: positive C{float}
617            """
618            return self.__inc
619    
620      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))        def derivative(self, F0, x0, v, v_is_normalised=True):
621            """
622            Returns the directional derivative at C{x0} in the direction of C{v}.
623    
624  # Modified Gram-Schmidt          @param F0: value of this defect at x0
625      for j in range(iter+1):          @param x0: value at which derivative is calculated
626        h[j][iter]=bilinearform(v[j],v[iter+1])            @param v: direction
627        v[iter+1]-=h[j][iter]*v[j]          @param v_is_normalised: True to indicate that C{v} is nomalized
628                                          (self.norm(v)=0)
629      h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))          @return: derivative of this defect at x0 in the direction of C{v}
630      v_norm2=h[iter+1][iter]          @note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is
631                   used but this method maybe overwritten to use exact evaluation.
632            """
633            normx=self.norm(x0)
634            if normx>0:
635                 epsnew = self.getDerivativeIncrementLength() * normx
636            else:
637                 epsnew = self.getDerivativeIncrementLength()
638            if not v_is_normalised:
639                normv=self.norm(v)
640                if normv<=0:
641                   return F0*0
642                else:
643                   epsnew /= normv
644            F1=self.eval(x0 + epsnew * v)
645            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):
649       """
650       Solves a non-linear problem M{F(x)=0} for unknown M{x} using the stopping
651       criterion:
652    
653  # Reorthogonalize if needed     M{norm(F(x) <= atol + rtol * norm(F(x0)}
     if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)  
      for j in range(iter+1):  
         hr=bilinearform(v[j],v[iter+1])  
             h[j][iter]=h[j][iter]+hr #vhat  
             v[iter+1] -= hr*v[j]  
654    
655       v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))       where M{x0} is the initial guess.
      h[iter+1][iter]=v_norm2  
656    
657  #   watch out for happy breakdown     @param defect: object defining the function M{F}. C{defect.norm} defines the
658          if not v_norm2 == 0:                    M{norm} used in the stopping criterion.
659           v[iter+1]=v[iter+1]/h[iter+1][iter]     @type defect: L{Defect}
660       @param x: initial guess for the solution, C{x} is altered.
661       @type x: any object type allowing basic operations such as
662                C{numpy.ndarray}, L{Data}
663       @param iter_max: maximum number of iteration steps
664       @type iter_max: positive C{int}
665       @param sub_iter_max: maximum number of inner iteration steps
666       @type sub_iter_max: positive C{int}
667       @param atol: absolute tolerance for the solution
668       @type atol: positive C{float}
669       @param rtol: relative tolerance for the solution
670       @type rtol: positive C{float}
671       @param gamma: tolerance safety factor for inner iteration
672       @type gamma: positive C{float}, less than 1
673       @param sub_tol_max: upper bound for inner tolerance
674       @type sub_tol_max: positive C{float}, less than 1
675       @return: an approximation of the solution with the desired accuracy
676       @rtype: same type as the initial guess
677       """
678       lmaxit=iter_max
679       if atol<0: raise ValueError,"atol needs to be non-negative."
680       if rtol<0: raise ValueError,"rtol needs to be non-negative."
681       if rtol+atol<=0: raise ValueError,"rtol or atol needs to be non-negative."
682       if gamma<=0 or gamma>=1: raise ValueError,"tolerance safety factor for inner iteration (gamma =%s) needs to be positive and less than 1."%gamma
683       if sub_tol_max<=0 or sub_tol_max>=1: raise ValueError,"upper bound for inner tolerance for inner iteration (sub_tol_max =%s) needs to be positive and less than 1."%sub_tol_max
684    
685       F=defect(x)
686       fnrm=defect.norm(F)
687       stop_tol=atol + rtol*fnrm
688       sub_tol=sub_tol_max
689       if verbose: print "NewtonGMRES: initial residual = %e."%fnrm
690       if verbose: print "             tolerance = %e."%sub_tol
691       iter=1
692       #
693       # main iteration loop
694       #
695       while not fnrm<=stop_tol:
696                if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
697                #
698            #   adjust sub_tol_
699            #
700                if iter > 1:
701               rat=fnrm/fnrmo
702                   sub_tol_old=sub_tol
703               sub_tol=gamma*rat**2
704               if gamma*sub_tol_old**2 > .1: sub_tol=max(sub_tol,gamma*sub_tol_old**2)
705               sub_tol=max(min(sub_tol,sub_tol_max), .5*stop_tol/fnrm)
706            #
707            # calculate newton increment xc
708                #     if iter_max in __FDGMRES is reached MaxIterReached is thrown
709                #     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
711                #
712                #
713                if verbose: print "             subiteration (GMRES) is called with relative tolerance %e."%sub_tol
714                try:
715                   xc, sub_iter=__FDGMRES(F, defect, x, sub_tol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)
716                except MaxIterReached:
717                   raise MaxIterReached,"maximum number of %s steps reached."%iter_max
718                if sub_iter<0:
719                   iter+=sub_iter_max
720                else:
721                   iter+=sub_iter
722                # ====
723            x+=xc
724                F=defect(x)
725            iter+=1
726                fnrmo, fnrm=fnrm, defect.norm(F)
727                if verbose: print "             step %s: residual %e."%(iter,fnrm)
728       if verbose: print "NewtonGMRES: completed after %s steps."%iter
729       return x
730    
731  #   Form and store the information for the new Givens rotation  def __givapp(c,s,vin):
732      if iter > 0 :      """
733          hhat=numarray.zeros(iter+1,numarray.Float64)      Applies a sequence of Givens rotations (c,s) recursively to the vector
734          for i in range(iter+1) : hhat[i]=h[i][iter]      C{vin}
         hhat=givapp(c[0:iter],s[0:iter],hhat);  
             for i in range(iter+1) : h[i][iter]=hhat[i]  
735    
736      mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])      @warning: C{vin} is altered.
737      if mu!=0 :      """
738          c[iter]=h[iter][iter]/mu      vrot=vin
739          s[iter]=-h[iter+1][iter]/mu      if isinstance(c,float):
740          h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]          vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
741          h[iter+1][iter]=0.0      else:
742          g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])          for i in range(len(c)):
743                w1=c[i]*vrot[i]-s[i]*vrot[i+1]
744            w2=s[i]*vrot[i]+c[i]*vrot[i+1]
745                vrot[i]=w1
746                vrot[i+1]=w2
747        return vrot
748    
749  # Update the residual norm  def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
750          rho=abs(g[iter+1])     h=numpy.zeros((iter_restart,iter_restart),numpy.float64)
751      iter+=1     c=numpy.zeros(iter_restart,numpy.float64)
752       s=numpy.zeros(iter_restart,numpy.float64)
753       g=numpy.zeros(iter_restart,numpy.float64)
754       v=[]
755    
756  # At this point either iter > iter_max or rho < tol.     rho=defect.norm(F0)
757  # It's time to compute x and leave.             if rho<=0.: return x0*0
758    
759       v.append(-F0/rho)
760       g[0]=rho
761       iter=0
762       while rho > atol and iter<iter_restart-1:
763            if iter  >= iter_max:
764                raise MaxIterReached,"maximum number of %s steps reached."%iter_max
765    
766            p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
767            v.append(p)
768    
769            v_norm1=defect.norm(v[iter+1])
770    
771            # 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            h[iter+1,iter]=defect.norm(v[iter+1])
777            v_norm2=h[iter+1,iter]
778    
779            # Reorthogonalize if needed
780            if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
781                for j in range(iter+1):
782                    hr=defect.bilinearform(v[j],v[iter+1])
783                    h[j,iter]=h[j,iter]+hr
784                    v[iter+1] -= hr*v[j]
785    
786                v_norm2=defect.norm(v[iter+1])
787                h[iter+1,iter]=v_norm2
788            #   watch out for happy breakdown
789            if not v_norm2 == 0:
790                v[iter+1]=v[iter+1]/h[iter+1,iter]
791    
792            #   Form and store the information for the new Givens rotation
793            if iter > 0 :
794                hhat=numpy.zeros(iter+1,numpy.float64)
795                for i in range(iter+1) : hhat[i]=h[i,iter]
796                hhat=__givapp(c[0:iter],s[0:iter],hhat);
797                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])
800    
801            if mu!=0 :
802                c[iter]=h[iter,iter]/mu
803                s[iter]=-h[iter+1,iter]/mu
804                h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
805                h[iter+1,iter]=0.0
806                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
811            rho=abs(g[iter+1])
812            iter+=1
813    
814     if iter > 0 :     # At this point either iter > iter_max or rho < tol.
815       y=numarray.zeros(iter,numarray.Float64)         # It's time to compute x and leave.
816       y[iter-1] = g[iter-1] / h[iter-1][iter-1]     if iter > 0 :
817       if iter > 1 :         y=numpy.zeros(iter,numpy.float64)
818          i=iter-2         y[iter-1] = g[iter-1] / h[iter-1,iter-1]
819         if iter > 1 :
820            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 : xhat=v[0]     else :
828              xhat=v[0] * 0
    x += xhat  
    if iter<iter_restart-1:  
       stopped=True  
    else:  
       stopped=False  
   
    return x,stopped  
829    
830       if iter<iter_restart-1:
831          stopped=iter
832       else:
833          stopped=-1
834    
835  ######################################################     return xhat,stopped
 def dirder(x, w, bilinearform, Aprod, Msolve, f0, b ):  
 ######################################################  
836    
837  # DIRDER estimates the directional derivative of a function F.  def GMRES(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100, iter_restart=20, verbose=False,P_R=None):
838       """
839       Solver for
840    
841       M{Ax=b}
842    
843  # Hardwired difference increment.     with a general operator A (more details required!).
844  #     It uses the generalized minimum residual method (GMRES).
   epsnew = 1.0e-07  
 #  
 #  Scale the step.  
 #  
   norm_w=math.sqrt(bilinearform(w,w))  
   if norm_w== 0.0:  
     return x/x  
845    
846    epsnew = epsnew / norm_w     The iteration is terminated if
847    
848    if norm_w > 0.0:     M{|r| <= atol+rtol*|r0|}
     epsnew = epsnew * math.sqrt(bilinearform(x,x))  
 #  
 #  DEL and F1 could share the same space if storage  
 #  is more important than clarity.  
 #  
849    
850    DEL = x + epsnew * w     where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
   f1 = -Msolve(Aprod(DEL))  
   z = ( f1 - f0 ) / epsnew  
   return z  
   
 ######################################################  
 def FDGMRES(f0, Aprod, Msolve, bilinearform, stoppingcriterium, xc=None, x=None, iter_max=100, iter_restart=20,TOL=None):  
 ######################################################  
    b=-f0  
    b_dot_b = bilinearform(b, b)  
    if b_dot_b<0: raise NegativeNorm,"negative norm."  
    norm_b=math.sqrt(b_dot_b)  
851    
852     r=b     M{|r| = sqrt( bilinearform(r,r))}
853    
854     if x==None:     @param r: initial residual M{r=b-Ax}. C{r} is altered.
855        x=0*f0     @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
881       restarted=False
882       iter=0
883       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:     else:
889        r=-dirder(xc,x,bilinearform,Aprod,Msolve,f0,b)-f0          atol2=atol
890                if verbose: print "GMRES: absolute tolerance = %e"%atol2
891     r_dot_r = bilinearform(r, r)     if atol2<=0:
892     if r_dot_r<0: raise NegativeNorm,"negative norm."        raise ValueError,"Non-positive tolarance."
893      
894     h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)     while True:
895     c=numarray.zeros(iter_restart,numarray.Float64)        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
896     s=numarray.zeros(iter_restart,numarray.Float64)        if restarted:
897     g=numarray.zeros(iter_restart,numarray.Float64)           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
904          if verbose: print "GMRES: restart."
905          restarted=True
906       if verbose: print "GMRES: tolerance has been reached."
907       return x
908    
909    def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False, P_R=None):
910       iter=0
911    
912       h=numpy.zeros((iter_restart+1,iter_restart),numpy.float64)
913       c=numpy.zeros(iter_restart,numpy.float64)
914       s=numpy.zeros(iter_restart,numpy.float64)
915       g=numpy.zeros(iter_restart+1,numpy.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
    iter=0  
924    
925     while not (stoppingcriterium(rho,norm_b,solver="FDGMRES",TOL=TOL) 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                if P_R!=None:
931          p=dirder(xc, v[iter], bilinearform,Aprod,Msolve,f0,b)              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  # Modified Gram-Schmidt
939      for j in range(iter+1):      for j in range(iter+1):
940        h[j][iter]=bilinearform(v[j],v[iter+1])          h[j,iter]=bilinearform(v[j],v[iter+1])
941        v[iter+1]+=(-1.)*h[j][iter]*v[j]        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]  
942    
943        h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
944        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 #vhat              h[j,iter]=h[j,iter]+hr
951              v[iter+1] +=(-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 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=[]      mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
         for i in range(iter+1) : hhat.append(h[i][iter])  
         hhat=givapp(c[0:iter],s[0:iter],hhat);  
             for i in range(iter+1) : h[i][iter]=hhat[i]  
963    
     mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])  
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     if iter > 0 :  # At this point either iter > iter_max or rho < tol.
979       y=numarray.zeros(iter,numarray.Float64)      # It's time to compute x and leave.
980       y[iter-1] = g[iter-1] / h[iter-1][iter-1]  
981       if iter > 1 :       if verbose: print "GMRES: iteration stopped after %s step."%iter
982          i=iter-2       if iter > 0 :
983         y=numpy.zeros(iter,numpy.float64)
984         y[iter-1] = g[iter-1] / h[iter-1,iter-1]
985         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 891  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 918  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 935  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 943  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 953  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 961  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 978  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(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1185  def NewtonGMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium,x=None, iter_max=100,iter_restart=20,atol=1.e-2,rtol=1.e-4):    """
1186  #####################################    Solver for
     gamma=.9  
     lmaxit=100  
     etamax=.5  
   
     n = 1 #len(x)  
     iter=0  
       
     # evaluate f at the initial iterate  
     # compute the stop tolerance  
     #  
     r=b  
     if x==None:  
       x=0*b  
     else:  
       r += (-1)*Aprod(x)  
   
     f0=-Msolve(r)  
     fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)  
     fnrmo=1  
     stop_tol=atol + rtol*fnrm  
     #  
     # main iteration loop  
     #  
     while not stoppingcriterium(fnrm*1,stop_tol,'NewtonGMRES',TOL=1.):  
   
             if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max  
         #  
         # keep track of the ratio (rat = fnrm/frnmo)  
         # of successive residual norms and  
         # the iteration counter (iter)  
         #  
         #rat=fnrm/fnrmo  
         fnrmo=fnrm  
         iter+=1  
         #  
             # compute the step using a GMRES(m) routine especially designed for this purpose  
         #  
             initer=0  
             while True:  
                xc,stopped=FDGMRES(f0*1, Aprod, Msolve, bilinearform, stoppingcriterium, xc=x*1, iter_max=lmaxit-initer, iter_restart=iter_restart, TOL=etamax)  
                if stopped: break  
                initer+=iter_restart  
         x+=xc  
         f0=-Msolve(Aprod(x))  
         fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)  
         rat=fnrm/fnrmo  
   
   
     #   adjust eta  
     #  
         if etamax > 0:  
             etaold=etamax  
             etanew=gamma*rat*rat  
             if gamma*etaold*etaold > .1 :  
                 etanew=max(etanew,gamma*etaold*etaold)  
               
             etamax=min(etanew,etamax)  
             etamax=max(etamax,.5*stop_tol/fnrm)  
   
     return x  
1187    
1188  def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):    M{Ax=b}
1189    
1190  # TFQMR solver for linear systems    with a general operator A (more details required!).
1191  #    It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).
 #  
 # 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  
1192    
1193    r=Msolve(r)    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 1109  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 1127  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 1155  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 1186  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 1196  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 1204  class ArithmeticTuple(object): Line 1338  class ArithmeticTuple(object):
1338         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1339         """         """
1340         out=[]         out=[]
1341         for i in range(len(self)):         try:
1342             out.append(self[i]*other)             l=len(other)
1343               if l!=len(self):
1344                   raise ValueError,"length of arguments don't match."
1345               for i in range(l): out.append(self[i]*other[i])
1346           except TypeError:
1347               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 1218  class ArithmeticTuple(object): Line 1357  class ArithmeticTuple(object):
1357         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1358         """         """
1359         out=[]         out=[]
1360         for i in range(len(self)):         try:
1361             out.append(other*self[i])             l=len(other)
1362               if l!=len(self):
1363                   raise ValueError,"length of arguments don't match."
1364               for i in range(l): out.append(other[i]*self[i])
1365           except TypeError:
1366               for i in range(len(self)): out.append(other*self[i])
1367         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1368    
 #########################  
 # Added by Artak  
 #########################  
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}
1375         @return: itemwise self/other         @return: itemwise self/other
1376         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1377         """         """
1378         out=[]         return self*(1/other)
        for i in range(len(self)):  
            out.append(self[i]*(1./other))  
        return ArithmeticTuple(*tuple(out))  
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 1249  class ArithmeticTuple(object): Line 1387  class ArithmeticTuple(object):
1387         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1388         """         """
1389         out=[]         out=[]
1390         for i in range(len(self)):         try:
1391             out.append(other/self[i])             l=len(other)
1392               if l!=len(self):
1393                   raise ValueError,"length of arguments don't match."
1394               for i in range(l): out.append(other[i]/self[i])
1395           except TypeError:
1396               for i in range(len(self)): out.append(other/self[i])
1397         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
     
 ##########################################33  
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         if not isinstance(other,float):         out=[]
1420          if len(self) != len(other):         try:
1421             raise ValueError,"tuple length must match."             l=len(other)
1422          for i in range(len(self)):             if l!=len(self):
1423             self.__items[i]+=other[i]                 raise ValueError,"length of arguments don't match."
1424         else:             for i in range(l): out.append(self[i]+other[i])
1425          for i in range(len(self)):         except TypeError:
1426             self.__items[i]+=other             for i in range(len(self)): out.append(self[i]+other)
1427           return ArithmeticTuple(*tuple(out))
        return self  
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}
1435           """
1436           out=[]
1437           try:
1438               l=len(other)
1439               if l!=len(self):
1440                   raise ValueError,"length of arguments don't match."
1441               for i in range(l): out.append(self[i]-other[i])
1442           except TypeError:
1443               for i in range(len(self)): out.append(self[i]-other)
1444           return ArithmeticTuple(*tuple(out))
1445    
1446       def __isub__(self,other):
1447           """
1448           Inplace subtraction of C{other} from self.
1449    
1450           @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 1301  class ArithmeticTuple(object): Line 1458  class ArithmeticTuple(object):
1458    
1459     def __neg__(self):     def __neg__(self):
1460         """         """
1461         negate         Negates values.
   
1462         """         """
1463           out=[]
1464         for i in range(len(self)):         for i in range(len(self)):
1465             self.__items[i]=-self.__items[i]             out.append(-self[i])
1466         return self         return ArithmeticTuple(*tuple(out))
1467    
1468    
1469  class HomogeneousSaddlePointProblem(object):  class HomogeneousSaddlePointProblem(object):
1470        """        """
1471        This provides a framwork for solving 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
       def B(self,v):  
          """  
          returns Bv (overwrite)  
          @rtype: equal to the type of p  
1496    
1497           @note: boundary conditions on p should be zero!        def inner_pBv(self,p,Bv):
1498           """           """
1499           pass           Returns inner product of element p and Bv (overwrite).
1500    
1501        def inner(self,p0,p1):           @param p: a pressure increment
1502           """           @param v: a residual
1503           returns inner product of two element p0 and p1  (overwrite)           @return: inner product of element p and Bv
           
          @type p0: equal to the type of p  
          @type p1: equal to the type of p  
1504           @rtype: C{float}           @rtype: C{float}
1505             @note: used if PCG is applied.
          @rtype: equal to the type of p  
1506           """           """
1507           pass           raise NotImplementedError,"no inner product for p and Bv implemented."
1508    
1509        def solve_A(self,u,p):        def inner_p(self,p0,p1):
1510           """           """
1511           solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)           Returns inner product of p0 and p1 (overwrite).
1512    
1513           @rtype: equal to the type of v           @param p0: a pressure
1514           @note: boundary conditions on v should be zero!           @param p1: a pressure
1515             @return: inner product of p0 and p1
1516             @rtype: C{float}
1517           """           """
1518           pass           raise NotImplementedError,"no inner product for p implemented."
1519      
1520        def solve_prec(self,p):        def norm_v(self,v):
1521           """           """
1522           provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)           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):
       def stoppingcriterium(self,Bv,v,p):  
1530           """           """
1531           returns a True if iteration is terminated. (overwrite)           return the value for v for a given p (overwrite)
1532    
1533           @rtype: C{bool}           @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           pass           raise NotImplementedError,"no v calculation implemented."
               
       def __inner(self,p,r):  
          return self.inner(p,r[1])  
   
       def __inner_p(self,p1,p2):  
          return self.inner(p1,p2)  
         
       def __inner_a(self,a1,a2):  
          return self.inner_a(a1,a2)  
   
       def __inner_a1(self,a1,a2):  
          return self.inner(a1[1],a2[1])  
1538    
1539        def __stoppingcriterium(self,norm_r,r,p):          
1540            return self.stoppingcriterium(r[1],r[0],p)        def Bv(self,v):
1541            """
1542        def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):          Returns Bv (overwrite).
           return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)  
1543    
1544        def setTolerance(self,tolerance=1.e-8):          @rtype: equal to the type of p
1545                self.__tol=tolerance          @note: boundary conditions on p should be zero!
1546        def getTolerance(self):          """
1547                return self.__tol          raise NotImplementedError, "no operator B implemented."
       def setToleranceReductionFactor(self,reduction=0.01):  
               self.__reduction=reduction  
       def getSubProblemTolerance(self):  
               return self.__reduction*self.getTolerance()  
1548    
1549        def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):        def norm_Bv(self,Bv):
1550                """          """
1551                solves the saddle point problem using initial guesses v and p.          Returns the norm of Bv (overwrite).
   
               @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=='NewtonGMRES':      
                 if self.verbose: print "enter NewtonGMRES method (iter_max=%s)"%max_iter  
                 p=NewtonGMRES(Bz,self.__Aprod_Newton,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,atol=0,rtol=self.getTolerance())  
                 # 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_prec(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)  
1552    
1553                print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)          @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            return u,p        def solve_AinvBt(self,p):
1559             """
1560             Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
1561             (overwrite).
1562    
1563        def __Msolve(self,r):           @param p: a pressure increment
1564            return self.solve_prec(r[1])           @return: the solution of M{Av=B^*p}
1565             @note: boundary conditions on v should be zero!
1566             """
1567             raise NotImplementedError,"no operator A implemented."
1568    
1569        def __Msolve2(self,r):        def solve_prec(self,Bv):
1570            return self.solve_prec(r*1)           """
1571             Provides a preconditioner for M{BA^{-1}B^*} applied to Bv with accuracy
1572             L{self.getSubProblemTolerance()} (overwrite).
1573    
1574        def __Mempty(self,r):           @rtype: equal to the type of p
1575            return r           @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
1584          #=============================================================
1585          def __Aprod_PCG(self,p):
1586              dv=self.solve_AinvBt(p)
1587              return ArithmeticTuple(dv,self.Bv(dv))
1588    
1589          def __inner_PCG(self,p,r):
1590             return self.inner_pBv(p,r[1])
1591    
1592        def __Aprod(self,p):        def __Msolve_PCG(self,r):
1593            # return BA^-1B*p            return self.solve_prec(r[1])
1594            #solve Av =B^*p as Av =f-Az-B^*(-p)        #=============================================================
1595            v=self.solve_A(self.__z,-p)        def __Aprod_GMRES(self,p):
1596            return ArithmeticTuple(v, self.B(v))            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             Solves the saddle point problem using initial guesses v and p.
1621    
1622        def __Anext(self,x):           @param v: initial guess for velocity
1623            # return next v,p           @param p: initial guess for pressure
1624            #solve Av =-B^*p as Av =f-Az-B^*p           @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             self.verbose=verbose
1640             rtol=self.getTolerance()
1641             atol=self.getAbsoluteTolerance()
1642         if self.adaptSubTolerance(): self.setSubProblemTolerance()
1643             correction_step=0
1644             converged=False
1645             while not converged:
1646                  # calculate velocity for current pressure:
1647                  v=self.getV(p,v)
1648                  Bv=self.Bv(v)
1649                  norm_v=self.norm_v(v)
1650                  norm_Bv=self.norm_Bv(Bv)
1651                  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        pc=x[1]        #========================================================================
1677            v=self.solve_A(self.__z,-pc)        def setTolerance(self,tolerance=1.e-4):
1678        p=self.solve_prec(self.B(v))           """
1679             Sets the relative tolerance for (v,p).
1680    
1681            return ArithmeticTuple(v,p)           @param tolerance: tolerance to be used
1682             @type tolerance: non-negative C{float}
1683             """
1684             if tolerance<0:
1685                 raise ValueError,"tolerance must be positive."
1686             self.__rtol=tolerance
1687    
1688          def getTolerance(self):
1689             """
1690             Returns the relative tolerance.
1691    
1692        def __Aprod2(self,p):           @return: relative tolerance
1693            # return BA^-1B*p           @rtype: C{float}
1694            #solve Av =B^*p as Av =f-Az-B^*(-p)           """
1695        v=self.solve_A(self.__z,-p)           return self.__rtol
           return self.B(v)  
1696    
1697        def __Aprod_Newton(self,p):        def setAbsoluteTolerance(self,tolerance=0.):
1698            # return BA^-1B*p - Bz           """
1699            #solve Av =-B^*p as Av =f-Az-B^*p           Sets the absolute tolerance.
       v=self.solve_A(self.__z,-p)  
           return self.B(v-self.__z)  
1700    
1701        def __Aprod_Newton2(self,x):           @param tolerance: tolerance to be used
1702            # return BA^-1B*p - Bz           @type tolerance: non-negative C{float}
1703            #solve Av =-B^*p as Av =f-Az-B^*p           """
1704            pc=x[1]           if tolerance<0:
1705        v=self.solve_A(self.__z,-pc)               raise ValueError,"tolerance must be non-negative."
1706            p=self.solve_prec(self.B(v-self.__z))           self.__atol=tolerance
           return ArithmeticTuple(v,p)  
1707    
1708  class SaddlePointProblem(object):        def getAbsoluteTolerance(self):
1709     """           """
1710     This implements a solver for a saddlepoint problem           Returns the absolute tolerance.
1711    
1712     M{f(u,p)=0}           @return: absolute tolerance
1713     M{g(u)=0}           @rtype: C{float}
1714             """
1715             return self.__atol
1716    
1717     for u and p. The problem is solved with an inexact Uszawa scheme for p:        def getSubProblemTolerance(self):
1718             """
1719             Sets the relative tolerance to solve the subproblem(s).
1720    
1721     M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}           @param rtol: relative tolerence
1722     M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}           @type rtol: positive C{float}
1723             """
1724             return max(200.*util.EPSILON,self.getTolerance()**2)
1725    
1726     where Q_f is an approximation of the Jacobiean A_f of f with respect to u  and Q_f is an approximation of  def MaskFromBoundaryTag(domain,*tags):
    A_g A_f^{-1} A_g with A_g is the jacobiean of g with respect to p. As a the construction of a 'proper'  
    Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays  
    in fact the role of a preconditioner.  
1727     """     """
1728     def __init__(self,verbose=False,*args):     Creates a mask on the Solution(domain) function space where the value is
1729         """     one for samples that touch the boundary 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  
        """  
        if not isinstance(verbose,bool):  
             raise TypeError("verbose needs to be of type bool.")  
        self.__verbose=verbose  
        self.relaxation=1.  
   
    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  
   
        A_f du = f(u,p)  
1730    
1731         with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.     Usage: m=MaskFromBoundaryTag(domain, "left", "right")
   
        @param u: current approximation of u  
        @type u: L{escript.Data}  
        @param p: current approximation of p  
        @type p: L{escript.Data}  
        @param tol: tolerance expected for du  
        @type tol: C{float}  
        @return: increment du  
        @rtype: L{escript.Data}  
        @note: this method has to be overwritten by a particular saddle point problem  
        """  
        pass  
   
    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  
1732    
1733     subiter_max=3     @param domain: domain to be used
1734     def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):     @type domain: L{escript.Domain}
1735          """     @param tags: boundary tags
1736          runs the solver     @type tags: C{str}
1737       @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
1740       """
1741       pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1742       d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
1743       for t in tags: d.setTaggedValue(t,1.)
1744       pde.setValue(y=d)
1745       return util.whereNonZero(pde.getRightHandSide())
1746    
1747          @param u0: initial guess for C{u}  def MaskFromTag(domain,*tags):
         @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  
 #   def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):  
 #      tol=1.e-2  
 #      iter=0  
 #      converged=False  
 #      u=u0*1.  
 #      p=p0*1.  
 #      while not converged and iter<iter_max:  
 #          du=self.solve_f(u,p,tol)  
 #          u-=du  
 #          norm_du=util.Lsup(du)  
 #          norm_u=util.Lsup(u)  
 #          
 #          dp=self.relaxation*self.solve_g(u,tol)  
 #          p+=dp  
 #          norm_dp=util.sqrt(self.inner(dp,dp))  
 #          norm_p=util.sqrt(self.inner(p,p))  
 #          print iter,"-th step rel. errror u,p= (%s,%s) (%s,%s)(%s,%s)"%(norm_du,norm_dp,norm_du/norm_u,norm_dp/norm_p,norm_u,norm_p)  
 #          iter+=1  
 #  
 #          converged = (norm_du <= tolerance*norm_u) and  (norm_dp <= tolerance*norm_p)  
 #      if converged:  
 #          print "convergence after %s steps."%iter  
 #      else:  
 #          raise ArithmeticError("no convergence after %s steps."%iter)  
 #  
 #      return u,p  
             
 def MaskFromBoundaryTag(function_space,*tags):  
1748     """     """
1749     create a mask on the given function space which one for samples     Creates a mask on the Solution(domain) function space where the value is
1750     that touch the boundary tagged by tags.     one for samples that touch regions tagged by tags.
1751    
1752     usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")     Usage: m=MaskFromTag(domain, "ham")
1753    
1754     @param function_space: a given function space     @param domain: domain to be used
1755     @type function_space: L{escript.FunctionSpace}     @type domain: L{escript.Domain}
1756     @param tags: boundray tags     @param tags: boundary tags
1757     @type tags: C{str}     @type tags: C{str}
1758     @return: a mask which marks samples used by C{function_space} that are touching the     @return: a mask which marks samples that are touching the boundary tagged
1759              boundary tagged by any of the given tags.              by any of the given tags
1760     @rtype: L{escript.Data} of rank 0     @rtype: L{escript.Data} of rank 0
1761     """     """
1762     pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1763     d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))     d=escript.Scalar(0.,escript.Function(domain))
1764     for t in tags: d.setTaggedValue(t,1.)     for t in tags: d.setTaggedValue(t,1.)
1765     pde.setValue(y=d)     pde.setValue(Y=d)
1766     out=util.whereNonZero(pde.getRightHandSide())     return util.whereNonZero(pde.getRightHandSide())
    if out.getFunctionSpace() == function_space:  
       return out  
    else:  
       return util.whereNonZero(util.interpolate(out,function_space))  
   
1767    
1768    

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

  ViewVC Help
Powered by ViewVC 1.1.26