/[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 2455 by jfenwick, Wed Jun 3 03:29:07 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.setSolverMethod(linearPDEs.LinearPDE.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.)
# Line 151  class Projector: Line 153  class Projector:
153    
154    def __call__(self, input_data):    def __call__(self, input_data):
155      """      """
156      Projects input_data onto a continuous function      Projects C{input_data} onto a continuous function.
157    
158      @param input_data: The input_data to be projected.      @param input_data: the data to be projected
159      """      """
160      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
161      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 188  class Projector:
188    
189  class NoPDE:  class NoPDE:
190       """       """
191       solves the following problem for u:       Solves the following problem for u:
192    
193       M{kronecker[i,j]*D[j]*u[j]=Y[i]}       M{kronecker[i,j]*D[j]*u[j]=Y[i]}
194    
195       with constraint       with constraint
196    
197       M{u[j]=r[j]}  where M{q[j]>0}       M{u[j]=r[j]}  where M{q[j]>0}
198    
199       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.
200    
201       In the case of scalars this takes the form       In the case of scalars this takes the form
202    
203       M{D*u=Y}       M{D*u=Y}
204    
205       with constraint       with constraint
206    
207       M{u=r}  where M{q>0}       M{u=r} where M{q>0}
208    
209       where D, Y, r and q are given scalar functions.       where M{D}, M{Y}, M{r} and M{q} are given scalar functions.
210    
211       The constraint is overwriting any other condition.       The constraint overwrites any other condition.
212    
213       @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
214              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
215              thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.              in L{Solution} or L{ReducedSolution}.
216       """       """
217         # The whole thing is a bit strange and I blame Rob Woodcock (CSIRO) for
218         # this.
219       def __init__(self,domain,D=None,Y=None,q=None,r=None):       def __init__(self,domain,D=None,Y=None,q=None,r=None):
220           """           """
221           initialize the problem           Initializes the problem.
222    
223           @param domain: domain of the PDE.           @param domain: domain of the PDE
224           @type domain: L{Domain}           @type domain: L{Domain}
225           @param D: coefficient of the solution.           @param D: coefficient of the solution
226           @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}           @type D: C{float}, C{int}, C{numpy.ndarray}, L{Data}
227           @param Y: right hand side           @param Y: right hand side
228           @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}           @type Y: C{float}, C{int}, C{numpy.ndarray}, L{Data}
229           @param q: location of constraints           @param q: location of constraints
230           @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}           @type q: C{float}, C{int}, C{numpy.ndarray}, L{Data}
231           @param r: value of solution at locations of constraints           @param r: value of solution at locations of constraints
232           @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}           @type r: C{float}, C{int}, C{numpy.ndarray}, L{Data}
233           """           """
234           self.__domain=domain           self.__domain=domain
235           self.__D=D           self.__D=D
# Line 234  class NoPDE: Line 238  class NoPDE:
238           self.__r=r           self.__r=r
239           self.__u=None           self.__u=None
240           self.__function_space=escript.Solution(self.__domain)           self.__function_space=escript.Solution(self.__domain)
241    
242       def setReducedOn(self):       def setReducedOn(self):
243           """           """
244           sets the L{FunctionSpace} of the solution to L{ReducedSolution}           Sets the L{FunctionSpace} of the solution to L{ReducedSolution}.
245           """           """
246           self.__function_space=escript.ReducedSolution(self.__domain)           self.__function_space=escript.ReducedSolution(self.__domain)
247           self.__u=None           self.__u=None
248    
249       def setReducedOff(self):       def setReducedOff(self):
250           """           """
251           sets the L{FunctionSpace} of the solution to L{Solution}           Sets the L{FunctionSpace} of the solution to L{Solution}.
252           """           """
253           self.__function_space=escript.Solution(self.__domain)           self.__function_space=escript.Solution(self.__domain)
254           self.__u=None           self.__u=None
255            
256       def setValue(self,D=None,Y=None,q=None,r=None):       def setValue(self,D=None,Y=None,q=None,r=None):
257           """           """
258           assigns values to the parameters.           Assigns values to the parameters.
259    
260           @param D: coefficient of the solution.           @param D: coefficient of the solution
261           @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}           @type D: C{float}, C{int}, C{numpy.ndarray}, L{Data}
262           @param Y: right hand side           @param Y: right hand side
263           @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}           @type Y: C{float}, C{int}, C{numpy.ndarray}, L{Data}
264           @param q: location of constraints           @param q: location of constraints
265           @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}           @type q: C{float}, C{int}, C{numpy.ndarray}, L{Data}
266           @param r: value of solution at locations of constraints           @param r: value of solution at locations of constraints
267           @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}           @type r: C{float}, C{int}, C{numpy.ndarray}, L{Data}
268           """           """
269           if not D==None:           if not D==None:
270              self.__D=D              self.__D=D
# Line 276  class NoPDE: Line 281  class NoPDE:
281    
282       def getSolution(self):       def getSolution(self):
283           """           """
284           returns the solution           Returns the solution.
285            
286           @return: the solution of the problem           @return: the solution of the problem
287           @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.           @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or
288                     L{ReducedSolution}
289           """           """
290           if self.__u==None:           if self.__u==None:
291              if self.__D==None:              if self.__D==None:
# Line 296  class NoPDE: Line 302  class NoPDE:
302                  self.__u*=(1.-q)                  self.__u*=(1.-q)
303                  if not self.__r==None: self.__u+=q*self.__r                  if not self.__r==None: self.__u+=q*self.__r
304           return self.__u           return self.__u
305                
306  class Locator:  class Locator:
307       """       """
308       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
309       spatial coordinate x.         coordinate x.
310        
311       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
312       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.  
313       """       """
314    
315       def __init__(self,where,x=numarray.zeros((3,))):       def __init__(self,where,x=numpy.zeros((3,))):
316         """         """
317         Initializes a Locator to access values in Data objects on the Doamin         Initializes a Locator to access values in Data objects on the Doamin
318         or FunctionSpace where for the sample point which         or FunctionSpace for the sample point which is closest to the given
319         closest to the given point x.         point x.
320    
321         @param where: function space         @param where: function space
322         @type where: L{escript.FunctionSpace}         @type where: L{escript.FunctionSpace}
323         @param x: coefficient of the solution.         @param x: coefficient of the solution
324         @type x: L{numarray.NumArray} or C{list} of L{numarray.NumArray}         @type x: C{numpy.ndarray} or C{list} of C{numpy.ndarray}
325         """         """
326         if isinstance(where,escript.FunctionSpace):         if isinstance(where,escript.FunctionSpace):
327            self.__function_space=where            self.__function_space=where
# Line 350  class Locator: Line 355  class Locator:
355    
356       def getX(self):       def getX(self):
357          """          """
358      Returns the exact coordinates of the Locator.          Returns the exact coordinates of the Locator.
359      """          """
360          return self(self.getFunctionSpace().getX())          return self(self.getFunctionSpace().getX())
361    
362       def getFunctionSpace(self):       def getFunctionSpace(self):
363          """          """
364      Returns the function space of the Locator.          Returns the function space of the Locator.
365      """          """
366          return self.__function_space          return self.__function_space
367    
368       def getId(self,item=None):       def getId(self,item=None):
369          """          """
370      Returns the identifier of the location.          Returns the identifier of the location.
371      """          """
372          if item == None:          if item == None:
373             return self.__id             return self.__id
374          else:          else:
# Line 375  class Locator: Line 380  class Locator:
380    
381       def __call__(self,data):       def __call__(self,data):
382          """          """
383      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.
384      the object is returned.          """
     """  
