/[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 2445 by gross, Fri May 29 03:23:25 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 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{numarray.NumArray}, 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{numarray.NumArray}, 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{numarray.NumArray}, 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{numarray.NumArray}, 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{numarray.NumArray}, 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{numarray.NumArray}, 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{numarray.NumArray}, 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{numarray.NumArray}, 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=numarray.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{numarray.NumArray} or C{list} of C{numarray.NumArray}
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 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    
# Line 557  type like argument C{x}. Line 532  type like argument C{x}.
532    
533         rhat_dot_r = rhat_dot_r_new         rhat_dot_r = rhat_dot_r_new
534         if rhat_dot_r<0: raise NegativeNorm,"negative norm."         if rhat_dot_r<0: raise NegativeNorm,"negative norm."
535           if verbose: print "PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))
536     return x,r     if verbose: print "PCG: tolerance reached after %s steps."%iter
537       return x,r,math.sqrt(rhat_dot_r)
538    
539  class Defect(object):  class Defect(object):
540      """      """
541      defines a non-linear defect F(x) of a variable x      Defines a non-linear defect F(x) of a variable x.
542      """      """
543      def __init__(self):      def __init__(self):
544          """          """
545          initialize defect          Initializes defect.
546          """          """
547          self.setDerivativeIncrementLength()          self.setDerivativeIncrementLength()
548    
549      def bilinearform(self, x0, x1):      def bilinearform(self, x0, x1):
550          """          """
551          returns the inner product of x0 and x1          Returns the inner product of x0 and x1
552          @param x0: a value for x  
553          @param x1: a value for x          @param x0: value for x0
554            @param x1: value for x1
555          @return: the inner product of x0 and x1          @return: the inner product of x0 and x1
556          @rtype: C{float}          @rtype: C{float}
557          """          """
558          return 0          return 0
559          
560      def norm(self,x):      def norm(self,x):
561          """          """
562          the norm of argument C{x}          Returns the norm of argument C{x}.
563    
564          @param x: a value for x          @param x: a value
565          @return: norm of argument x          @return: norm of argument x
566          @rtype: C{float}          @rtype: C{float}
567          @note: by default C{sqrt(self.bilinearform(x,x)} is retrurned.          @note: by default C{sqrt(self.bilinearform(x,x)} is returned.
568          """          """
569          s=self.bilinearform(x,x)          s=self.bilinearform(x,x)
570          if s<0: raise NegativeNorm,"negative norm."          if s<0: raise NegativeNorm,"negative norm."
571          return math.sqrt(s)          return math.sqrt(s)
572    
   
573      def eval(self,x):      def eval(self,x):
574          """          """
575          returns the value F of a given x          Returns the value F of a given C{x}.
576    
577          @param x: value for which the defect C{F} is evalulated.          @param x: value for which the defect C{F} is evaluated
578          @return: value of the defect at C{x}          @return: value of the defect at C{x}
579          """          """
580          return 0          return 0
# Line 608  class Defect(object): Line 584  class Defect(object):
584    
585      def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):      def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):
586          """          """
587          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
588          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
589            direction of v with x as a starting point.
590    
591          @param inc: relative increment length          @param inc: relative increment length
592          @type inc: positive C{float}          @type inc: positive C{float}
# Line 619  class Defect(object): Line 596  class Defect(object):
596    
597      def getDerivativeIncrementLength(self):      def getDerivativeIncrementLength(self):
598          """          """
599          returns the relative increment length used to approximate the derivative of the defect          Returns the relative increment length used to approximate the
600            derivative of the defect.
601          @return: value of the defect at C{x}          @return: value of the defect at C{x}
602          @rtype: positive C{float}          @rtype: positive C{float}
603          """          """
# Line 627  class Defect(object): Line 605  class Defect(object):
605    
606      def derivative(self, F0, x0, v, v_is_normalised=True):      def derivative(self, F0, x0, v, v_is_normalised=True):
607          """          """
608          returns the directional derivative at x0 in the direction of v          Returns the directional derivative at C{x0} in the direction of C{v}.
609    
610          @param F0: value of this defect at x0          @param F0: value of this defect at x0
611          @param x0: value at which derivative is calculated.          @param x0: value at which derivative is calculated
612          @param v: direction          @param v: direction
613          @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
614                                    (self.norm(v)=0)
615          @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}
616          @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
617          maybe oepsnew verwritten to use exact evalution.                 used but this method maybe overwritten to use exact evaluation.
618          """          """
619          normx=self.norm(x0)          normx=self.norm(x0)
620          if normx>0:          if normx>0:
# Line 651  class Defect(object): Line 630  class Defect(object):
630          F1=self.eval(x0 + epsnew * v)          F1=self.eval(x0 + epsnew * v)
631          return (F1-F0)/epsnew          return (F1-F0)/epsnew
632    
633  ######################################      ######################################
634  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):
635     """     """
636     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
637       criterion:
638    
639     M{norm(F(x) <= atol + rtol * norm(F(x0)}     M{norm(F(x) <= atol + rtol * norm(F(x0)}
640      
641     where M{x0} is the initial guess.     where M{x0} is the initial guess.
642    
643     @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
644                      M{norm} used in the stopping criterion.
645     @type defect: L{Defect}     @type defect: L{Defect}
646     @param x: initial guess for the solution, C{x} is altered.     @param x: initial guess for the solution, C{x} is altered.
647     @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
648     @param iter_max: maximum number of iteration steps              C{numarray.NumArray}, L{Data}
649       @param iter_max: maximum number of iteration steps
650     @type iter_max: positive C{int}     @type iter_max: positive C{int}
651     @param sub_iter_max:     @param sub_iter_max: maximum number of inner iteration steps
652     @type sub_iter_max:     @type sub_iter_max: positive C{int}
653     @param atol: absolute tolerance for the solution     @param atol: absolute tolerance for the solution
654     @type atol: positive C{float}     @type atol: positive C{float}
655     @param rtol: relative tolerance for the solution     @param rtol: relative tolerance for the solution
656     @type rtol: positive C{float}     @type rtol: positive C{float}
657     @param gamma: tolerance safety factor for inner iteration     @param gamma: tolerance safety factor for inner iteration
658     @type gamma: positive C{float}, less than 1     @type gamma: positive C{float}, less than 1
659     @param sub_tol_max: upper bound for inner tolerance.     @param sub_tol_max: upper bound for inner tolerance
660     @type sub_tol_max: positive C{float}, less than 1     @type sub_tol_max: positive C{float}, less than 1
661     @return: an approximation of the solution with the desired accuracy     @return: an approximation of the solution with the desired accuracy
662     @rtype: same type as the initial guess.     @rtype: same type as the initial guess
663     """     """
664     lmaxit=iter_max     lmaxit=iter_max
665     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 679  def NewtonGMRES(defect, x, iter_max=100,
679     # main iteration loop     # main iteration loop
680     #     #
681     while not fnrm<=stop_tol:     while not fnrm<=stop_tol:
682              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
683              #              #
684          #   adjust sub_tol_          #   adjust sub_tol_
685          #          #
686              if iter > 1:              if iter > 1:
687             rat=fnrm/fnrmo             rat=fnrm/fnrmo
# Line 712  def NewtonGMRES(defect, x, iter_max=100, Line 694  def NewtonGMRES(defect, x, iter_max=100,
694              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown
695              #     if iter_restart -1 is returned as sub_iter              #     if iter_restart -1 is returned as sub_iter
696              #     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
697              #              #
698              #                #
699              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
700              try:              try:
701                 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 716  def NewtonGMRES(defect, x, iter_max=100,
716    
717  def __givapp(c,s,vin):  def __givapp(c,s,vin):
718      """      """
719      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
720        C{vin}
721    
722      @warning: C{vin} is altered.      @warning: C{vin} is altered.
723      """      """
724      vrot=vin      vrot=vin
725      if isinstance(c,float):      if isinstance(c,float):
726          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]]
727      else:      else:
728          for i in range(len(c)):          for i in range(len(c)):
729              w1=c[i]*vrot[i]-s[i]*vrot[i+1]              w1=c[i]*vrot[i]-s[i]*vrot[i+1]
730          w2=s[i]*vrot[i]+c[i]*vrot[i+1]          w2=s[i]*vrot[i]+c[i]*vrot[i+1]
731              vrot[i:i+2]=w1,w2              vrot[i]=w1
732                vrot[i+1]=w2
733      return vrot      return vrot
734    
735  def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):  def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
# Line 756  def __FDGMRES(F0, defect, x0, atol, iter Line 741  def __FDGMRES(F0, defect, x0, atol, iter
741    
742     rho=defect.norm(F0)     rho=defect.norm(F0)
743     if rho<=0.: return x0*0     if rho<=0.: return x0*0
744      
745     v.append(-F0/rho)     v.append(-F0/rho)
746     g[0]=rho     g[0]=rho
747     iter=0     iter=0
748     while rho > atol and iter<iter_restart-1:     while rho > atol and iter<iter_restart-1:
749            if iter  >= iter_max:
750      if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max              raise MaxIterReached,"maximum number of %s steps reached."%iter_max
751    
752          p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)          p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
753      v.append(p)          v.append(p)
754    
755      v_norm1=defect.norm(v[iter+1])          v_norm1=defect.norm(v[iter+1])
756    
757          # Modified Gram-Schmidt          # Modified Gram-Schmidt
758      for j in range(iter+1):          for j in range(iter+1):
759           h[j][iter]=defect.bilinearform(v[j],v[iter+1])                h[j,iter]=defect.bilinearform(v[j],v[iter+1])
760           v[iter+1]-=h[j][iter]*v[j]              v[iter+1]-=h[j,iter]*v[j]
761          
762      h[iter+1][iter]=defect.norm(v[iter+1])          h[iter+1,iter]=defect.norm(v[iter+1])
763      v_norm2=h[iter+1][iter]          v_norm2=h[iter+1,iter]
764    
765          # Reorthogonalize if needed          # Reorthogonalize if needed
766      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)
767          for j in range(iter+1):                for j in range(iter+1):
768             hr=defect.bilinearform(v[j],v[iter+1])                  hr=defect.bilinearform(v[j],v[iter+1])
769                 h[j][iter]=h[j][iter]+hr                  h[j,iter]=h[j,iter]+hr
770                 v[iter+1] -= hr*v[j]                  v[iter+1] -= hr*v[j]
771    
772          v_norm2=defect.norm(v[iter+1])              v_norm2=defect.norm(v[iter+1])
773          h[iter+1][iter]=v_norm2              h[iter+1,iter]=v_norm2
774          #   watch out for happy breakdown          #   watch out for happy breakdown
775          if not v_norm2 == 0:          if not v_norm2 == 0:
776                  v[iter+1]=v[iter+1]/h[iter+1][iter]              v[iter+1]=v[iter+1]/h[iter+1,iter]
777    
778          #   Form and store the information for the new Givens rotation          #   Form and store the information for the new Givens rotation
779      if iter > 0 :          if iter > 0 :
780          hhat=numarray.zeros(iter+1,numarray.Float64)              hhat=numarray.zeros(iter+1,numarray.Float64)
781          for i in range(iter+1) : hhat[i]=h[i][iter]              for i in range(iter+1) : hhat[i]=h[i,iter]
782          hhat=__givapp(c[0:iter],s[0:iter],hhat);              hhat=__givapp(c[0:iter],s[0:iter],hhat);
783              for i in range(iter+1) : h[i][iter]=hhat[i]              for i in range(iter+1) : h[i,iter]=hhat[i]
784    
785      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])
786    
787      if mu!=0 :          if mu!=0 :
788          c[iter]=h[iter][iter]/mu              c[iter]=h[iter,iter]/mu
789          s[iter]=-h[iter+1][iter]/mu              s[iter]=-h[iter+1,iter]/mu
790          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]
791          h[iter+1][iter]=0.0              h[iter+1,iter]=0.0
792          g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])              gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
793                g[iter]=gg[0]
794                g[iter+1]=gg[1]
795    
796          # Update the residual norm          # Update the residual norm
797          rho=abs(g[iter+1])          rho=abs(g[iter+1])
798      iter+=1          iter+=1
799    
800     # At this point either iter > iter_max or rho < tol.     # At this point either iter > iter_max or rho < tol.
801     # It's time to compute x and leave.             # It's time to compute x and leave.
802     if iter > 0 :     if iter > 0 :
803       y=numarray.zeros(iter,numarray.Float64)           y=numarray.zeros(iter,numarray.Float64)
804       y[iter-1] = g[iter-1] / h[iter-1][iter-1]       y[iter-1] = g[iter-1] / h[iter-1,iter-1]
805       if iter > 1 :         if iter > 1 :
806          i=iter-2            i=iter-2
807          while i>=0 :          while i>=0 :
808            y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]            y[i] = ( g[i] - numarray.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
809            i=i-1            i=i-1
810       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
811       for i in range(iter-1):       for i in range(iter-1):
812      xhat += v[i]*y[i]      xhat += v[i]*y[i]
813     else :     else :
814        xhat=v[0] * 0        xhat=v[0] * 0
815    
816     if iter<iter_restart-1:     if iter<iter_restart-1:
817        stopped=iter        stopped=iter
818     else:     else:
819        stopped=-1        stopped=-1
820    
821     return xhat,stopped     return xhat,stopped
822    
823  ##############################################  def GMRES(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100, iter_restart=20, verbose=False,P_R=None):
824  def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):     """
825  ################################################     Solver for
826    
827       M{Ax=b}
828    
829       with a general operator A (more details required!).
830       It uses the generalized minimum residual method (GMRES).
831    
832       The iteration is terminated if
833    
834       M{|r| <= atol+rtol*|r0|}
835    
836       where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
837    
838       M{|r| = sqrt( bilinearform(r,r))}
839    
840       @param r: initial residual M{r=b-Ax}. C{r} is altered.
841       @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
842       @param x: an initial guess for the solution
843       @type x: same like C{r}
844       @param Aprod: returns the value Ax
845       @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
846                    argument C{x}. The returned object needs to be of the same
847                    type like argument C{r}.
848       @param bilinearform: inner product C{<x,r>}
849       @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
850                           type like argument C{x} and C{r}. The returned value is
851                           a C{float}.
852       @param atol: absolute tolerance
853       @type atol: non-negative C{float}
854       @param rtol: relative tolerance
855       @type rtol: non-negative C{float}
856       @param iter_max: maximum number of iteration steps
857       @type iter_max: C{int}
858       @param iter_restart: in order to save memory the orthogonalization process
859                            is terminated after C{iter_restart} steps and the
860                            iteration is restarted.
861       @type iter_restart: C{int}
862       @return: the solution approximation and the corresponding residual
863       @rtype: C{tuple}
864       @warning: C{r} and C{x} are altered.
865       """
866     m=iter_restart     m=iter_restart
867       restarted=False
868     iter=0     iter=0
869     xc=x     if rtol>0:
870          r_dot_r = bilinearform(r, r)
871          if r_dot_r<0: raise NegativeNorm,"negative norm."
872          atol2=atol+rtol*math.sqrt(r_dot_r)
873          if verbose: print "GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)
874       else:
875          atol2=atol
876          if verbose: print "GMRES: absolute tolerance = %e"%atol2
877       if atol2<=0:
878          raise ValueError,"Non-positive tolarance."
879    
880     while True:     while True:
881        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
882        xc,stopped=__GMRESm(b*1, Aprod, Msolve, bilinearform, stoppingcriterium, x=xc*1, iter_max=iter_max-iter, iter_restart=m)        if restarted:
883             r2 = r-Aprod(x-x2)
884          else:
885             r2=1*r
886          x2=x*1.
887          x,stopped=_GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose,P_R=P_R)
888          iter+=iter_restart
889        if stopped: break        if stopped: break
890        iter+=iter_restart            if verbose: print "GMRES: restart."
891     return xc        restarted=True
892       if verbose: print "GMRES: tolerance has been reached."
893       return x
894    
895  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):
896     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)  
897    
898     if x==None:     h=numarray.zeros((iter_restart+1,iter_restart),numarray.Float64)
       x=0*b  
    else:  
       r=Msolve(b-Aprod(x))  
       r_dot_r = bilinearform(r, r)  
       if r_dot_r<0: raise NegativeNorm,"negative norm."  
     
    h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)  
