/[escript]/trunk/escript/py_src/pdetools.py
ViewVC logotype

Diff of /trunk/escript/py_src/pdetools.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

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

  ViewVC Help
Powered by ViewVC 1.1.26