385          return self.getValue(data)          return self.getValue(data)
386    
387       def getValue(self,data):       def getValue(self,data):
388          """          """
389      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}
390      otherwise the object is returned.          object otherwise the object is returned.
391      """          """
392          if isinstance(data,escript.Data):          if isinstance(data,escript.Data):
393             if data.getFunctionSpace()==self.getFunctionSpace():             if data.getFunctionSpace()==self.getFunctionSpace():
394               dat=data               dat=data
# Line 395  class Locator: Line 399  class Locator:
399             if isinstance(id,list):             if isinstance(id,list):
400                 out=[]                 out=[]
401                 for i in id:                 for i in id:
402                    o=data.getValueOfGlobalDataPoint(*i)                    o=numpy.array(data.getTupleForGlobalDataPoint(*i))
403                    if data.getRank()==0:                    if data.getRank()==0:
404                       out.append(o[0])                       out.append(o[0])
405                    else:                    else:
406                       out.append(o)                       out.append(o)
407                 return out                 return out
408             else:             else:
409               out=data.getValueOfGlobalDataPoint(*id)               out=numpy.array(data.getTupleForGlobalDataPoint(*id))
410               if data.getRank()==0:               if data.getRank()==0:
411                  return out[0]                  return out[0]
412               else:               else:
# Line 412  class Locator: Line 416  class Locator:
416    
417  class SolverSchemeException(Exception):  class SolverSchemeException(Exception):
418     """     """
419     exceptions thrown by solvers     This is a generic exception thrown by solvers.
420     """     """
421     pass     pass
422    
423  class IndefinitePreconditioner(SolverSchemeException):  class IndefinitePreconditioner(SolverSchemeException):
424     """     """
425     the preconditioner is not positive definite.     Exception thrown if the preconditioner is not positive definite.
426     """     """
427     pass     pass
428    
429  class MaxIterReached(SolverSchemeException):  class MaxIterReached(SolverSchemeException):
430     """     """
431     maxium number of iteration steps is reached.     Exception thrown if the maximum number of iteration steps is reached.
432     """     """
433     pass     pass
434  class IterationBreakDown(SolverSchemeException):  
435    class CorrectionFailed(SolverSchemeException):
436     """     """
437     iteration scheme econouters an incurable breakdown.     Exception thrown if no convergence has been achieved in the solution
438       correction scheme.
439     """     """
440     pass     pass
441  class NegativeNorm(SolverSchemeException):  
442    class IterationBreakDown(SolverSchemeException):
443     """     """
444     a norm calculation returns a negative norm.     Exception thrown if the iteration scheme encountered an incurable breakdown.
445     """     """
446     pass     pass
447    
448  class IterationHistory(object):  class NegativeNorm(SolverSchemeException):
449     """     """
450     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.  
451     """     """
452     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  
453    
454  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):
455     """     """
456     Solver for     Solver for
457    
458     M{Ax=b}     M{Ax=b}
459    
460     with a symmetric and positive definite operator A (more details required!).     with a symmetric and positive definite operator A (more details required!).
461     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
462       approximation of A.
463    
464     The iteration is terminated if the C{stoppingcriterium} function return C{True}.     The iteration is terminated if
465    
466     For details on the preconditioned conjugate gradient method see the book:     M{|r| <= atol+rtol*|r0|}
467    
468     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.  
469    
470     @param b: the right hand side of the liner system. C{b} is altered.     M{|r| = sqrt( bilinearform(Msolve(r),r))}
471     @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)  
472     @param Aprod: returns the value Ax     For details on the preconditioned conjugate gradient method see the book:
473     @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}.  
474     @param Msolve: solves Mx=r     I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
475     @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,
476  type like argument C{x}.     C. Romine, and H. van der Vorst}.
477     @param bilinearform: inner product C{<x,r>}  
478     @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.
479     @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)
480     @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.  
481     @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)
482     @param iter_max: maximum number of iteration steps.     @param Aprod: returns the value Ax
483       @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
484                    argument C{x}. The returned object needs to be of the same type
485                    like argument C{r}.
486       @param Msolve: solves Mx=r
487       @type Msolve: function C{Msolve(r)} where C{r} is of the same type like
488                     argument C{r}. The returned object needs to be of the same
489                     type like argument C{x}.
490       @param bilinearform: inner product C{<x,r>}
491       @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
492                           type like argument C{x} and C{r} is. The returned value
493                           is a C{float}.
494       @param atol: absolute tolerance
495       @type atol: non-negative C{float}
496       @param rtol: relative tolerance
497       @type rtol: non-negative C{float}
498       @param iter_max: maximum number of iteration steps
499     @type iter_max: C{int}     @type iter_max: C{int}
500     @return: the solution approximation and the corresponding residual     @return: the solution approximation and the corresponding residual
501     @rtype: C{tuple}     @rtype: C{tuple}
502     @warning: C{b} and C{x} are altered.     @warning: C{r} and C{x} are altered.
503     """     """
504     iter=0     iter=0
    if x==None:  
       x=0*b  
    else:  
       b += (-1)*Aprod(x)  
    r=b  