899     c=numarray.zeros(iter_restart,numarray.Float64)     c=numarray.zeros(iter_restart,numarray.Float64)
900     s=numarray.zeros(iter_restart,numarray.Float64)     s=numarray.zeros(iter_restart,numarray.Float64)
901     g=numarray.zeros(iter_restart,numarray.Float64)     g=numarray.zeros(iter_restart+1,numarray.Float64)
902     v=[]     v=[]
903    
904       r_dot_r = bilinearform(r, r)
905       if r_dot_r<0: raise NegativeNorm,"negative norm."
906     rho=math.sqrt(r_dot_r)     rho=math.sqrt(r_dot_r)
907      
908     v.append(r/rho)     v.append(r/rho)
909     g[0]=rho     g[0]=rho
910    
911     while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):     if verbose: print "GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)
912       while not (rho<=atol or iter==iter_restart):
913    
914      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
915    
916      p=Msolve(Aprod(v[iter]))          if P_R!=None:
917                p=Aprod(P_R(v[iter]))
918            else:
919            p=Aprod(v[iter])
920      v.append(p)      v.append(p)
921    
922      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
923    
924    # Modified Gram-Schmidt
925        for j in range(iter+1):
926          h[j,iter]=bilinearform(v[j],v[iter+1])
927          v[iter+1]-=h[j,iter]*v[j]
928    
929  # Modified Gram-Schmidt      h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
930      for j in range(iter+1):      v_norm2=h[iter+1,iter]
       h[j][iter]=bilinearform(v[j],v[iter+1])    
       v[iter+1]-=h[j][iter]*v[j]  
         
     h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))  
     v_norm2=h[iter+1][iter]  