505     rhat=Msolve(r)     rhat=Msolve(r)
506     d = rhat     d = rhat
507     rhat_dot_r = bilinearform(rhat, r)     rhat_dot_r = bilinearform(rhat, r)
508     if rhat_dot_r<0: raise NegativeNorm,"negative norm."     if rhat_dot_r<0: raise NegativeNorm,"negative norm."
509       norm_r0=math.sqrt(rhat_dot_r)
510       atol2=atol+rtol*norm_r0
511       if atol2<=0:
512          raise ValueError,"Non-positive tolarance."
513       atol2=max(atol2, 100. * util.EPSILON * norm_r0)
514    
515       if verbose: print "PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_r0, atol2)
516    
517    
518     while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):     while not math.sqrt(rhat_dot_r) <= atol2:
519         iter+=1         iter+=1
520         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
521    
522         q=Aprod(d)         q=Aprod(d)
523         alpha = rhat_dot_r / bilinearform(d, q)         alpha = rhat_dot_r / bilinearform(d, q)
524         x += alpha * d         x += alpha * d
525         r += (-alpha) * q         if isinstance(q,ArithmeticTuple):
526           r += q * (-alpha)      # Doing it the other way calls the float64.__mul__ not AT.__rmul__
527           else:
528               r += (-alpha) * q
529         rhat=Msolve(r)         rhat=Msolve(r)
530         rhat_dot_r_new = bilinearform(rhat, r)         rhat_dot_r_new = bilinearform(rhat, r)
531         beta = rhat_dot_r_new / rhat_dot_r         beta = rhat_dot_r_new / rhat_dot_r
# Line 557  type like argument C{x}. Line 534  type like argument C{x}.
534    
535         rhat_dot_r = rhat_dot_r_new         rhat_dot_r = rhat_dot_r_new
536         if rhat_dot_r<0: raise NegativeNorm,"negative norm."         if rhat_dot_r<0: raise NegativeNorm,"negative norm."
537           if verbose: print "PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))
538     return x,r     if verbose: print "PCG: tolerance reached after %s steps."%iter
539       return x,r,math.sqrt(rhat_dot_r)
540    
541  class Defect(object):  class Defect(object):
542      """      """
543      defines a non-linear defect F(x) of a variable x      Defines a non-linear defect F(x) of a variable x.
544      """      """
545      def __init__(self):      def __init__(self):
546          """          """
547          initialize defect          Initializes defect.
548          """          """
549          self.setDerivativeIncrementLength()          self.setDerivativeIncrementLength()
550    
551      def bilinearform(self, x0, x1):      def bilinearform(self, x0, x1):
552          """          """
553          returns the inner product of x0 and x1          Returns the inner product of x0 and x1
554          @param x0: a value for x  
555          @param x1: a value for x          @param x0: value for x0
556            @param x1: value for x1
557          @return: the inner product of x0 and x1          @return: the inner product of x0 and x1
558          @rtype: C{float}          @rtype: C{float}
559          """          """
560          return 0          return 0
561          
562      def norm(self,x):      def norm(self,x):
563          """          """
564          the norm of argument C{x}          Returns the norm of argument C{x}.
565    
566          @param x: a value for x          @param x: a value
567          @return: norm of argument x          @return: norm of argument x
568          @rtype: C{float}          @rtype: C{float}
569          @note: by default C{sqrt(self.bilinearform(x,x)} is retrurned.          @note: by default C{sqrt(self.bilinearform(x,x)} is returned.
570          """          """
571          s=self.bilinearform(x,x)          s=self.bilinearform(x,x)
572          if s<0: raise NegativeNorm,"negative norm."          if s<0: raise NegativeNorm,"negative norm."
573          return math.sqrt(s)          return math.sqrt(s)
574    
   
575      def eval(self,x):      def eval(self,x):
576          """          """
577          returns the value F of a given x          Returns the value F of a given C{x}.
578    
579          @param x: value for which the defect C{F} is evalulated.          @param x: value for which the defect C{F} is evaluated
580          @return: value of the defect at C{x}          @return: value of the defect at C{x}
581          """          """
582          return 0          return 0
# Line 608  class Defect(object): Line 586  class Defect(object):
586    
587      def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):      def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):
588          """          """
589          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
590          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
591            direction of v with x as a starting point.
592    
593          @param inc: relative increment length          @param inc: relative increment length
594          @type inc: positive C{float}          @type inc: positive C{float}
# Line 619  class Defect(object): Line 598  class Defect(object):
598    
599      def getDerivativeIncrementLength(self):      def getDerivativeIncrementLength(self):
600          """          """
601          returns the relative increment length used to approximate the derivative of the defect          Returns the relative increment length used to approximate the
602            derivative of the defect.
603          @return: value of the defect at C{x}          @return: value of the defect at C{x}
604          @rtype: positive C{float}          @rtype: positive C{float}
605          """          """
# Line 627  class Defect(object): Line 607  class Defect(object):
607    
608      def derivative(self, F0, x0, v, v_is_normalised=True):      def derivative(self, F0, x0, v, v_is_normalised=True):
609          """          """
610          returns the directional derivative at x0 in the direction of v          Returns the directional derivative at C{x0} in the direction of C{v}.
611    
612          @param F0: value of this defect at x0          @param F0: value of this defect at x0
613          @param x0: value at which derivative is calculated.          @param x0: value at which derivative is calculated
614          @param v: direction          @param v: direction
615          @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
616                                    (self.norm(v)=0)
617          @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}
618          @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
619          maybe oepsnew verwritten to use exact evalution.                 used but this method maybe overwritten to use exact evaluation.
620          """          """
621          normx=self.norm(x0)          normx=self.norm(x0)
622          if normx>0:          if normx>0:
# Line 651  class Defect(object): Line 632  class Defect(object):
632          F1=self.eval(x0 + epsnew * v)          F1=self.eval(x0 + epsnew * v)
633          return (F1-F0)/epsnew          return (F1-F0)/epsnew
634    
635  ######################################      ######################################
636  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):
637     """     """
638     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
639       criterion:
640    
641     M{norm(F(x) <= atol + rtol * norm(F(x0)}     M{norm(F(x) <= atol + rtol * norm(F(x0)}
642      
643     where M{x0} is the initial guess.     where M{x0} is the initial guess.
644    
645     @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
646                      M{norm} used in the stopping criterion.
647     @type defect: L{Defect}     @type defect: L{Defect}
648     @param x: initial guess for the solution, C{x} is altered.     @param x: initial guess for the solution, C{x} is altered.
649     @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
650     @param iter_max: maximum number of iteration steps              C{numpy.ndarray}, L{Data}
651       @param iter_max: maximum number of iteration steps
652     @type iter_max: positive C{int}     @type iter_max: positive C{int}
653     @param sub_iter_max:     @param sub_iter_max: maximum number of inner iteration steps
654     @type sub_iter_max:     @type sub_iter_max: positive C{int}
655     @param atol: absolute tolerance for the solution     @param atol: absolute tolerance for the solution
656     @type atol: positive C{float}     @type atol: positive C{float}
657     @param rtol: relative tolerance for the solution     @param rtol: relative tolerance for the solution
658     @type rtol: positive C{float}     @type rtol: positive C{float}
659     @param gamma: tolerance safety factor for inner iteration     @param gamma: tolerance safety factor for inner iteration
660     @type gamma: positive C{float}, less than 1     @type gamma: positive C{float}, less than 1
661     @param sub_tol_max: upper bound for inner tolerance.     @param sub_tol_max: upper bound for inner tolerance
662     @type sub_tol_max: positive C{float}, less than 1     @type sub_tol_max: positive C{float}, less than 1
663     @return: an approximation of the solution with the desired accuracy     @return: an approximation of the solution with the desired accuracy
664     @rtype: same type as the initial guess.     @rtype: same type as the initial guess
665     """     """
666     lmaxit=iter_max     lmaxit=iter_max
667     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 681  def NewtonGMRES(defect, x, iter_max=100,
681     # main iteration loop     # main iteration loop
682     #     #
683     while not fnrm<=stop_tol:     while not fnrm<=stop_tol:
684              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
685              #              #
686          #   adjust sub_tol_          #   adjust sub_tol_
687          #          #
688              if iter > 1:              if iter > 1:
689             rat=fnrm/fnrmo             rat=fnrm/fnrmo
# Line 712  def NewtonGMRES(defect, x, iter_max=100, Line 696  def NewtonGMRES(defect, x, iter_max=100,
696              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown
697              #     if iter_restart -1 is returned as sub_iter              #     if iter_restart -1 is returned as sub_iter
698              #     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
699              #              #
700              #                #
701              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
702              try:              try:
703                 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 718  def NewtonGMRES(defect, x, iter_max=100,
718    
719  def __givapp(c,s,vin):  def __givapp(c,s,vin):
720      """      """
721      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
722        C{vin}
723    
724      @warning: C{vin} is altered.      @warning: C{vin} is altered.
725      """      """
726      vrot=vin      vrot=vin
727      if isinstance(c,float):      if isinstance(c,float):
728          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]]
729      else:      else:
730          for i in range(len(c)):          for i in range(len(c)):
731              w1=c[i]*vrot[i]-s[i]*vrot[i+1]              w1=c[i]*vrot[i]-s[i]*vrot[i+1]
732          w2=s[i]*vrot[i]+c[i]*vrot[i+1]          w2=s[i]*vrot[i]+c[i]*vrot[i+1]
733              vrot[i:i+2]=w1,w2              vrot[i]=w1
734                vrot[i+1]=w2
735      return vrot      return vrot
736    
737  def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):  def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
738     h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)     h=numpy.zeros((iter_restart,iter_restart),numpy.float64)
739     c=numarray.zeros(iter_restart,numarray.Float64)     c=numpy.zeros(iter_restart,numpy.float64)
740     s=numarray.zeros(iter_restart,numarray.Float64)     s=numpy.zeros(iter_restart,numpy.float64)
741     g=numarray.zeros(iter_restart,numarray.Float64)     g=numpy.zeros(iter_restart,numpy.float64)
742     v=[]     v=[]
743    
744     rho=defect.norm(F0)     rho=defect.norm(F0)
745     if rho<=0.: return x0*0     if rho<=0.: return x0*0
746      
747     v.append(-F0/rho)     v.append(-F0/rho)
748     g[0]=rho     g[0]=rho
749     iter=0     iter=0
750     while rho > atol and iter<iter_restart-1:     while rho > atol and iter<iter_restart-1:
751            if iter  >= iter_max:
752      if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max              raise MaxIterReached,"maximum number of %s steps reached."%iter_max
753    
754          p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)          p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
755      v.append(p)          v.append(p)
756    
757      v_norm1=defect.norm(v[iter+1])          v_norm1=defect.norm(v[iter+1])
758    
759          # Modified Gram-Schmidt          # Modified Gram-Schmidt
760      for j in range(iter+1):          for j in range(iter+1):
761           h[j][iter]=defect.bilinearform(v[j],v[iter+1])                h[j,iter]=defect.bilinearform(v[j],v[iter+1])
762           v[iter+1]-=h[j][iter]*v[j]              v[iter+1]-=h[j,iter]*v[j]
763          
764      h[iter+1][iter]=defect.norm(v[iter+1])          h[iter+1,iter]=defect.norm(v[iter+1])
765      v_norm2=h[iter+1][iter]          v_norm2=h[iter+1,iter]
766    
767          # Reorthogonalize if needed          # Reorthogonalize if needed
768      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)
769          for j in range(iter+1):                for j in range(iter+1):
770             hr=defect.bilinearform(v[j],v[iter+1])                  hr=defect.bilinearform(v[j],v[iter+1])
771                 h[j][iter]=h[j][iter]+hr                  h[j,iter]=h[j,iter]+hr
772                 v[iter+1] -= hr*v[j]                  v[iter+1] -= hr*v[j]
773    
774          v_norm2=defect.norm(v[iter+1])              v_norm2=defect.norm(v[iter+1])
775          h[iter+1][iter]=v_norm2              h[iter+1,iter]=v_norm2
776          #   watch out for happy breakdown          #   watch out for happy breakdown
777          if not v_norm2 == 0:          if not v_norm2 == 0:
778                  v[iter+1]=v[iter+1]/h[iter+1][iter]              v[iter+1]=v[iter+1]/h[iter+1,iter]
779    
780          #   Form and store the information for the new Givens rotation          #   Form and store the information for the new Givens rotation
781      if iter > 0 :          if iter > 0 :
782          hhat=numarray.zeros(iter+1,numarray.Float64)              hhat=numpy.zeros(iter+1,numpy.float64)
783          for i in range(iter+1) : hhat[i]=h[i][iter]              for i in range(iter+1) : hhat[i]=h[i,iter]
784          hhat=__givapp(c[0:iter],s[0:iter],hhat);              hhat=__givapp(c[0:iter],s[0:iter],hhat);
785              for i in range(iter+1) : h[i][iter]=hhat[i]              for i in range(iter+1) : h[i,iter]=hhat[i]
786    
787      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])
788    
789      if mu!=0 :          if mu!=0 :
790          c[iter]=h[iter][iter]/mu              c[iter]=h[iter,iter]/mu
791          s[iter]=-h[iter+1][iter]/mu              s[iter]=-h[iter+1,iter]/mu
792          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]
793          h[iter+1][iter]=0.0              h[iter+1,iter]=0.0
794          g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])              gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
795                g[iter]=gg[0]
796                g[iter+1]=gg[1]
797    
798          # Update the residual norm          # Update the residual norm
799          rho=abs(g[iter+1])          rho=abs(g[iter+1])
800      iter+=1          iter+=1
801    
802     # At this point either iter > iter_max or rho < tol.     # At this point either iter > iter_max or rho < tol.
803     # It's time to compute x and leave.             # It's time to compute x and leave.
804     if iter > 0 :     if iter > 0 :
805       y=numarray.zeros(iter,numarray.Float64)           y=numpy.zeros(iter,numpy.float64)
806       y[iter-1] = g[iter-1] / h[iter-1][iter-1]       y[iter-1] = g[iter-1] / h[iter-1,iter-1]
807       if iter > 1 :         if iter > 1 :
808          i=iter-2            i=iter-2
809          while i>=0 :          while i>=0 :
810            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]
811            i=i-1            i=i-1
812       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
813       for i in range(iter-1):       for i in range(iter-1):
814      xhat += v[i]*y[i]      xhat += v[i]*y[i]
815     else :     else :
816        xhat=v[0] * 0        xhat=v[0] * 0
817    
818     if iter<iter_restart-1:     if iter<iter_restart-1:
819        stopped=iter        stopped=iter
820     else:     else:
821        stopped=-1        stopped=-1
822    
823     return xhat,stopped     return xhat,stopped
824    
825  ##############################################  def GMRES(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100, iter_restart=20, verbose=False,P_R=None):
826  def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):     """
827  ################################################     Solver for
828    
829       M{Ax=b}
830    
831       with a general operator A (more details required!).
832       It uses the generalized minimum residual method (GMRES).
833    
834       The iteration is terminated if
835    
836       M{|r| <= atol+rtol*|r0|}
837    
838       where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
839    
840       M{|r| = sqrt( bilinearform(r,r))}
841    
842       @param r: initial residual M{r=b-Ax}. C{r} is altered.
843       @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
844       @param x: an initial guess for the solution
845       @type x: same like C{r}
846       @param Aprod: returns the value Ax
847       @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
848                    argument C{x}. The returned object needs to be of the same
849                    type like argument C{r}.
850       @param bilinearform: inner product C{<x,r>}
851       @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
852                           type like argument C{x} and C{r}. The returned value is
853                           a C{float}.
854       @param atol: absolute tolerance
855       @type atol: non-negative C{float}
856       @param rtol: relative tolerance
857       @type rtol: non-negative C{float}
858       @param iter_max: maximum number of iteration steps
859       @type iter_max: C{int}
860       @param iter_restart: in order to save memory the orthogonalization process
861                            is terminated after C{iter_restart} steps and the
862                            iteration is restarted.
863       @type iter_restart: C{int}
864       @return: the solution approximation and the corresponding residual
865       @rtype: C{tuple}
866       @warning: C{r} and C{x} are altered.
867       """
868     m=iter_restart     m=iter_restart
869       restarted=False
870     iter=0     iter=0
871     xc=x     if rtol>0:
872          r_dot_r = bilinearform(r, r)
873          if r_dot_r<0: raise NegativeNorm,"negative norm."
874          atol2=atol+rtol*math.sqrt(r_dot_r)
875          if verbose: print "GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)
876       else:
877          atol2=atol
878          if verbose: print "GMRES: absolute tolerance = %e"%atol2
879       if atol2<=0:
880          raise ValueError,"Non-positive tolarance."
881    
882     while True:     while True:
883        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
884        xc,stopped=__GMRESm(b*1, Aprod, Msolve, bilinearform, stoppingcriterium, x=xc*1, iter_max=iter_max-iter, iter_restart=m)        if restarted:
885             r2 = r-Aprod(x-x2)
886          else:
887             r2=1*r
888          x2=x*1.
889          x,stopped=_GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose,P_R=P_R)
890          iter+=iter_restart
891        if stopped: break        if stopped: break
892        iter+=iter_restart            if verbose: print "GMRES: restart."
893     return xc        restarted=True
894       if verbose: print "GMRES: tolerance has been reached."
895       return x
896    
897  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):
898     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)  
899    
900     if x==None:     h=numpy.zeros((iter_restart+1,iter_restart),numpy.float64)
901        x=0*b     c=numpy.zeros(iter_restart,numpy.float64)
902     else:     s=numpy.zeros(iter_restart,numpy.float64)
903        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)  
904     v=[]     v=[]
905    
906       r_dot_r = bilinearform(r, r)
907       if r_dot_r<0: raise NegativeNorm,"negative norm."
908     rho=math.sqrt(r_dot_r)     rho=math.sqrt(r_dot_r)
909      
910     v.append(r/rho)     v.append(r/rho)
911     g[0]=rho     g[0]=rho
912    
913     while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):     if verbose: print "GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)
914       while not (rho<=atol or iter==iter_restart):
915    
916      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
917    
918      p=Msolve(Aprod(v[iter]))          if P_R!=None:
919                p=Aprod(P_R(v[iter]))
920            else:
921            p=Aprod(v[iter])
922      v.append(p)      v.append(p)
923    
924      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
925    
926  # Modified Gram-Schmidt  # Modified Gram-Schmidt
927      for j in range(iter+1):      for j in range(iter+1):
928        h[j][iter]=bilinearform(v[j],v[iter+1])          h[j,iter]=bilinearform(v[j],v[iter+1])
929        v[iter+1]-=h[j][iter]*v[j]        v[iter+1]-=h[j,iter]*v[j]
930          
931      h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))      h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
932      v_norm2=h[iter+1][iter]      v_norm2=h[iter+1,iter]
933    
934  # Reorthogonalize if needed  # Reorthogonalize if needed
935      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)
936       for j in range(iter+1):         for j in range(iter+1):
937          hr=bilinearform(v[j],v[iter+1])          hr=bilinearform(v[j],v[iter+1])
938              h[j][iter]=h[j][iter]+hr              h[j,iter]=h[j,iter]+hr
939              v[iter+1] -= hr*v[j]              v[iter+1] -= hr*v[j]
940    
941       v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
942       h[iter+1][iter]=v_norm2       h[iter+1,iter]=v_norm2
943    
944  #   watch out for happy breakdown  #   watch out for happy breakdown
945          if not v_norm2 == 0:          if not v_norm2 == 0:
946           v[iter+1]=v[iter+1]/h[iter+1][iter]           v[iter+1]=v[iter+1]/h[iter+1,iter]
947    
948  #   Form and store the information for the new Givens rotation  #   Form and store the information for the new Givens rotation
949      if iter > 0 :      if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])
950          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])  
951    
952      if mu!=0 :      if mu!=0 :
953          c[iter]=h[iter][iter]/mu          c[iter]=h[iter,iter]/mu
954          s[iter]=-h[iter+1][iter]/mu          s[iter]=-h[iter+1,iter]/mu
955          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]
956          h[iter+1][iter]=0.0          h[iter+1,iter]=0.0
957          g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])                  gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
958                    g[iter]=gg[0]
959                    g[iter+1]=gg[1]
960  # Update the residual norm  # Update the residual norm
961                  
962          rho=abs(g[iter+1])          rho=abs(g[iter+1])
963            if verbose: print "GMRES: iteration step %s: residual %e"%(iter,rho)
964      iter+=1      iter+=1
965    
966  # At this point either iter > iter_max or rho < tol.  # At this point either iter > iter_max or rho < tol.
967  # It's time to compute x and leave.          # It's time to compute x and leave.
968    
969     if iter > 0 :     if verbose: print "GMRES: iteration stopped after %s step."%iter
970       y=numarray.zeros(iter,numarray.Float64)         if iter > 0 :
971       y[iter-1] = g[iter-1] / h[iter-1][iter-1]       y=numpy.zeros(iter,numpy.float64)
972       if iter > 1 :         y[iter-1] = g[iter-1] / h[iter-1,iter-1]
973          i=iter-2         if iter > 1 :
974            i=iter-2
975          while i>=0 :          while i>=0 :
976            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]
977            i=i-1            i=i-1
978       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
979       for i in range(iter-1):       for i in range(iter-1):
980      xhat += v[i]*y[i]      xhat += v[i]*y[i]
981     else : xhat=v[0]     else:
982         xhat=v[0] * 0
983     x += xhat     if P_R!=None:
984     if iter<iter_restart-1:        x += P_R(xhat)
985        stopped=True     else:
986     else:        x += xhat
987       if iter<iter_restart-1:
988          stopped=True
989       else:
990        stopped=False        stopped=False
991    
992     return x,stopped     return x,stopped
993    
994  #################################################  def MINRES(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
995  def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):      """
996  #################################################      Solver for
997      #  
998      #  minres solves the system of linear equations Ax = b      M{Ax=b}
999      #  where A is a symmetric matrix (possibly indefinite or singular)  
1000      #  and b is a given vector.      with a symmetric and positive definite operator A (more details required!).
1001      #        It uses the minimum residual method (MINRES) with preconditioner M
1002      #  "A" may be a dense or sparse matrix (preferably sparse!)      providing an approximation of A.
1003      #  or the name of a function such that  
1004      #               y = A(x)      The iteration is terminated if
1005      #  returns the product y = Ax for any given vector x.  
1006      #      M{|r| <= atol+rtol*|r0|}
1007      #  "M" defines a positive-definite preconditioner M = C C'.  
1008      #  "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
1009      #  or the name of a function such that  
1010      #  solves the system My = x for any given vector x.      M{|r| = sqrt( bilinearform(Msolve(r),r))}
1011      #  
1012      #      For details on the preconditioned conjugate gradient method see the book:
1013        
1014        I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
1015        T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
1016        C. Romine, and H. van der Vorst}.
1017    
1018        @param r: initial residual M{r=b-Ax}. C{r} is altered.
1019        @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1020        @param x: an initial guess for the solution
1021        @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1022        @param Aprod: returns the value Ax
1023        @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
1024                     argument C{x}. The returned object needs to be of the same
1025                     type like argument C{r}.
1026        @param Msolve: solves Mx=r
1027        @type Msolve: function C{Msolve(r)} where C{r} is of the same type like
1028                      argument C{r}. The returned object needs to be of the same
1029                      type like argument C{x}.
1030        @param bilinearform: inner product C{<x,r>}
1031        @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
1032                            type like argument C{x} and C{r} is. The returned value
1033                            is a C{float}.
1034        @param atol: absolute tolerance
1035        @type atol: non-negative C{float}
1036        @param rtol: relative tolerance
1037        @type rtol: non-negative C{float}
1038        @param iter_max: maximum number of iteration steps
1039        @type iter_max: C{int}
1040        @return: the solution approximation and the corresponding residual
1041        @rtype: C{tuple}
1042        @warning: C{r} and C{x} are altered.
1043        """
1044      #------------------------------------------------------------------      #------------------------------------------------------------------
1045      # Set up y and v for the first Lanczos vector v1.      # Set up y and v for the first Lanczos vector v1.
1046      # y  =  beta1 P' v1,  where  P = C**(-1).      # y  =  beta1 P' v1,  where  P = C**(-1).
1047      # v is really P' v1.      # v is really P' v1.
1048      #------------------------------------------------------------------      #------------------------------------------------------------------
1049      if x==None:      r1    = r
1050        x=0*b      y = Msolve(r)
1051      else:      beta1 = bilinearform(y,r)
       b += (-1)*Aprod(x)  