931    
932  # Reorthogonalize if needed  # Reorthogonalize if needed
933      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)
934       for j in range(iter+1):         for j in range(iter+1):
935          hr=bilinearform(v[j],v[iter+1])          hr=bilinearform(v[j],v[iter+1])
936              h[j][iter]=h[j][iter]+hr              h[j,iter]=h[j,iter]+hr
937              v[iter+1] -= hr*v[j]              v[iter+1] -= hr*v[j]
938    
939       v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
940       h[iter+1][iter]=v_norm2       h[iter+1,iter]=v_norm2
941    
942  #   watch out for happy breakdown  #   watch out for happy breakdown
943          if not v_norm2 == 0:          if not v_norm2 == 0:
944           v[iter+1]=v[iter+1]/h[iter+1][iter]           v[iter+1]=v[iter+1]/h[iter+1,iter]
945    
946  #   Form and store the information for the new Givens rotation  #   Form and store the information for the new Givens rotation
947      if iter > 0 :      if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])
948          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])  
949    
950      if mu!=0 :      if mu!=0 :
951          c[iter]=h[iter][iter]/mu          c[iter]=h[iter,iter]/mu
952          s[iter]=-h[iter+1][iter]/mu          s[iter]=-h[iter+1,iter]/mu
953          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]
954          h[iter+1][iter]=0.0          h[iter+1,iter]=0.0
955          g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])                  gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
956                    g[iter]=gg[0]
957                    g[iter+1]=gg[1]
958  # Update the residual norm  # Update the residual norm
959                  
960          rho=abs(g[iter+1])          rho=abs(g[iter+1])
961            if verbose: print "GMRES: iteration step %s: residual %e"%(iter,rho)
962      iter+=1      iter+=1
963    
964  # At this point either iter > iter_max or rho < tol.  # At this point either iter > iter_max or rho < tol.
965  # It's time to compute x and leave.          # It's time to compute x and leave.
966    
967     if iter > 0 :     if verbose: print "GMRES: iteration stopped after %s step."%iter
968       y=numarray.zeros(iter,numarray.Float64)         if iter > 0 :
969       y[iter-1] = g[iter-1] / h[iter-1][iter-1]       y=numarray.zeros(iter,numarray.Float64)
970       if iter > 1 :         y[iter-1] = g[iter-1] / h[iter-1,iter-1]
971          i=iter-2         if iter > 1 :
972            i=iter-2
973          while i>=0 :          while i>=0 :
974            y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]            y[i] = ( g[i] - numarray.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
975            i=i-1            i=i-1
976       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
977       for i in range(iter-1):       for i in range(iter-1):
978      xhat += v[i]*y[i]      xhat += v[i]*y[i]
979     else : xhat=v[0]     else:
980         xhat=v[0] * 0
981     x += xhat     if P_R!=None:
982     if iter<iter_restart-1:        x += P_R(xhat)
983        stopped=True     else:
984     else:        x += xhat
985       if iter<iter_restart-1:
986          stopped=True
987       else:
988        stopped=False        stopped=False
989    
990     return x,stopped     return x,stopped
991    
992  #################################################  def MINRES(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
993  def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):      """
994  #################################################      Solver for
995      #  
996      #  minres solves the system of linear equations Ax = b      M{Ax=b}
997      #  where A is a symmetric matrix (possibly indefinite or singular)  
998      #  and b is a given vector.      with a symmetric and positive definite operator A (more details required!).
999      #        It uses the minimum residual method (MINRES) with preconditioner M
1000      #  "A" may be a dense or sparse matrix (preferably sparse!)      providing an approximation of A.
1001      #  or the name of a function such that  
1002      #               y = A(x)      The iteration is terminated if
1003      #  returns the product y = Ax for any given vector x.  
1004      #      M{|r| <= atol+rtol*|r0|}
1005      #  "M" defines a positive-definite preconditioner M = C C'.  
1006      #  "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
1007      #  or the name of a function such that  
1008      #  solves the system My = x for any given vector x.      M{|r| = sqrt( bilinearform(Msolve(r),r))}
1009      #  
1010      #      For details on the preconditioned conjugate gradient method see the book:
1011        
1012        I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
1013        T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
1014        C. Romine, and H. van der Vorst}.
1015    
1016        @param r: initial residual M{r=b-Ax}. C{r} is altered.
1017        @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1018        @param x: an initial guess for the solution
1019        @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1020        @param Aprod: returns the value Ax
1021        @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
1022                     argument C{x}. The returned object needs to be of the same
1023                     type like argument C{r}.
1024        @param Msolve: solves Mx=r
1025        @type Msolve: function C{Msolve(r)} where C{r} is of the same type like
1026                      argument C{r}. The returned object needs to be of the same
1027                      type like argument C{x}.
1028        @param bilinearform: inner product C{<x,r>}
1029        @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
1030                            type like argument C{x} and C{r} is. The returned value
1031                            is a C{float}.
1032        @param atol: absolute tolerance
1033        @type atol: non-negative C{float}
1034        @param rtol: relative tolerance
1035        @type rtol: non-negative C{float}
1036        @param iter_max: maximum number of iteration steps
1037        @type iter_max: C{int}
1038        @return: the solution approximation and the corresponding residual
1039        @rtype: C{tuple}
1040        @warning: C{r} and C{x} are altered.
1041        """
1042      #------------------------------------------------------------------      #------------------------------------------------------------------
1043      # Set up y and v for the first Lanczos vector v1.      # Set up y and v for the first Lanczos vector v1.
1044      # y  =  beta1 P' v1,  where  P = C**(-1).      # y  =  beta1 P' v1,  where  P = C**(-1).
1045      # v is really P' v1.      # v is really P' v1.
1046      #------------------------------------------------------------------      #------------------------------------------------------------------
1047      if x==None:      r1    = r
1048        x=0*b      y = Msolve(r)
1049      else:      beta1 = bilinearform(y,r)
       b += (-1)*Aprod(x)  