1052    
     r1    = b  
     y = Msolve(b)  
     beta1 = bilinearform(y,b)  
   
1053      if beta1< 0: raise NegativeNorm,"negative norm."      if beta1< 0: raise NegativeNorm,"negative norm."
1054    
1055      #  If b = 0 exactly, stop with x = 0.      #  If r = 0 exactly, stop with x
1056      if beta1==0: return x*0.      if beta1==0: return x
1057    
1058      if beta1> 0:      if beta1> 0: beta1  = math.sqrt(beta1)
       beta1  = math.sqrt(beta1)        
1059    
1060      #------------------------------------------------------------------      #------------------------------------------------------------------
1061      # Initialize quantities.      # Initialize quantities.
# Line 1008  def MINRES(b, Aprod, Msolve, bilinearfor Line 1075  def MINRES(b, Aprod, Msolve, bilinearfor
1075      ynorm2 = 0      ynorm2 = 0
1076      cs     = -1      cs     = -1
1077      sn     = 0      sn     = 0
1078      w      = b*0.      w      = r*0.
1079      w2     = b*0.      w2     = r*0.
1080      r2     = r1      r2     = r1
1081      eps    = 0.0001      eps    = 0.0001
1082    
1083      #---------------------------------------------------------------------      #---------------------------------------------------------------------
1084      # Main iteration loop.      # Main iteration loop.
1085      # --------------------------------------------------------------------      # --------------------------------------------------------------------
1086      while not stoppingcriterium(rnorm,Anorm*ynorm,'MINRES'):    #  checks ||r|| < (||A|| ||x||) * TOL      while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL
1087    
1088      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
1089          iter    = iter  +  1          iter    = iter  +  1
# Line 1035  def MINRES(b, Aprod, Msolve, bilinearfor Line 1102  def MINRES(b, Aprod, Msolve, bilinearfor
1102          #-----------------------------------------------------------------          #-----------------------------------------------------------------
1103          s = 1/beta                 # Normalize previous vector (in y).          s = 1/beta                 # Normalize previous vector (in y).
1104          v = s*y                    # v = vk if P = I          v = s*y                    # v = vk if P = I
1105        
1106          y      = Aprod(v)          y      = Aprod(v)
1107        
1108          if iter >= 2:          if iter >= 2:
1109            y = y - (beta/oldb)*r1            y = y - (beta/oldb)*r1
1110    
1111          alfa   = bilinearform(v,y)              # alphak          alfa   = bilinearform(v,y)              # alphak
1112          y      += (- alfa/beta)*r2          y      += (- alfa/beta)*r2
1113          r1     = r2          r1     = r2
1114          r2     = y          r2     = y
1115          y = Msolve(r2)          y = Msolve(r2)
# Line 1052  def MINRES(b, Aprod, Msolve, bilinearfor Line 1119  def MINRES(b, Aprod, Msolve, bilinearfor
1119    
1120          beta   = math.sqrt( beta )          beta   = math.sqrt( beta )
1121          tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta          tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
1122            
1123          if iter==1:                 # Initialize a few things.          if iter==1:                 # Initialize a few things.
1124            gmax   = abs( alfa )      # alpha1            gmax   = abs( alfa )      # alpha1
1125            gmin   = gmax             # alpha1            gmin   = gmax             # alpha1
# Line 1060  def MINRES(b, Aprod, Msolve, bilinearfor Line 1127  def MINRES(b, Aprod, Msolve, bilinearfor
1127          # Apply previous rotation Qk-1 to get          # Apply previous rotation Qk-1 to get
1128          #   [deltak epslnk+1] = [cs  sn][dbark    0   ]          #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
1129          #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].          #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
1130        
1131          oldeps = epsln          oldeps = epsln
1132          delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak          delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
1133          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 1137  def MINRES(b, Aprod, Msolve, bilinearfor
1137          # Compute the next plane rotation Qk          # Compute the next plane rotation Qk
1138    
1139          gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak          gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
1140          gamma  = max(gamma,eps)          gamma  = max(gamma,eps)
1141          cs     = gbar / gamma             # ck          cs     = gbar / gamma             # ck
1142          sn     = beta / gamma             # sk          sn     = beta / gamma             # sk
1143          phi    = cs * phibar              # phik          phi    = cs * phibar              # phik
# Line 1078  def MINRES(b, Aprod, Msolve, bilinearfor Line 1145  def MINRES(b, Aprod, Msolve, bilinearfor
1145    
1146          # Update  x.          # Update  x.
1147    
1148          denom = 1/gamma          denom = 1/gamma
1149          w1    = w2          w1    = w2
1150          w2    = w          w2    = w
1151          w     = (v - oldeps*w1 - delta*w2) * denom          w     = (v - oldeps*w1 - delta*w2) * denom
1152          x     +=  phi*w          x     +=  phi*w
1153    
# Line 1095  def MINRES(b, Aprod, Msolve, bilinearfor Line 1162  def MINRES(b, Aprod, Msolve, bilinearfor
1162    
1163          # Estimate various norms and test for convergence.          # Estimate various norms and test for convergence.
1164    
1165          Anorm  = math.sqrt( tnorm2 )          Anorm  = math.sqrt( tnorm2 )
1166          ynorm  = math.sqrt( ynorm2 )          ynorm  = math.sqrt( ynorm2 )
1167    
1168          rnorm  = phibar          rnorm  = phibar
1169    
1170      return x      return x
1171    
1172  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):
1173      """
1174      Solver for
1175    
1176  # TFQMR solver for linear systems    M{Ax=b}
 #  
 #  
 # initialization  
 #  
   errtol = math.sqrt(bilinearform(b,b))  
   norm_b=errtol  
   kmax  = iter_max  
   error = []  
   
   if math.sqrt(bilinearform(x,x)) != 0.0:  
     r = b - Aprod(x)  
   else:  
     r = b  
1177    
1178    r=Msolve(r)    with a general operator A (more details required!).
1179      It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).
1180    
1181      The iteration is terminated if
1182    
1183      M{|r| <= atol+rtol*|r0|}
1184    
1185      where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
1186    
1187      M{|r| = sqrt( bilinearform(r,r))}
1188    
1189      @param r: initial residual M{r=b-Ax}. C{r} is altered.
1190      @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1191      @param x: an initial guess for the solution
1192      @type x: same like C{r}
1193      @param Aprod: returns the value Ax
1194      @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
1195                   argument C{x}. The returned object needs to be of the same type
1196                   like argument C{r}.
1197      @param bilinearform: inner product C{<x,r>}
1198      @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
1199                          type like argument C{x} and C{r}. The returned value is
1200                          a C{float}.
1201      @param atol: absolute tolerance
1202      @type atol: non-negative C{float}
1203      @param rtol: relative tolerance
1204      @type rtol: non-negative C{float}
1205      @param iter_max: maximum number of iteration steps
1206      @type iter_max: C{int}
1207      @rtype: C{tuple}
1208      @warning: C{r} and C{x} are altered.
1209      """
1210    u1=0    u1=0
1211    u2=0    u2=0
1212    y1=0    y1=0
1213    y2=0    y2=0
1214    
1215    w = r    w = r
1216    y1 = r    y1 = r
1217    iter = 0    iter = 0
1218    d = 0    d = 0
1219        v = Aprod(y1)
   v = Msolve(Aprod(y1))  
1220    u1 = v    u1 = v
1221      
1222    theta = 0.0;    theta = 0.0;
1223    eta = 0.0;    eta = 0.0;
1224    tau = math.sqrt(bilinearform(r,r))    rho=bilinearform(r,r)
1225    error = [ error, tau ]    if rho < 0: raise NegativeNorm,"negative norm."
1226    rho = tau * tau    tau = math.sqrt(rho)
1227    m=1    norm_r0=tau
1228  #    while tau>atol+rtol*norm_r0:
 #  TFQMR iteration  
 #  
 #  while ( iter < kmax-1 ):  
     
   while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):  
1229      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
1230    
1231      sigma = bilinearform(r,v)      sigma = bilinearform(r,v)
1232        if sigma == 0.0: raise IterationBreakDown,'TFQMR breakdown, sigma=0'
     if ( sigma == 0.0 ):  
       raise 'TFQMR breakdown, sigma=0'  
       
1233    
1234      alpha = rho / sigma      alpha = rho / sigma
1235    
# Line 1162  def TFQMR(b, Aprod, Msolve, bilinearform Line 1239  def TFQMR(b, Aprod, Msolve, bilinearform
1239  #  #
1240        if ( j == 1 ):        if ( j == 1 ):
1241          y2 = y1 - alpha * v          y2 = y1 - alpha * v
1242          u2 = Msolve(Aprod(y2))          u2 = Aprod(y2)
1243    
1244        m = 2 * (iter+1) - 2 + (j+1)        m = 2 * (iter+1) - 2 + (j+1)
1245        if j==0:        if j==0:
1246           w = w - alpha * u1           w = w - alpha * u1
1247           d = y1 + ( theta * theta * eta / alpha ) * d           d = y1 + ( theta * theta * eta / alpha ) * d
1248        if j==1:        if j==1:
# Line 1180  def TFQMR(b, Aprod, Msolve, bilinearform Line 1257  def TFQMR(b, Aprod, Msolve, bilinearform
1257  #  #
1258  #  Try to terminate the iteration at each pass through the loop  #  Try to terminate the iteration at each pass through the loop
1259  #  #
1260       # 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'  
       