1050    
     r1    = b  
     y = Msolve(b)  
     beta1 = bilinearform(y,b)  
   
1051      if beta1< 0: raise NegativeNorm,"negative norm."      if beta1< 0: raise NegativeNorm,"negative norm."
1052    
1053      #  If b = 0 exactly, stop with x = 0.      #  If r = 0 exactly, stop with x
1054      if beta1==0: return x*0.      if beta1==0: return x
1055    
1056      if beta1> 0:      if beta1> 0: beta1  = math.sqrt(beta1)
       beta1  = math.sqrt(beta1)        
1057    
1058      #------------------------------------------------------------------      #------------------------------------------------------------------
1059      # Initialize quantities.      # Initialize quantities.
# Line 1008  def MINRES(b, Aprod, Msolve, bilinearfor Line 1073  def MINRES(b, Aprod, Msolve, bilinearfor
1073      ynorm2 = 0      ynorm2 = 0
1074      cs     = -1      cs     = -1
1075      sn     = 0      sn     = 0
1076      w      = b*0.      w      = r*0.
1077      w2     = b*0.      w2     = r*0.
1078      r2     = r1      r2     = r1
1079      eps    = 0.0001      eps    = 0.0001
1080    
1081      #---------------------------------------------------------------------      #---------------------------------------------------------------------
1082      # Main iteration loop.      # Main iteration loop.
1083      # --------------------------------------------------------------------      # --------------------------------------------------------------------
1084      while not stoppingcriterium(rnorm,Anorm*ynorm,'MINRES'):    #  checks ||r|| < (||A|| ||x||) * TOL      while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL
1085    
1086      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
1087          iter    = iter  +  1          iter    = iter  +  1
# Line 1035  def MINRES(b, Aprod, Msolve, bilinearfor Line 1100  def MINRES(b, Aprod, Msolve, bilinearfor
1100          #-----------------------------------------------------------------          #-----------------------------------------------------------------
1101          s = 1/beta                 # Normalize previous vector (in y).          s = 1/beta                 # Normalize previous vector (in y).
1102          v = s*y                    # v = vk if P = I          v = s*y                    # v = vk if P = I
1103        
1104          y      = Aprod(v)          y      = Aprod(v)
1105        
1106          if iter >= 2:          if iter >= 2:
1107            y = y - (beta/oldb)*r1            y = y - (beta/oldb)*r1
1108    
1109          alfa   = bilinearform(v,y)              # alphak          alfa   = bilinearform(v,y)              # alphak
1110          y      += (- alfa/beta)*r2          y      += (- alfa/beta)*r2
1111          r1     = r2          r1     = r2
1112          r2     = y          r2     = y
1113          y = Msolve(r2)          y = Msolve(r2)
# Line 1052  def MINRES(b, Aprod, Msolve, bilinearfor Line 1117  def MINRES(b, Aprod, Msolve, bilinearfor
1117    
1118          beta   = math.sqrt( beta )          beta   = math.sqrt( beta )
1119          tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta          tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
1120            
1121          if iter==1:                 # Initialize a few things.          if iter==1:                 # Initialize a few things.
1122            gmax   = abs( alfa )      # alpha1            gmax   = abs( alfa )      # alpha1
1123            gmin   = gmax             # alpha1            gmin   = gmax             # alpha1
# Line 1060  def MINRES(b, Aprod, Msolve, bilinearfor Line 1125  def MINRES(b, Aprod, Msolve, bilinearfor
1125          # Apply previous rotation Qk-1 to get          # Apply previous rotation Qk-1 to get
1126          #   [deltak epslnk+1] = [cs  sn][dbark    0   ]          #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
1127          #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].          #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
1128        
1129          oldeps = epsln          oldeps = epsln
1130          delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak          delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
1131          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 1135  def MINRES(b, Aprod, Msolve, bilinearfor
1135          # Compute the next plane rotation Qk          # Compute the next plane rotation Qk
1136    
1137          gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak          gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
1138          gamma  = max(gamma,eps)          gamma  = max(gamma,eps)
1139          cs     = gbar / gamma             # ck          cs     = gbar / gamma             # ck
1140          sn     = beta / gamma             # sk          sn     = beta / gamma             # sk
1141          phi    = cs * phibar              # phik          phi    = cs * phibar              # phik
# Line 1078  def MINRES(b, Aprod, Msolve, bilinearfor Line 1143  def MINRES(b, Aprod, Msolve, bilinearfor
1143    
1144          # Update  x.          # Update  x.
1145    
1146          denom = 1/gamma          denom = 1/gamma
1147          w1    = w2          w1    = w2
1148          w2    = w          w2    = w
1149          w     = (v - oldeps*w1 - delta*w2) * denom          w     = (v - oldeps*w1 - delta*w2) * denom
1150          x     +=  phi*w          x     +=  phi*w
1151    
# Line 1095  def MINRES(b, Aprod, Msolve, bilinearfor Line 1160  def MINRES(b, Aprod, Msolve, bilinearfor
1160    
1161          # Estimate various norms and test for convergence.          # Estimate various norms and test for convergence.
1162    
1163          Anorm  = math.sqrt( tnorm2 )          Anorm  = math.sqrt( tnorm2 )
1164          ynorm  = math.sqrt( ynorm2 )          ynorm  = math.sqrt( ynorm2 )
1165    
1166          rnorm  = phibar          rnorm  = phibar
1167    
1168      return x      return x
1169    
1170  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):
1171      """
1172      Solver for
1173    
1174  # TFQMR solver for linear systems    M{Ax=b}
1175  #  
1176  #    with a general operator A (more details required!).
1177  # initialization    It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).
 #  
   errtol = math.sqrt(bilinearform(b,b))  
   norm_b=errtol  
   kmax  = iter_max  
   error = []  
   
   if math.sqrt(bilinearform(x,x)) != 0.0:  
     r = b - Aprod(x)  
   else:  
     r = b  
1178    
1179    r=Msolve(r)    The iteration is terminated if
1180    
1181      M{|r| <= atol+rtol*|r0|}
1182    
1183      where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
1184    
1185      M{|r| = sqrt( bilinearform(r,r))}
1186    
1187      @param r: initial residual M{r=b-Ax}. C{r} is altered.
1188      @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1189      @param x: an initial guess for the solution
1190      @type x: same like C{r}
1191      @param Aprod: returns the value Ax
1192      @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
1193                   argument C{x}. The returned object needs to be of the same type
1194                   like argument C{r}.
1195      @param bilinearform: inner product C{<x,r>}
1196      @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
1197                          type like argument C{x} and C{r}. The returned value is
1198                          a C{float}.
1199      @param atol: absolute tolerance
1200      @type atol: non-negative C{float}
1201      @param rtol: relative tolerance
1202      @type rtol: non-negative C{float}
1203      @param iter_max: maximum number of iteration steps
1204      @type iter_max: C{int}
1205      @rtype: C{tuple}
1206      @warning: C{r} and C{x} are altered.
1207      """
1208    u1=0    u1=0
1209    u2=0    u2=0
1210    y1=0    y1=0
1211    y2=0    y2=0
1212    
1213    w = r    w = r
1214    y1 = r    y1 = r
1215    iter = 0    iter = 0
1216    d = 0    d = 0
1217        v = Aprod(y1)
   v = Msolve(Aprod(y1))  