1261    
1262      rhon = bilinearform(r,w)      rhon = bilinearform(r,w)
1263      beta = rhon / rho;      beta = rhon / rho;
1264      rho = rhon;      rho = rhon;
1265      y1 = w + beta * y2;      y1 = w + beta * y2;
1266      u1 = Msolve(Aprod(y1))      u1 = Aprod(y1)
1267      v = u1 + beta * ( u2 + beta * v )      v = u1 + beta * ( u2 + beta * v )
1268      error = [ error, tau ]  
1269      total_iters = iter      iter += 1
       
     iter = iter + 1  
1270    
1271    return x    return x
1272    
# Line 1208  def TFQMR(b, Aprod, Msolve, bilinearform Line 1275  def TFQMR(b, Aprod, Msolve, bilinearform
1275    
1276  class ArithmeticTuple(object):  class ArithmeticTuple(object):
1277     """     """
1278     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
1279       ArithmeticTuple and C{a} is a float.
1280    
1281     example of usage:     Example of usage::
1282    
1283     from esys.escript import Data         from esys.escript import Data
1284     from numarray import array         from numpy import array
1285     a=Data(...)         a=Data(...)
1286     b=array([1.,4.])         b=array([1.,4.])
1287     x=ArithmeticTuple(a,b)         x=ArithmeticTuple(a,b)
1288     y=5.*x         y=5.*x
1289    
1290     """     """
1291     def __init__(self,*args):     def __init__(self,*args):
1292         """         """
1293         initialize object with elements args.         Initializes object with elements C{args}.
1294    
1295         @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
1296                        scaling (x=a*y)
1297         """         """
1298         self.__items=list(args)         self.__items=list(args)
1299    
1300     def __len__(self):     def __len__(self):
1301         """         """
1302         number of items         Returns the number of items.
1303    
1304         @return: number of items         @return: number of items
1305         @rtype: C{int}         @rtype: C{int}
# Line 1239  class ArithmeticTuple(object): Line 1308  class ArithmeticTuple(object):
1308    
1309     def __getitem__(self,index):     def __getitem__(self,index):
1310         """         """
1311         get an item         Returns item at specified position.
1312    
1313         @param index: item to be returned         @param index: index of item to be returned
1314         @type index: C{int}         @type index: C{int}
1315         @return: item with index C{index}         @return: item with index C{index}
1316         """         """
# Line 1249  class ArithmeticTuple(object): Line 1318  class ArithmeticTuple(object):
1318    
1319     def __mul__(self,other):     def __mul__(self,other):
1320         """         """
1321         scaling from the right         Scales by C{other} from the right.
1322    
1323         @param other: scaling factor         @param other: scaling factor
1324         @type other: C{float}         @type other: C{float}
# Line 1257  class ArithmeticTuple(object): Line 1326  class ArithmeticTuple(object):
1326         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1327         """         """
1328         out=[]         out=[]
1329         try:           try:
1330             l=len(other)             l=len(other)
1331             if l!=len(self):             if l!=len(self):
1332                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1333             for i in range(l): out.append(self[i]*other[i])             for i in range(l): out.append(self[i]*other[i])
1334         except TypeError:         except TypeError:
1335         for i in range(len(self)): out.append(self[i]*other)             for i in range(len(self)): out.append(self[i]*other)
1336         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1337    
1338     def __rmul__(self,other):     def __rmul__(self,other):
1339         """         """
1340         scaling from the left         Scales by C{other} from the left.
1341    
1342         @param other: scaling factor         @param other: scaling factor
1343         @type other: C{float}         @type other: C{float}
# Line 1276  class ArithmeticTuple(object): Line 1345  class ArithmeticTuple(object):
1345         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1346         """         """
1347         out=[]         out=[]
1348         try:           try:
1349             l=len(other)             l=len(other)
1350             if l!=len(self):             if l!=len(self):
1351                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1352             for i in range(l): out.append(other[i]*self[i])             for i in range(l): out.append(other[i]*self[i])
1353         except TypeError:         except TypeError:
1354         for i in range(len(self)): out.append(other*self[i])             for i in range(len(self)): out.append(other*self[i])
1355         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1356    
1357     def __div__(self,other):     def __div__(self,other):
1358         """         """
1359         dividing from the right         Scales by (1/C{other}) from the right.
1360    
1361         @param other: scaling factor         @param other: scaling factor
1362         @type other: C{float}         @type other: C{float}
# Line 1298  class ArithmeticTuple(object): Line 1367  class ArithmeticTuple(object):
1367    
1368     def __rdiv__(self,other):     def __rdiv__(self,other):
1369         """         """
1370         dividing from the left         Scales by (1/C{other}) from the left.
1371    
1372         @param other: scaling factor         @param other: scaling factor
1373         @type other: C{float}         @type other: C{float}
# Line 1306  class ArithmeticTuple(object): Line 1375  class ArithmeticTuple(object):
1375         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1376         """         """
1377         out=[]         out=[]
1378         try:           try:
1379             l=len(other)             l=len(other)
1380             if l!=len(self):             if l!=len(self):
1381                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1382             for i in range(l): out.append(other[i]/self[i])             for i in range(l): out.append(other[i]/self[i])
1383         except TypeError:         except TypeError:
1384         for i in range(len(self)): out.append(other/self[i])             for i in range(len(self)): out.append(other/self[i])
1385         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1386      
1387     def __iadd__(self,other):     def __iadd__(self,other):
1388         """         """
1389         in-place add of other to self         Inplace addition of C{other} to self.
1390    
1391         @param other: increment         @param other: increment
1392         @type other: C{ArithmeticTuple}         @type other: C{ArithmeticTuple}
1393         """         """
1394         if len(self) != len(other):         if len(self) != len(other):
1395             raise ValueError,"tuple length must match."             raise ValueError,"tuple lengths must match."
1396         for i in range(len(self)):         for i in range(len(self)):
1397             self.__items[i]+=other[i]             self.__items[i]+=other[i]
1398         return self         return self
1399    
1400     def __add__(self,other):     def __add__(self,other):
1401         """         """
1402         add other to self         Adds C{other} to self.
1403    
1404         @param other: increment         @param other: increment
1405         @type other: C{ArithmeticTuple}         @type other: C{ArithmeticTuple}
1406         """         """
1407         out=[]         out=[]
1408         try:           try:
1409             l=len(other)             l=len(other)
1410             if l!=len(self):             if l!=len(self):
1411                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1412             for i in range(l): out.append(self[i]+other[i])             for i in range(l): out.append(self[i]+other[i])
1413         except TypeError:         except TypeError:
1414         for i in range(len(self)): out.append(self[i]+other)             for i in range(len(self)): out.append(self[i]+other)
1415         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1416    
1417     def __sub__(self,other):     def __sub__(self,other):
1418         """         """
1419         subtract other from self         Subtracts C{other} from self.
1420    
1421         @param other: increment         @param other: decrement
1422         @type other: C{ArithmeticTuple}         @type other: C{ArithmeticTuple}
1423         """         """
1424         out=[]         out=[]
1425         try:           try:
1426             l=len(other)             l=len(other)
1427             if l!=len(self):             if l!=len(self):
1428                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1429             for i in range(l): out.append(self[i]-other[i])             for i in range(l): out.append(self[i]-other[i])
1430         except TypeError:         except TypeError:
1431         for i in range(len(self)): out.append(self[i]-other)             for i in range(len(self)): out.append(self[i]-other)
1432         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1433      
1434     def __isub__(self,other):     def __isub__(self,other):
1435         """         """
1436         subtract other from self         Inplace subtraction of C{other} from self.
1437    
1438         @param other: increment         @param other: decrement
1439         @type other: C{ArithmeticTuple}         @type other: C{ArithmeticTuple}
1440         """         """
1441         if len(self) != len(other):         if len(self) != len(other):
# Line 1377  class ArithmeticTuple(object): Line 1446  class ArithmeticTuple(object):
1446    
1447     def __neg__(self):     def __neg__(self):
1448         """         """
1449         negate         Negates values.
   
1450         """         """
1451         out=[]         out=[]
1452         for i in range(len(self)):         for i in range(len(self)):
# Line 1388  class ArithmeticTuple(object): Line 1456  class ArithmeticTuple(object):
1456    
1457  class HomogeneousSaddlePointProblem(object):  class HomogeneousSaddlePointProblem(object):
1458        """        """
1459        This provides a framwork for solving linear homogeneous saddle point problem of the form        This class provides a framework for solving linear homogeneous saddle
1460          point problems of the form::
1461    
1462               Av+B^*p=f            M{Av+B^*p=f}
1463               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.  
1464    
1465          for the unknowns M{v} and M{p} and given operators M{A} and M{B} and
1466          given right hand side M{f}. M{B^*} is the adjoint operator of M{B}.
1467        """        """
1468        def __init__(self,**kwargs):        def __init__(self,**kwargs):
1469          self.setTolerance()          self.setTolerance()
1470          self.setToleranceReductionFactor()          self.setAbsoluteTolerance()
1471            self.setSubProblemTolerance()
1472    
1473          #=============================================================
1474        def initialize(self):        def initialize(self):
1475          """          """
1476          initialize the problem (overwrite)          Initializes the problem (overwrite).
1477          """          """
1478          pass          pass
1479        def B(self,v):  
1480          def inner_pBv(self,p,Bv):
1481           """           """
1482           returns Bv (overwrite)           Returns inner product of element p and Bv (overwrite).
          @rtype: equal to the type of p  
1483    
1484           @note: boundary conditions on p should be zero!           @param p: a pressure increment
1485             @param v: a residual
1486             @return: inner product of element p and Bv
1487             @rtype: C{float}
1488             @note: used if PCG is applied.
1489           """           """
1490           pass           raise NotImplementedError,"no inner product for p and Bv implemented."
1491    
1492        def inner(self,p0,p1):        def inner_p(self,p0,p1):
1493           """           """
1494           returns inner product of two element p0 and p1  (overwrite)           Returns inner product of p0 and p1 (overwrite).
1495            
1496           @type p0: equal to the type of p           @param p0: a pressure
1497           @type p1: equal to the type of p           @param p1: a pressure
1498             @return: inner product of p0 and p1
1499           @rtype: C{float}           @rtype: C{float}
1500             """
1501             raise NotImplementedError,"no inner product for p implemented."
1502      
1503          def norm_v(self,v):
1504             """
1505             Returns the norm of v (overwrite).
1506    
1507           @rtype: equal to the type of p           @param v: a velovity
1508             @return: norm of v
1509             @rtype: non-negative C{float}
1510             """
1511             raise NotImplementedError,"no norm of v implemented."
1512          def getV(self, p, v0):
1513           """           """
1514           pass           return the value for v for a given p (overwrite)
1515    
1516             @param p: a pressure
1517             @param v0: a initial guess for the value v to return.
1518             @return: v given as M{v= A^{-1} (f-B^*p)}
1519             """
1520             raise NotImplementedError,"no v calculation implemented."
1521    
1522            
1523          def Bv(self,v):
1524            """
1525            Returns Bv (overwrite).
1526    
1527        def solve_A(self,u,p):          @rtype: equal to the type of p
1528            @note: boundary conditions on p should be zero!
1529            """
1530            raise NotImplementedError, "no operator B implemented."
1531    
1532          def norm_Bv(self,Bv):
1533            """
1534            Returns the norm of Bv (overwrite).
1535    
1536            @rtype: equal to the type of p
1537            @note: boundary conditions on p should be zero!
1538            """
1539            raise NotImplementedError, "no norm of Bv implemented."
1540    
1541          def solve_AinvBt(self,p):
1542           """           """
1543           solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
1544             (overwrite).
1545    
1546           @rtype: equal to the type of v           @param p: a pressure increment
1547             @return: the solution of M{Av=B^*p}
1548           @note: boundary conditions on v should be zero!           @note: boundary conditions on v should be zero!
1549           """           """
1550           pass           raise NotImplementedError,"no operator A implemented."
1551    
1552        def solve_prec(self,p):        def solve_prec(self,Bv):
1553           """           """
1554           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
1555             L{self.getSubProblemTolerance()} (overwrite).
1556    
1557           @rtype: equal to the type of p           @rtype: equal to the type of p
1558             @note: boundary conditions on p should be zero!
1559           """           """
1560           pass           raise NotImplementedError,"no preconditioner for Schur complement implemented."
1561          #=============================================================
1562          def __Aprod_PCG(self,p):
1563              dv=self.solve_AinvBt(p)
1564              return ArithmeticTuple(dv,self.Bv(dv))
1565    
1566        def stoppingcriterium(self,Bv,v,p):        def __inner_PCG(self,p,r):
1567           """           return self.inner_pBv(p,r[1])
          returns a True if iteration is terminated. (overwrite)  
1568    
1569           @rtype: C{bool}        def __Msolve_PCG(self,r):
1570           """            return self.solve_prec(r[1])
1571           pass        #=============================================================
1572                # rename solve_prec and change argument v to Bv
1573        def __inner(self,p,r):  # chnage the argument of inner_pBv to v->Bv
1574           return self.inner(p,r[1])  # add Bv
1575    # inner p still needed?
1576    # change norm_Bv argument to Bv
1577          def __Aprod_GMRES(self,p):
1578              return self.solve_prec(self.Bv(self.solve_AinvBt(p)))
1579          def __inner_GMRES(self,p0,p1):
1580             return self.inner_p(p0,p1)
1581          #=============================================================
1582          def norm_p(self,p):
1583              """
1584              calculates the norm of C{p}
1585              
1586              @param p: a pressure
1587              @return: the norm of C{p} using the inner product for pressure
1588              @rtype: C{float}
1589              """
1590              f=self.inner_p(p,p)
1591              if f<0: raise ValueError,"negative pressure norm."
1592              return math.sqrt(f)
1593              
1594    
1595        def __inner_p(self,p1,p2):        def solve(self,v,p,max_iter=20, verbose=False, show_details=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1596           return self.inner(p1,p2)           """
1597                   Solves the saddle point problem using initial guesses v and p.
       def __inner_a(self,a1,a2):  
          return self.inner_a(a1,a2)  
1598    
1599        def __inner_a1(self,a1,a2):           @param v: initial guess for velocity
1600           return self.inner(a1[1],a2[1])           @param p: initial guess for pressure
1601             @type v: L{Data}
1602             @type p: L{Data}
1603             @param usePCG: indicates the usage of the PCG rather than GMRES scheme.
1604             @param max_iter: maximum number of iteration steps per correction
1605                              attempt
1606             @param verbose: if True, shows information on the progress of the
1607                             saddlepoint problem solver.
1608             @param show_details: if True, shows details of the sub problem solver
1609             @param iter_restart: restart the iteration after C{iter_restart} steps
1610                                  (only used if useUzaw=False)
1611             @type usePCG: C{bool}
1612             @type max_iter: C{int}
1613             @type verbose: C{bool}
1614             @type show_details: C{bool}
1615             @type iter_restart: C{int}
1616             @rtype: C{tuple} of L{Data} objects
1617             """
1618             self.verbose=verbose
1619             self.show_details=show_details and self.verbose
1620             rtol=self.getTolerance()
1621             atol=self.getAbsoluteTolerance()
1622             correction_step=0
1623             converged=False
1624             while not converged:
1625                  # calculate velocity for current pressure:
1626                  v=self.getV(p,v)
1627                  Bv=self.Bv(v)
1628                  norm_v=self.norm_v(v)
1629                  norm_Bv=self.norm_Bv(Bv)
1630                  ATOL=norm_v*rtol+atol
1631                  if self.verbose: print "HomogeneousSaddlePointProblem: norm v= %e, norm_Bv= %e, tolerance = %e."%(norm_v, norm_Bv,ATOL)
1632                  if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."
1633                  if norm_Bv <= ATOL:
1634                     converged=True
1635                  else:
1636                     correction_step+=1
1637                     if correction_step>max_correction_steps:
1638                          raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1639                     dp=self.solve_prec(Bv)
1640                     if usePCG:
1641                       norm2=self.inner_pBv(dp,Bv)
1642                       if norm2<0: raise ValueError,"negative PCG norm."
1643                       norm2=math.sqrt(norm2)
1644                     else:
1645                       norm2=self.norm_p(dp)
1646                     ATOL_ITER=ATOL/norm_Bv*norm2*0.5
1647                     if self.verbose: print "HomogeneousSaddlePointProblem: tolerance for solver: %e"%ATOL_ITER
1648                     if usePCG:
1649                           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)
1650                     else:
1651                           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)
1652             if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached."
1653         return v,p
1654    
1655        def __stoppingcriterium(self,norm_r,r,p):        #========================================================================
1656            return self.stoppingcriterium(r[1],r[0],p)        def setTolerance(self,tolerance=1.e-4):
1657             """
1658             Sets the relative tolerance for (v,p).
1659    
1660        def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):           @param tolerance: tolerance to be used
1661            return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)           @type tolerance: non-negative C{float}
1662             """
1663             if tolerance<0:
1664                 raise ValueError,"tolerance must be positive."
1665             self.__rtol=tolerance
1666             self.setSubProblemTolerance()
1667    
       def setTolerance(self,tolerance=1.e-8):  
               self.__tol=tolerance  