1218    u1 = v    u1 = v
1219      
1220    theta = 0.0;    theta = 0.0;
1221    eta = 0.0;    eta = 0.0;
1222    tau = math.sqrt(bilinearform(r,r))    rho=bilinearform(r,r)
1223    error = [ error, tau ]    if rho < 0: raise NegativeNorm,"negative norm."
1224    rho = tau * tau    tau = math.sqrt(rho)
1225    m=1    norm_r0=tau
1226  #    while tau>atol+rtol*norm_r0:
 #  TFQMR iteration  
 #  
 #  while ( iter < kmax-1 ):  
     
   while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):  
1227      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
1228    
1229      sigma = bilinearform(r,v)      sigma = bilinearform(r,v)
1230        if sigma == 0.0: raise IterationBreakDown,'TFQMR breakdown, sigma=0'
     if ( sigma == 0.0 ):  
       raise 'TFQMR breakdown, sigma=0'  
       
1231    
1232      alpha = rho / sigma      alpha = rho / sigma
1233    
# Line 1162  def TFQMR(b, Aprod, Msolve, bilinearform Line 1237  def TFQMR(b, Aprod, Msolve, bilinearform
1237  #  #
1238        if ( j == 1 ):        if ( j == 1 ):
1239          y2 = y1 - alpha * v          y2 = y1 - alpha * v
1240          u2 = Msolve(Aprod(y2))          u2 = Aprod(y2)
1241    
1242        m = 2 * (iter+1) - 2 + (j+1)        m = 2 * (iter+1) - 2 + (j+1)
1243        if j==0:        if j==0:
1244           w = w - alpha * u1           w = w - alpha * u1
1245           d = y1 + ( theta * theta * eta / alpha ) * d           d = y1 + ( theta * theta * eta / alpha ) * d
1246        if j==1:        if j==1:
# Line 1180  def TFQMR(b, Aprod, Msolve, bilinearform Line 1255  def TFQMR(b, Aprod, Msolve, bilinearform
1255  #  #
1256  #  Try to terminate the iteration at each pass through the loop  #  Try to terminate the iteration at each pass through the loop
1257  #  #
1258       # 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'  
       
1259    
1260      rhon = bilinearform(r,w)      rhon = bilinearform(r,w)
1261      beta = rhon / rho;      beta = rhon / rho;
1262      rho = rhon;      rho = rhon;
1263      y1 = w + beta * y2;      y1 = w + beta * y2;
1264      u1 = Msolve(Aprod(y1))      u1 = Aprod(y1)
1265      v = u1 + beta * ( u2 + beta * v )      v = u1 + beta * ( u2 + beta * v )
1266      error = [ error, tau ]  
1267      total_iters = iter      iter += 1
       
     iter = iter + 1  