1668        def getTolerance(self):        def getTolerance(self):
1669                return self.__tol           """
1670        def setToleranceReductionFactor(self,reduction=0.01):           Returns the relative tolerance.
               self.__reduction=reduction  
       def getSubProblemTolerance(self):  
               return self.__reduction*self.getTolerance()  
1671    
1672        def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):           @return: relative tolerance
1673                """           @rtype: C{float}
1674                solves the saddle point problem using initial guesses v and p.           """
1675             return self.__rtol
               @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)  
1676    
1677                # print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)        def setAbsoluteTolerance(self,tolerance=0.):
1678             """
1679             Sets the absolute tolerance.
1680    
1681            return u,p           @param tolerance: tolerance to be used
1682             @type tolerance: non-negative C{float}
1683             """
1684             if tolerance<0:
1685                 raise ValueError,"tolerance must be non-negative."
1686             self.__atol=tolerance
1687    
1688        def __Msolve(self,r):        def getAbsoluteTolerance(self):
1689            return self.solve_prec(r[1])           """
1690             Returns the absolute tolerance.
1691    
1692        def __Msolve2(self,r):           @return: absolute tolerance
1693            return self.solve_prec(r*1)           @rtype: C{float}
1694             """
1695             return self.__atol
1696    
1697        def __Mempty(self,r):        def setSubProblemTolerance(self,rtol=None):
1698            return r           """
1699             Sets the relative tolerance to solve the subproblem(s).
1700    
1701             @param rtol: relative tolerence
1702             @type rtol: positive C{float}
1703             """
1704             if rtol == None:
1705                  rtol=max(200.*util.EPSILON,self.getTolerance()**2)
1706             if rtol<=0:
1707                 raise ValueError,"tolerance must be positive."
1708             self.__sub_tol=rtol
1709    
1710        def __Aprod(self,p):        def getSubProblemTolerance(self):
1711            # return BA^-1B*p           """
1712            #solve Av =B^*p as Av =f-Az-B^*(-p)           Returns the subproblem reduction factor.
           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)  
1713    
1714             @return: subproblem reduction factor
1715             @rtype: C{float}
1716             """
1717             return self.__sub_tol
1718    
1719  def MaskFromBoundaryTag(domain,*tags):  def MaskFromBoundaryTag(domain,*tags):
1720     """     """
1721     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
1722     that touch the boundary tagged by tags.     one for samples that touch the boundary tagged by tags.
1723    
1724     usage: m=MaskFromBoundaryTag(domain,"left", "right")     Usage: m=MaskFromBoundaryTag(domain, "left", "right")
1725    
1726     @param domain: a given domain     @param domain: domain to be used
1727     @type domain: L{escript.Domain}     @type domain: L{escript.Domain}
1728     @param tags: boundray tags     @param tags: boundary tags
1729     @type tags: C{str}     @type tags: C{str}
1730     @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
1731                by any of the given tags
1732     @rtype: L{escript.Data} of rank 0     @rtype: L{escript.Data} of rank 0
1733     """     """
1734     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
# Line 1621  def MaskFromBoundaryTag(domain,*tags): Line 1736  def MaskFromBoundaryTag(domain,*tags):
1736     for t in tags: d.setTaggedValue(t,1.)     for t in tags: d.setTaggedValue(t,1.)
1737     pde.setValue(y=d)     pde.setValue(y=d)
1738     return util.whereNonZero(pde.getRightHandSide())     return util.whereNonZero(pde.getRightHandSide())
1739  #============================================================================================================================================  
1740    #==============================================================================
1741  class SaddlePointProblem(object):  class SaddlePointProblem(object):
1742     """     """
1743     This implements a solver for a saddlepoint problem     This implements a solver for a saddle point problem
1744    
1745     M{f(u,p)=0}     M{f(u,p)=0}
1746     M{g(u)=0}     M{g(u)=0}
# Line 1634  class SaddlePointProblem(object): Line 1750  class SaddlePointProblem(object):
1750     M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}     M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
1751     M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}     M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}
1752    
1753     where Q_f is an approximation of the Jacobiean A_f of f with respect to u  and Q_f is an approximation of     where Q_f is an approximation of the Jacobian A_f of f with respect to u and
1754     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_f is an approximation of A_g A_f^{-1} A_g with A_g is the Jacobian of g
1755     Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays     with respect to p. As a the construction of a 'proper' Q_g can be difficult,
1756       non-linear conjugate gradient method is applied to solve for p, so Q_g plays
1757     in fact the role of a preconditioner.     in fact the role of a preconditioner.
1758     """     """
1759     def __init__(self,verbose=False,*args):     def __init__(self,verbose=False,*args):
1760         """         """
1761         initializes the problem         Initializes the problem.
1762    
1763         @param verbose: switches on the printing out some information         @param verbose: if True, some information is printed in the process
1764         @type verbose: C{bool}         @type verbose: C{bool}
1765         @note: this method may be overwritten by a particular saddle point problem         @note: this method may be overwritten by a particular saddle point
1766                  problem
1767         """         """
1768         print "SaddlePointProblem should not be used anymore!"         print "SaddlePointProblem should not be used anymore!"
1769         if not isinstance(verbose,bool):         if not isinstance(verbose,bool):
# Line 1656  class SaddlePointProblem(object): Line 1774  class SaddlePointProblem(object):
1774    
1775     def trace(self,text):     def trace(self,text):
1776         """         """
1777         prints text if verbose has been set         Prints C{text} if verbose has been set.
1778    
1779         @param text: a text message         @param text: a text message
1780         @type text: C{str}         @type text: C{str}
# Line 1665  class SaddlePointProblem(object): Line 1783  class SaddlePointProblem(object):
1783    
1784     def solve_f(self,u,p,tol=1.e-8):     def solve_f(self,u,p,tol=1.e-8):
1785         """         """
1786         solves         Solves
1787    
1788         A_f du = f(u,p)         A_f du = f(u,p)
1789    
1790         with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.         with tolerance C{tol} and returns du. A_f is Jacobian of f with respect
1791           to u.
1792    
1793         @param u: current approximation of u         @param u: current approximation of u
1794         @type u: L{escript.Data}         @type u: L{escript.Data}
# Line 1679  class SaddlePointProblem(object): Line 1798  class SaddlePointProblem(object):
1798         @type tol: C{float}         @type tol: C{float}
1799         @return: increment du         @return: increment du
1800         @rtype: L{escript.Data}         @rtype: L{escript.Data}
1801         @note: this method has to be overwritten by a particular saddle point problem         @note: this method has to be overwritten by a particular saddle point
1802                  problem
1803         """         """
1804         pass         pass
1805    
1806     def solve_g(self,u,tol=1.e-8):     def solve_g(self,u,tol=1.e-8):
1807         """         """
1808         solves         Solves
1809    
1810         Q_g dp = g(u)         Q_g dp = g(u)
1811    
1812         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.         where Q_g is a preconditioner for A_g A_f^{-1} A_g with A_g is the
1813           Jacobian of g with respect to p.
1814    
1815         @param u: current approximation of u         @param u: current approximation of u
1816         @type u: L{escript.Data}         @type u: L{escript.Data}
# Line 1697  class SaddlePointProblem(object): Line 1818  class SaddlePointProblem(object):
1818         @type tol: C{float}         @type tol: C{float}
1819         @return: increment dp         @return: increment dp
1820         @rtype: L{escript.Data}         @rtype: L{escript.Data}
1821         @note: this method has to be overwritten by a particular saddle point problem         @note: this method has to be overwritten by a particular saddle point
1822                  problem
1823         """         """
1824         pass         pass
1825    
1826     def inner(self,p0,p1):     def inner(self,p0,p1):
1827         """         """
1828         inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)         Inner product of p0 and p1 approximating p. Typically this returns
1829           C{integrate(p0*p1)}.
1830         @return: inner product of p0 and p1         @return: inner product of p0 and p1
1831         @rtype: C{float}         @rtype: C{float}
1832         """         """
# Line 1712  class SaddlePointProblem(object): Line 1835  class SaddlePointProblem(object):
1835     subiter_max=3     subiter_max=3
1836     def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):     def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
1837          """          """
1838          runs the solver          Runs the solver.
1839    
1840          @param u0: initial guess for C{u}          @param u0: initial guess for C{u}
1841          @type u0: L{esys.escript.Data}          @type u0: L{esys.escript.Data}
# Line 1720  class SaddlePointProblem(object): Line 1843  class SaddlePointProblem(object):
1843          @type p0: L{esys.escript.Data}          @type p0: L{esys.escript.Data}
1844          @param tolerance: tolerance for relative error in C{u} and C{p}          @param tolerance: tolerance for relative error in C{u} and C{p}
1845          @type tolerance: positive C{float}          @type tolerance: positive C{float}
1846          @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}          @param tolerance_u: tolerance for relative error in C{u} if different
1847                                from C{tolerance}
1848          @type tolerance_u: positive C{float}          @type tolerance_u: positive C{float}
1849          @param iter_max: maximum number of iteration steps.          @param iter_max: maximum number of iteration steps
1850          @type iter_max: C{int}          @type iter_max: C{int}
1851          @param accepted_reduction: if the norm  g cannot be reduced by C{accepted_reduction} backtracking to adapt the          @param accepted_reduction: if the norm g cannot be reduced by
1852                                     relaxation factor. If C{accepted_reduction=None} no backtracking is used.                                     C{accepted_reduction} backtracking to adapt
1853                                       the relaxation factor. If
1854                                       C{accepted_reduction=None} no backtracking
1855                                       is used.
1856          @type accepted_reduction: positive C{float} or C{None}          @type accepted_reduction: positive C{float} or C{None}
1857          @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.          @param relaxation: initial relaxation factor. If C{relaxation==None},
1858                               the last relaxation factor is used.
1859          @type relaxation: C{float} or C{None}          @type relaxation: C{float} or C{None}
1860          """          """
1861          tol=1.e-2          tol=1.e-2
# Line 1826  class SaddlePointProblem(object): Line 1954  class SaddlePointProblem(object):
1954              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              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
1955          self.trace("convergence after %s steps."%self.iter)          self.trace("convergence after %s steps."%self.iter)
1956          return u,p          return u,p
1957    

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

  ViewVC Help
Powered by ViewVC 1.1.26