1268    
1269    return x    return x
1270    
# Line 1208  def TFQMR(b, Aprod, Msolve, bilinearform Line 1273  def TFQMR(b, Aprod, Msolve, bilinearform
1273    
1274  class ArithmeticTuple(object):  class ArithmeticTuple(object):
1275     """     """
1276     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
1277       ArithmeticTuple and C{a} is a float.
1278    
1279     example of usage:     Example of usage::
1280    
1281     from esys.escript import Data         from esys.escript import Data
1282     from numarray import array         from numarray import array
1283     a=Data(...)         a=Data(...)
1284     b=array([1.,4.])         b=array([1.,4.])
1285     x=ArithmeticTuple(a,b)         x=ArithmeticTuple(a,b)
1286     y=5.*x         y=5.*x
1287    
1288     """     """
1289     def __init__(self,*args):     def __init__(self,*args):
1290         """         """
1291         initialize object with elements args.         Initializes object with elements C{args}.
1292    
1293         @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
1294                        scaling (x=a*y)
1295         """         """
1296         self.__items=list(args)         self.__items=list(args)
1297    
1298     def __len__(self):     def __len__(self):
1299         """         """
1300         number of items         Returns the number of items.
1301    
1302         @return: number of items         @return: number of items
1303         @rtype: C{int}         @rtype: C{int}
# Line 1239  class ArithmeticTuple(object): Line 1306  class ArithmeticTuple(object):
1306    
1307     def __getitem__(self,index):     def __getitem__(self,index):
1308         """         """
1309         get an item         Returns item at specified position.
1310    
1311         @param index: item to be returned         @param index: index of item to be returned
1312         @type index: C{int}         @type index: C{int}
1313         @return: item with index C{index}         @return: item with index C{index}
1314         """         """
# Line 1249  class ArithmeticTuple(object): Line 1316  class ArithmeticTuple(object):
1316    
1317     def __mul__(self,other):     def __mul__(self,other):
1318         """         """
1319         scaling from the right         Scales by C{other} from the right.
1320    
1321         @param other: scaling factor         @param other: scaling factor
1322         @type other: C{float}         @type other: C{float}
# Line 1257  class ArithmeticTuple(object): Line 1324  class ArithmeticTuple(object):
1324         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1325         """         """
1326         out=[]         out=[]
1327         try:           try:
1328             l=len(other)             l=len(other)
1329             if l!=len(self):             if l!=len(self):
1330                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1331             for i in range(l): out.append(self[i]*other[i])             for i in range(l): out.append(self[i]*other[i])
1332         except TypeError:         except TypeError:
1333         for i in range(len(self)): out.append(self[i]*other)             for i in range(len(self)): out.append(self[i]*other)
1334         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1335    
1336     def __rmul__(self,other):     def __rmul__(self,other):
1337         """         """
1338         scaling from the left         Scales by C{other} from the left.
1339    
1340         @param other: scaling factor         @param other: scaling factor
1341         @type other: C{float}         @type other: C{float}
# Line 1276  class ArithmeticTuple(object): Line 1343  class ArithmeticTuple(object):
1343         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1344         """         """
1345         out=[]         out=[]
1346         try:           try:
1347             l=len(other)             l=len(other)
1348             if l!=len(self):             if l!=len(self):
1349                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1350             for i in range(l): out.append(other[i]*self[i])             for i in range(l): out.append(other[i]*self[i])
1351         except TypeError:         except TypeError:
1352         for i in range(len(self)): out.append(other*self[i])             for i in range(len(self)): out.append(other*self[i])
1353         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1354    
1355     def __div__(self,other):     def __div__(self,other):
1356         """         """
1357         dividing from the right         Scales by (1/C{other}) from the right.
1358    
1359         @param other: scaling factor         @param other: scaling factor
1360         @type other: C{float}         @type other: C{float}
# Line 1298  class ArithmeticTuple(object): Line 1365  class ArithmeticTuple(object):
1365    
1366     def __rdiv__(self,other):     def __rdiv__(self,other):
1367         """         """
1368         dividing from the left         Scales by (1/C{other}) from the left.
1369    
1370         @param other: scaling factor         @param other: scaling factor
1371         @type other: C{float}         @type other: C{float}
# Line 1306  class ArithmeticTuple(object): Line 1373  class ArithmeticTuple(object):
1373         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1374         """         """
1375         out=[]         out=[]
1376         try:           try:
1377             l=len(other)             l=len(other)
1378             if l!=len(self):             if l!=len(self):
1379                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1380             for i in range(l): out.append(other[i]/self[i])             for i in range(l): out.append(other[i]/self[i])
1381         except TypeError:         except TypeError:
1382         for i in range(len(self)): out.append(other/self[i])             for i in range(len(self)): out.append(other/self[i])
1383         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1384      
1385     def __iadd__(self,other):     def __iadd__(self,other):
1386         """         """
1387         in-place add of other to self         Inplace addition of C{other} to self.
1388    
1389         @param other: increment         @param other: increment
1390         @type other: C{ArithmeticTuple}         @type other: C{ArithmeticTuple}
1391         """         """
1392         if len(self) != len(other):         if len(self) != len(other):
1393             raise ValueError,"tuple length must match."             raise ValueError,"tuple lengths must match."
1394         for i in range(len(self)):         for i in range(len(self)):
1395             self.__items[i]+=other[i]             self.__items[i]+=other[i]
1396         return self         return self
1397    
1398     def __add__(self,other):     def __add__(self,other):
1399         """         """
1400         add other to self         Adds C{other} to self.
1401    
1402         @param other: increment         @param other: increment
1403         @type other: C{ArithmeticTuple}         @type other: C{ArithmeticTuple}
1404         """         """
1405         out=[]         out=[]
1406         try:           try:
1407             l=len(other)             l=len(other)
1408             if l!=len(self):             if l!=len(self):
1409                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1410             for i in range(l): out.append(self[i]+other[i])             for i in range(l): out.append(self[i]+other[i])
1411         except TypeError:         except TypeError:
1412         for i in range(len(self)): out.append(self[i]+other)             for i in range(len(self)): out.append(self[i]+other)
1413         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1414    
1415     def __sub__(self,other):     def __sub__(self,other):
1416         """         """
1417         subtract other from self         Subtracts C{other} from self.
1418    
1419         @param other: increment         @param other: decrement
1420         @type other: C{ArithmeticTuple}         @type other: C{ArithmeticTuple}
1421         """         """
1422         out=[]         out=[]
1423         try:           try:
1424             l=len(other)             l=len(other)
1425             if l!=len(self):             if l!=len(self):
1426                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1427             for i in range(l): out.append(self[i]-other[i])             for i in range(l): out.append(self[i]-other[i])
1428         except TypeError:         except TypeError:
1429         for i in range(len(self)): out.append(self[i]-other)             for i in range(len(self)): out.append(self[i]-other)
1430         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1431      
1432     def __isub__(self,other):     def __isub__(self,other):
1433         """         """
1434         subtract other from self         Inplace subtraction of C{other} from self.
1435    
1436         @param other: increment         @param other: decrement
1437         @type other: C{ArithmeticTuple}         @type other: C{ArithmeticTuple}
1438         """         """
1439         if len(self) != len(other):         if len(self) != len(other):
# Line 1377  class ArithmeticTuple(object): Line 1444  class ArithmeticTuple(object):
1444    
1445     def __neg__(self):     def __neg__(self):
1446         """         """
1447         negate         Negates values.
   
1448         """         """
1449         out=[]         out=[]
1450         for i in range(len(self)):         for i in range(len(self)):
# Line 1388  class ArithmeticTuple(object): Line 1454  class ArithmeticTuple(object):
1454    
1455  class HomogeneousSaddlePointProblem(object):  class HomogeneousSaddlePointProblem(object):
1456        """        """
1457        This provides a framwork for solving linear homogeneous saddle point problem of the form        This class provides a framework for solving linear homogeneous saddle
1458          point problems of the form::
              Av+B^*p=f  
              Bv    =0  
1459    
1460        for the unknowns v and p and given operators A and B and given right hand side f.            M{Av+B^*p=f}
1461        B^* is the adjoint operator of B is the given inner product.            M{Bv     =0}
1462    
1463          for the unknowns M{v} and M{p} and given operators M{A} and M{B} and
1464          given right hand side M{f}. M{B^*} is the adjoint operator of M{B}.
1465        """        """
1466        def __init__(self,**kwargs):        def __init__(self,**kwargs):
1467          self.setTolerance()          self.setTolerance()
1468          self.setToleranceReductionFactor()          self.setAbsoluteTolerance()
1469            self.setSubProblemTolerance()
1470    
1471    
1472          #=============================================================
1473        def initialize(self):        def initialize(self):
1474          """          """
1475          initialize the problem (overwrite)          Initializes the problem (overwrite).
1476          """          """
1477          pass          pass
1478        def B(self,v):  
1479          def inner_pBv(self,p,Bv):
1480           """           """
1481           returns Bv (overwrite)           Returns inner product of element p and Bv (overwrite).
          @rtype: equal to the type of p  
1482    
1483           @note: boundary conditions on p should be zero!           @param p: a pressure increment
1484             @param v: a residual
1485             @return: inner product of element p and Bv
1486             @rtype: C{float}
1487             @note: used if PCG is applied.
1488           """           """
1489           pass           raise NotImplementedError,"no inner product for p and Bv implemented."
1490    
1491        def inner(self,p0,p1):        def inner_p(self,p0,p1):
1492           """           """
1493           returns inner product of two element p0 and p1  (overwrite)           Returns inner product of p0 and p1 (overwrite).
1494            
1495           @type p0: equal to the type of p           @param p0: a pressure
1496           @type p1: equal to the type of p           @param p1: a pressure
1497             @return: inner product of p0 and p1
1498           @rtype: C{float}           @rtype: C{float}
1499             """
1500             raise NotImplementedError,"no inner product for p implemented."
1501      
1502          def norm_v(self,v):
1503             """
1504             Returns the norm of v (overwrite).
1505    
1506           @rtype: equal to the type of p           @param v: a velovity
1507             @return: norm of v
1508             @rtype: non-negative C{float}
1509             """
1510             raise NotImplementedError,"no norm of v implemented."
1511          def getV(self, p, v0):
1512           """           """
1513           pass           return the value for v for a given p (overwrite)
1514    
1515             @param p: a pressure
1516             @param v0: a initial guess for the value v to return.
1517             @return: v given as M{v= A^{-1} (f-B^*p)}
1518             """
1519             raise NotImplementedError,"no v calculation implemented."
1520    
1521            
1522          def Bv(self,v):
1523            """
1524            Returns Bv (overwrite).
1525    
1526        def solve_A(self,u,p):          @rtype: equal to the type of p
1527            @note: boundary conditions on p should be zero!
1528            """
1529            raise NotImplementedError, "no operator B implemented."
1530    
1531          def norm_Bv(self,Bv):
1532            """
1533            Returns the norm of Bv (overwrite).
1534    
1535            @rtype: equal to the type of p
1536            @note: boundary conditions on p should be zero!
1537            """
1538            raise NotImplementedError, "no norm of Bv implemented."
1539    
1540          def solve_AinvBt(self,p):
1541           """           """
1542           solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
1543             (overwrite).
1544    
1545           @rtype: equal to the type of v           @param p: a pressure increment
1546             @return: the solution of M{Av=B^*p}
1547           @note: boundary conditions on v should be zero!           @note: boundary conditions on v should be zero!
1548           """           """
1549           pass           raise NotImplementedError,"no operator A implemented."
1550    
1551        def solve_prec(self,p):        def solve_prec(self,Bv):
1552           """           """
1553           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
1554             L{self.getSubProblemTolerance()} (overwrite).
1555    
1556           @rtype: equal to the type of p           @rtype: equal to the type of p
1557             @note: boundary conditions on p should be zero!
1558           """           """
1559           pass           raise NotImplementedError,"no preconditioner for Schur complement implemented."
1560          #=============================================================
1561          def __Aprod_PCG(self,p):
1562              dv=self.solve_AinvBt(p)
1563              return ArithmeticTuple(dv,self.Bv(dv))
1564    
1565        def stoppingcriterium(self,Bv,v,p):        def __inner_PCG(self,p,r):
1566           """           return self.inner_pBv(p,r[1])
          returns a True if iteration is terminated. (overwrite)  
1567    
1568           @rtype: C{bool}        def __Msolve_PCG(self,r):
1569           """            return self.solve_prec(r[1])
1570           pass        #=============================================================
1571                # rename solve_prec and change argument v to Bv
1572        def __inner(self,p,r):  # chnage the argument of inner_pBv to v->Bv
1573           return self.inner(p,r[1])  # add Bv
1574    # inner p still needed?
1575    # change norm_Bv argument to Bv
1576          def __Aprod_GMRES(self,p):
1577              return self.solve_prec(self.Bv(self.solve_AinvBt(p)))
1578          def __inner_GMRES(self,p0,p1):
1579             return self.inner_p(p0,p1)
1580          #=============================================================
1581          def norm_p(self,p):
1582              """
1583              calculates the norm of C{p}
1584              
1585              @param p: a pressure
1586              @return: the norm of C{p} using the inner product for pressure
1587              @rtype: C{float}
1588              """
1589              f=self.inner_p(p,p)
1590              if f<0: raise ValueError,"negative pressure norm."
1591              return math.sqrt(f)
1592              
1593    
1594        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):
1595           return self.inner(p1,p2)           """
1596                   Solves the saddle point problem using initial guesses v and p.
       def __inner_a(self,a1,a2):  
          return self.inner_a(a1,a2)  
1597    
1598        def __inner_a1(self,a1,a2):           @param v: initial guess for velocity
1599           return self.inner(a1[1],a2[1])           @param p: initial guess for pressure
1600             @type v: L{Data}
1601             @type p: L{Data}
1602             @param usePCG: indicates the usage of the PCG rather than GMRES scheme.
1603             @param max_iter: maximum number of iteration steps per correction
1604                              attempt
1605             @param verbose: if True, shows information on the progress of the
1606                             saddlepoint problem solver.
1607             @param show_details: if True, shows details of the sub problem solver
1608             @param iter_restart: restart the iteration after C{iter_restart} steps
1609                                  (only used if useUzaw=False)
1610             @type usePCG: C{bool}
1611             @type max_iter: C{int}
1612             @type verbose: C{bool}
1613             @type show_details: C{bool}
1614             @type iter_restart: C{int}
1615             @rtype: C{tuple} of L{Data} objects
1616             """
1617             self.verbose=verbose
1618             self.show_details=show_details and self.verbose
1619             rtol=self.getTolerance()
1620             atol=self.getAbsoluteTolerance()
1621             correction_step=0
1622             converged=False
1623             while not converged:
1624                  # calculate velocity for current pressure:
1625                  v=self.getV(p,v)
1626                  Bv=self.Bv(v)
1627                  norm_v=self.norm_v(v)
1628                  norm_Bv=self.norm_Bv(Bv)
1629                  ATOL=norm_v*rtol+atol
1630                  if self.verbose: print "HomogeneousSaddlePointProblem: norm v= %e, norm_Bv= %e, tolerance = %e."%(norm_v, norm_Bv,ATOL)
1631                  if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."
1632                  if norm_Bv <= ATOL:
1633                     converged=True
1634                  else:
1635                     correction_step+=1
1636                     if correction_step>max_correction_steps:
1637                          raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1638                     dp=self.solve_prec(Bv)
1639                     if usePCG:
1640                       norm2=self.inner_pBv(dp,Bv)
1641                       if norm2<0: raise ValueError,"negative PCG norm."
1642                       norm2=math.sqrt(norm2)
1643                     else:
1644                       norm2=self.norm_p(dp)
1645                     ATOL_ITER=ATOL/norm_Bv*norm2*0.5
1646                     if self.verbose: print "HomogeneousSaddlePointProblem: tolerance for solver: %e"%ATOL_ITER
1647                     if usePCG:
1648                           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)
1649                     else:
1650                           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)
1651             if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached."
1652         return v,p
1653    
1654        def __stoppingcriterium(self,norm_r,r,p):        #========================================================================
1655            return self.stoppingcriterium(r[1],r[0],p)        def setTolerance(self,tolerance=1.e-4):
1656             """
1657             Sets the relative tolerance for (v,p).
1658    
1659        def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):           @param tolerance: tolerance to be used
1660            return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)           @type tolerance: non-negative C{float}
1661             """
1662             if tolerance<0:
1663                 raise ValueError,"tolerance must be positive."
1664             self.__rtol=tolerance
1665             self.setSubProblemTolerance()
1666    
       def setTolerance(self,tolerance=1.e-8):  
               self.__tol=tolerance  
1667        def getTolerance(self):        def getTolerance(self):
1668                return self.__tol           """
1669        def setToleranceReductionFactor(self,reduction=0.01):           Returns the relative tolerance.
               self.__reduction=reduction  
       def getSubProblemTolerance(self):  
               return self.__reduction*self.getTolerance()  
1670    
1671        def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):           @return: relative tolerance
1672                """           @rtype: C{float}
1673                solves the saddle point problem using initial guesses v and p.           """
1674             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)  
1675    
1676                # print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)        def setAbsoluteTolerance(self,tolerance=0.):
1677             """
1678             Sets the absolute tolerance.
1679    
1680            return u,p           @param tolerance: tolerance to be used
1681             @type tolerance: non-negative C{float}
1682             """
1683             if tolerance<0:
1684                 raise ValueError,"tolerance must be non-negative."
1685             self.__atol=tolerance
1686    
1687        def __Msolve(self,r):        def getAbsoluteTolerance(self):
1688            return self.solve_prec(r[1])           """
1689             Returns the absolute tolerance.
1690    
1691        def __Msolve2(self,r):           @return: absolute tolerance
1692            return self.solve_prec(r*1)           @rtype: C{float}
1693             """
1694             return self.__atol
1695    
1696        def __Mempty(self,r):        def setSubProblemTolerance(self,rtol=None):
1697            return r           """
1698             Sets the relative tolerance to solve the subproblem(s).
1699    
1700             @param rtol: relative tolerence
1701             @type rtol: positive C{float}
1702             """
1703             if rtol == None:
1704                  rtol=max(200.*util.EPSILON,self.getTolerance()**2)
1705             if rtol<=0:
1706                 raise ValueError,"tolerance must be positive."
1707             self.__sub_tol=rtol
1708    
1709        def __Aprod(self,p):        def getSubProblemTolerance(self):
1710            # return BA^-1B*p           """
1711            #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)  
1712    
1713             @return: subproblem reduction factor
1714             @rtype: C{float}
1715             """
1716             return self.__sub_tol
1717    
1718  def MaskFromBoundaryTag(domain,*tags):  def MaskFromBoundaryTag(domain,*tags):
1719     """     """
1720     creates a mask on the Solution(domain) function space which one for samples     Creates a mask on the Solution(domain) function space where the value is
1721     that touch the boundary tagged by tags.     one for samples that touch the boundary tagged by tags.
1722    
1723     usage: m=MaskFromBoundaryTag(domain,"left", "right")     Usage: m=MaskFromBoundaryTag(domain, "left", "right")
1724    
1725     @param domain: a given domain     @param domain: domain to be used
1726     @type domain: L{escript.Domain}     @type domain: L{escript.Domain}
1727     @param tags: boundray tags     @param tags: boundary tags
1728     @type tags: C{str}     @type tags: C{str}
1729     @return: a mask which marks samples that are touching the boundary tagged by any of the given tags.     @return: a mask which marks samples that are touching the boundary tagged
1730                by any of the given tags
1731     @rtype: L{escript.Data} of rank 0     @rtype: L{escript.Data} of rank 0
1732     """     """
1733     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
# Line 1621  def MaskFromBoundaryTag(domain,*tags): Line 1735  def MaskFromBoundaryTag(domain,*tags):
1735     for t in tags: d.setTaggedValue(t,1.)     for t in tags: d.setTaggedValue(t,1.)
1736     pde.setValue(y=d)     pde.setValue(y=d)
1737     return util.whereNonZero(pde.getRightHandSide())     return util.whereNonZero(pde.getRightHandSide())
1738  #============================================================================================================================================  
1739    #==============================================================================
1740  class SaddlePointProblem(object):  class SaddlePointProblem(object):
1741     """     """
1742     This implements a solver for a saddlepoint problem     This implements a solver for a saddle point problem
1743    
1744     M{f(u,p)=0}     M{f(u,p)=0}
1745     M{g(u)=0}     M{g(u)=0}
# Line 1634  class SaddlePointProblem(object): Line 1749  class SaddlePointProblem(object):
1749     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})}
1750     M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}     M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}
1751    
1752     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
1753     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
1754     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,
1755       non-linear conjugate gradient method is applied to solve for p, so Q_g plays
1756     in fact the role of a preconditioner.     in fact the role of a preconditioner.
1757     """     """
1758     def __init__(self,verbose=False,*args):     def __init__(self,verbose=False,*args):
1759         """         """
1760         initializes the problem         Initializes the problem.
1761    
1762         @param verbose: switches on the printing out some information         @param verbose: if True, some information is printed in the process
1763         @type verbose: C{bool}         @type verbose: C{bool}
1764         @note: this method may be overwritten by a particular saddle point problem         @note: this method may be overwritten by a particular saddle point
1765                  problem
1766         """         """
1767         print "SaddlePointProblem should not be used anymore!"         print "SaddlePointProblem should not be used anymore!"
1768         if not isinstance(verbose,bool):         if not isinstance(verbose,bool):
# Line 1656  class SaddlePointProblem(object): Line 1773  class SaddlePointProblem(object):
1773    
1774     def trace(self,text):     def trace(self,text):
1775         """         """
1776         prints text if verbose has been set         Prints C{text} if verbose has been set.
1777    
1778         @param text: a text message         @param text: a text message
1779         @type text: C{str}         @type text: C{str}
# Line 1665  class SaddlePointProblem(object): Line 1782  class SaddlePointProblem(object):
1782    
1783     def solve_f(self,u,p,tol=1.e-8):     def solve_f(self,u,p,tol=1.e-8):
1784         """         """
1785         solves         Solves
1786    
1787         A_f du = f(u,p)         A_f du = f(u,p)
1788    
1789         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
1790           to u.
1791    
1792         @param u: current approximation of u         @param u: current approximation of u
1793         @type u: L{escript.Data}         @type u: L{escript.Data}
# Line 1679  class SaddlePointProblem(object): Line 1797  class SaddlePointProblem(object):
1797         @type tol: C{float}         @type tol: C{float}
1798         @return: increment du         @return: increment du
1799         @rtype: L{escript.Data}         @rtype: L{escript.Data}
1800         @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
1801                  problem
1802         """         """
1803         pass         pass
1804    
1805     def solve_g(self,u,tol=1.e-8):     def solve_g(self,u,tol=1.e-8):
1806         """         """
1807         solves         Solves
1808    
1809         Q_g dp = g(u)         Q_g dp = g(u)
1810    
1811         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
1812           Jacobian of g with respect to p.
1813    
1814         @param u: current approximation of u         @param u: current approximation of u
1815         @type u: L{escript.Data}         @type u: L{escript.Data}
# Line 1697  class SaddlePointProblem(object): Line 1817  class SaddlePointProblem(object):
1817         @type tol: C{float}         @type tol: C{float}
1818         @return: increment dp         @return: increment dp
1819         @rtype: L{escript.Data}         @rtype: L{escript.Data}
1820         @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
1821                  problem
1822         """         """
1823         pass         pass
1824    
1825     def inner(self,p0,p1):     def inner(self,p0,p1):
1826         """         """
1827         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
1828           C{integrate(p0*p1)}.
1829         @return: inner product of p0 and p1         @return: inner product of p0 and p1
1830         @rtype: C{float}         @rtype: C{float}
1831         """         """
# Line 1712  class SaddlePointProblem(object): Line 1834  class SaddlePointProblem(object):
1834     subiter_max=3     subiter_max=3
1835     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):
1836          """          """
1837          runs the solver          Runs the solver.
1838    
1839          @param u0: initial guess for C{u}          @param u0: initial guess for C{u}
1840          @type u0: L{esys.escript.Data}          @type u0: L{esys.escript.Data}
# Line 1720  class SaddlePointProblem(object): Line 1842  class SaddlePointProblem(object):
1842          @type p0: L{esys.escript.Data}          @type p0: L{esys.escript.Data}
1843          @param tolerance: tolerance for relative error in C{u} and C{p}          @param tolerance: tolerance for relative error in C{u} and C{p}
1844          @type tolerance: positive C{float}          @type tolerance: positive C{float}
1845          @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
1846                                from C{tolerance}
1847          @type tolerance_u: positive C{float}          @type tolerance_u: positive C{float}
1848          @param iter_max: maximum number of iteration steps.          @param iter_max: maximum number of iteration steps
1849          @type iter_max: C{int}          @type iter_max: C{int}
1850          @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
1851                                     relaxation factor. If C{accepted_reduction=None} no backtracking is used.                                     C{accepted_reduction} backtracking to adapt
1852                                       the relaxation factor. If
1853                                       C{accepted_reduction=None} no backtracking
1854                                       is used.
1855          @type accepted_reduction: positive C{float} or C{None}          @type accepted_reduction: positive C{float} or C{None}
1856          @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.          @param relaxation: initial relaxation factor. If C{relaxation==None},
1857                               the last relaxation factor is used.
1858          @type relaxation: C{float} or C{None}          @type relaxation: C{float} or C{None}
1859          """          """
1860          tol=1.e-2          tol=1.e-2
# Line 1826  class SaddlePointProblem(object): Line 1953  class SaddlePointProblem(object):
1953              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
1954          self.trace("convergence after %s steps."%self.iter)          self.trace("convergence after %s steps."%self.iter)
1955          return u,p          return u,p
1956    

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

  ViewVC Help
Powered by ViewVC 1.1.26