/[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 2168 by gross, Mon Dec 15 05:09:02 2008 UTC revision 2169 by caltinay, Wed Dec 17 03:08:58 2008 UTC
# Line 20  http://www.opensource.org/licenses/osl-3 Line 20  http://www.opensource.org/licenses/osl-3
20  __url__="http://www.uq.edu.au/esscc/escript-finley"  __url__="http://www.uq.edu.au/esscc/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    
435  class CorrectionFailed(SolverSchemeException):  class CorrectionFailed(SolverSchemeException):
436     """     """
437     no convergence has been achieved in the solution correction scheme.     Exception thrown if no convergence has been achieved in the solution
438       correction scheme.
439     """     """
440     pass     pass
441    
442  class IterationBreakDown(SolverSchemeException):  class IterationBreakDown(SolverSchemeException):
443     """     """
444     iteration scheme econouters an incurable breakdown.     Exception thrown if the iteration scheme encountered an incurable breakdown.
445     """     """
446     pass     pass
447    
448  class NegativeNorm(SolverSchemeException):  class NegativeNorm(SolverSchemeException):
449     """     """
450     a norm calculation returns a negative norm.     Exception thrown if a norm calculation returns a negative norm.
451     """     """
452     pass     pass
453    
454  def PCG(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100, initial_guess=True, verbose=False):  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 iteration is terminated if
465    
466     M{|r| <= atol+rtol*|r0|}     M{|r| <= atol+rtol*|r0|}
467    
468     where M{r0} is the initial residual and M{|.|} is the energy norm. In fact     where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
469    
470     M{|r| =sqrt( bilinearform(Msolve(r),r))}     M{|r| = sqrt( bilinearform(Msolve(r),r))}
471    
472     For details on the preconditioned conjugate gradient method see the book:     For details on the preconditioned conjugate gradient method see the book:
473    
474     Templates for the Solution of Linear Systems by R. Barrett, M. Berry,     I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
475     T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,     T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
476     C. Romine, and H. van der Vorst.     C. Romine, and H. van der Vorst}.
477    
478     @param r: initial residual M{r=b-Ax}. C{r} is altered.     @param r: initial residual M{r=b-Ax}. C{r} is altered.
479     @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)     @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
480     @param x: an initial guess for the solution.     @param x: an initial guess for the solution
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 Aprod: returns the value Ax     @param Aprod: returns the value Ax
483     @type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}.     @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
484         The returned object needs to be of the same type like argument C{r}.                  argument C{x}. The returned object needs to be of the same type
485                    like argument C{r}.
486     @param Msolve: solves Mx=r     @param Msolve: solves Mx=r
487     @type Msolve: function C{Msolve(r)} where C{r} is of the same type like argument C{r}.     @type Msolve: function C{Msolve(r)} where C{r} is of the same type like
488        The returned object needs to be of the same type like argument C{x}.                   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>}     @param bilinearform: inner product C{<x,r>}
491     @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is.     @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
492         The returned value is a C{float}.                         type like argument C{x} and C{r} is. The returned value
493                           is a C{float}.
494     @param atol: absolute tolerance     @param atol: absolute tolerance
495     @type atol: non-negative C{float}     @type atol: non-negative C{float}
496     @param rtol: relative tolerance     @param rtol: relative tolerance
497     @type rtol: non-negative C{float}     @type rtol: non-negative C{float}
498     @param iter_max: maximum number of iteration steps.     @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{r} and C{x} are altered.     @warning: C{r} and C{x} are altered.
503     """     """
504     iter=0     iter=0
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)     norm_r0=math.sqrt(rhat_dot_r)
# Line 500  def PCG(r, Aprod, x, Msolve, bilinearfor Line 513  def PCG(r, Aprod, x, Msolve, bilinearfor
513     atol2=max(atol2, 100. * util.EPSILON * norm_r0)     atol2=max(atol2, 100. * util.EPSILON * norm_r0)
514    
515     if verbose: print "PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_r0, atol2)     if verbose: print "PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_r0, atol2)
516      
517    
518     while not math.sqrt(rhat_dot_r) <= atol2:     while not math.sqrt(rhat_dot_r) <= atol2:
519         iter+=1         iter+=1
# Line 525  def PCG(r, Aprod, x, Msolve, bilinearfor Line 538  def PCG(r, Aprod, x, Msolve, bilinearfor
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 571  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 582  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 590  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 614  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 660  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 675  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 697  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:
# Line 719  def __FDGMRES(F0, defect, x0, atol, iter Line 740  def __FDGMRES(F0, defect, x0, atol, iter
740    
741     rho=defect.norm(F0)     rho=defect.norm(F0)
742     if rho<=0.: return x0*0     if rho<=0.: return x0*0
743      
744     v.append(-F0/rho)     v.append(-F0/rho)
745     g[0]=rho     g[0]=rho
746     iter=0     iter=0
747     while rho > atol and iter<iter_restart-1:     while rho > atol and iter<iter_restart-1:
748            if iter  >= iter_max:
749      if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max              raise MaxIterReached,"maximum number of %s steps reached."%iter_max
750    
751          p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)          p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
752      v.append(p)          v.append(p)
753    
754      v_norm1=defect.norm(v[iter+1])          v_norm1=defect.norm(v[iter+1])
755    
756          # Modified Gram-Schmidt          # Modified Gram-Schmidt
757      for j in range(iter+1):          for j in range(iter+1):
758           h[j,iter]=defect.bilinearform(v[j],v[iter+1])                h[j,iter]=defect.bilinearform(v[j],v[iter+1])
759           v[iter+1]-=h[j,iter]*v[j]              v[iter+1]-=h[j,iter]*v[j]
760          
761      h[iter+1,iter]=defect.norm(v[iter+1])          h[iter+1,iter]=defect.norm(v[iter+1])
762      v_norm2=h[iter+1,iter]          v_norm2=h[iter+1,iter]
763    
764          # Reorthogonalize if needed          # Reorthogonalize if needed
765      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)
766          for j in range(iter+1):                for j in range(iter+1):
767             hr=defect.bilinearform(v[j],v[iter+1])                  hr=defect.bilinearform(v[j],v[iter+1])
768                 h[j,iter]=h[j,iter]+hr                  h[j,iter]=h[j,iter]+hr
769                 v[iter+1] -= hr*v[j]                  v[iter+1] -= hr*v[j]
770    
771          v_norm2=defect.norm(v[iter+1])              v_norm2=defect.norm(v[iter+1])
772          h[iter+1,iter]=v_norm2              h[iter+1,iter]=v_norm2
773          #   watch out for happy breakdown          #   watch out for happy breakdown
774          if not v_norm2 == 0:          if not v_norm2 == 0:
775                  v[iter+1]=v[iter+1]/h[iter+1,iter]              v[iter+1]=v[iter+1]/h[iter+1,iter]
776    
777          #   Form and store the information for the new Givens rotation          #   Form and store the information for the new Givens rotation
778      if iter > 0 :          if iter > 0 :
779          hhat=numarray.zeros(iter+1,numarray.Float64)              hhat=numarray.zeros(iter+1,numarray.Float64)
780          for i in range(iter+1) : hhat[i]=h[i,iter]              for i in range(iter+1) : hhat[i]=h[i,iter]
781          hhat=__givapp(c[0:iter],s[0:iter],hhat);              hhat=__givapp(c[0:iter],s[0:iter],hhat);
782              for i in range(iter+1) : h[i,iter]=hhat[i]              for i in range(iter+1) : h[i,iter]=hhat[i]
783    
784      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])
785    
786      if mu!=0 :          if mu!=0 :
787          c[iter]=h[iter,iter]/mu              c[iter]=h[iter,iter]/mu
788          s[iter]=-h[iter+1,iter]/mu              s[iter]=-h[iter+1,iter]/mu
789          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]
790          h[iter+1,iter]=0.0              h[iter+1,iter]=0.0
791          g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])              g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])
792    
793          # Update the residual norm          # Update the residual norm
794          rho=abs(g[iter+1])          rho=abs(g[iter+1])
795      iter+=1          iter+=1
796    
797     # At this point either iter > iter_max or rho < tol.     # At this point either iter > iter_max or rho < tol.
798     # It's time to compute x and leave.             # It's time to compute x and leave.
799     if iter > 0 :     if iter > 0 :
800       y=numarray.zeros(iter,numarray.Float64)           y=numarray.zeros(iter,numarray.Float64)
801       y[iter-1] = g[iter-1] / h[iter-1,iter-1]       y[iter-1] = g[iter-1] / h[iter-1,iter-1]
802       if iter > 1 :         if iter > 1 :
803          i=iter-2            i=iter-2
804          while i>=0 :          while i>=0 :
805            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]
806            i=i-1            i=i-1
807       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
808       for i in range(iter-1):       for i in range(iter-1):
809      xhat += v[i]*y[i]      xhat += v[i]*y[i]
810     else :     else :
811        xhat=v[0] * 0        xhat=v[0] * 0
812    
813     if iter<iter_restart-1:     if iter<iter_restart-1:
814        stopped=iter        stopped=iter
815     else:     else:
816        stopped=-1        stopped=-1
817    
818     return xhat,stopped     return xhat,stopped
819    
820  def GMRES(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100, iter_restart=20, verbose=False):  def GMRES(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100, iter_restart=20, verbose=False):
821     """     """
822     Solver for     Solver for
823    
824     M{Ax=b}     M{Ax=b}
825    
826     with a general operator A (more details required!).     with a general operator A (more details required!).
827     It uses the generalized minimum residual method (GMRES).     It uses the generalized minimum residual method (GMRES).
828    
829     The iteration is terminated if       The iteration is terminated if
830    
831     M{|r| <= atol+rtol*|r0|}     M{|r| <= atol+rtol*|r0|}
832    
833     where M{r0} is the initial residual and M{|.|} is the energy norm. In fact     where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
834    
835     M{|r| =sqrt( bilinearform(r,r))}     M{|r| = sqrt( bilinearform(r,r))}
836    
837     @param r: initial residual M{r=b-Ax}. C{r} is altered.     @param r: initial residual M{r=b-Ax}. C{r} is altered.
838     @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)     @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
839     @param x: an initial guess for the solution.     @param x: an initial guess for the solution
840     @type x: same like C{r}     @type x: same like C{r}
841     @param Aprod: returns the value Ax     @param Aprod: returns the value Ax
842     @type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}.     @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
843         The returned object needs to be of the same type like argument C{r}.                  argument C{x}. The returned object needs to be of the same
844                    type like argument C{r}.
845     @param bilinearform: inner product C{<x,r>}     @param bilinearform: inner product C{<x,r>}
846     @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is.     @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
847         The returned value is a C{float}.                         type like argument C{x} and C{r}. The returned value is
848                           a C{float}.
849     @param atol: absolute tolerance     @param atol: absolute tolerance
850     @type atol: non-negative C{float}     @type atol: non-negative C{float}
851     @param rtol: relative tolerance     @param rtol: relative tolerance
852     @type rtol: non-negative C{float}     @type rtol: non-negative C{float}
853     @param iter_max: maximum number of iteration steps.     @param iter_max: maximum number of iteration steps
854     @type iter_max: C{int}     @type iter_max: C{int}
855     @param iter_restart: in order to ssave memory the orthogonalization process is terminated after C{iter_restart} steps and the     @param iter_restart: in order to save memory the orthogonalization process
856                      iteration is restarted.                          is terminated after C{iter_restart} steps and the
857                            iteration is restarted.
858     @type iter_restart: C{int}     @type iter_restart: C{int}
859     @return: the solution approximation and the corresponding residual.     @return: the solution approximation and the corresponding residual
860     @rtype: C{tuple}     @rtype: C{tuple}
861     @warning: C{r} and C{x} are altered.     @warning: C{r} and C{x} are altered.
862     """     """
863     m=iter_restart     m=iter_restart
# Line 849  def GMRES(r, Aprod, x, bilinearform, ato Line 873  def GMRES(r, Aprod, x, bilinearform, ato
873        if verbose: print "GMRES: absolute tolerance = %e"%atol2        if verbose: print "GMRES: absolute tolerance = %e"%atol2
874     if atol2<=0:     if atol2<=0:
875        raise ValueError,"Non-positive tolarance."        raise ValueError,"Non-positive tolarance."
876      
877     while True:     while True:
878        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
879        if restarted:        if restarted:
880           r2 = r-Aprod(x-x2)           r2 = r-Aprod(x-x2)
881        else:        else:
882           r2=1*r           r2=1*r
883        x2=x*1.        x2=x*1.
884        x,stopped=__GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose)        x,stopped=__GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose)
885        iter+=iter_restart            iter+=iter_restart
886        if stopped: break        if stopped: break
887        if verbose: print "GMRES: restart."        if verbose: print "GMRES: restart."
888        restarted=True        restarted=True
# Line 867  def GMRES(r, Aprod, x, bilinearform, ato Line 891  def GMRES(r, Aprod, x, bilinearform, ato
891    
892  def __GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False):  def __GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False):
893     iter=0     iter=0
894      
895     h=numarray.zeros((iter_restart+1,iter_restart),numarray.Float64)     h=numarray.zeros((iter_restart+1,iter_restart),numarray.Float64)
896     c=numarray.zeros(iter_restart,numarray.Float64)     c=numarray.zeros(iter_restart,numarray.Float64)
897     s=numarray.zeros(iter_restart,numarray.Float64)     s=numarray.zeros(iter_restart,numarray.Float64)
# Line 877  def __GMRESm(r, Aprod, x, bilinearform, Line 901  def __GMRESm(r, Aprod, x, bilinearform,
901     r_dot_r = bilinearform(r, r)     r_dot_r = bilinearform(r, r)
902     if r_dot_r<0: raise NegativeNorm,"negative norm."     if r_dot_r<0: raise NegativeNorm,"negative norm."
903     rho=math.sqrt(r_dot_r)     rho=math.sqrt(r_dot_r)
904      
905     v.append(r/rho)     v.append(r/rho)
906     g[0]=rho     g[0]=rho
907    
# Line 889  def __GMRESm(r, Aprod, x, bilinearform, Line 913  def __GMRESm(r, Aprod, x, bilinearform,
913      p=Aprod(v[iter])      p=Aprod(v[iter])
914      v.append(p)      v.append(p)
915    
916      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
917    
918  # Modified Gram-Schmidt  # Modified Gram-Schmidt
919      for j in range(iter+1):      for j in range(iter+1):
920        h[j,iter]=bilinearform(v[j],v[iter+1])          h[j,iter]=bilinearform(v[j],v[iter+1])
921        v[iter+1]-=h[j,iter]*v[j]        v[iter+1]-=h[j,iter]*v[j]
922          
923      h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))      h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
924      v_norm2=h[iter+1,iter]      v_norm2=h[iter+1,iter]
925    
926  # Reorthogonalize if needed  # Reorthogonalize if needed
927      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)
928       for j in range(iter+1):         for j in range(iter+1):
929          hr=bilinearform(v[j],v[iter+1])          hr=bilinearform(v[j],v[iter+1])
930              h[j,iter]=h[j,iter]+hr              h[j,iter]=h[j,iter]+hr
931              v[iter+1] -= hr*v[j]              v[iter+1] -= hr*v[j]
932    
933       v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
934       h[iter+1,iter]=v_norm2       h[iter+1,iter]=v_norm2
935    
936  #   watch out for happy breakdown  #   watch out for happy breakdown
937          if not v_norm2 == 0:          if not v_norm2 == 0:
938           v[iter+1]=v[iter+1]/h[iter+1,iter]           v[iter+1]=v[iter+1]/h[iter+1,iter]
939    
# Line 924  def __GMRESm(r, Aprod, x, bilinearform, Line 948  def __GMRESm(r, Aprod, x, bilinearform,
948          h[iter+1,iter]=0.0          h[iter+1,iter]=0.0
949          g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])          g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])
950  # Update the residual norm  # Update the residual norm
951                  
952          rho=abs(g[iter+1])          rho=abs(g[iter+1])
953          if verbose: print "GMRES: iteration step %s: residual %e"%(iter,rho)          if verbose: print "GMRES: iteration step %s: residual %e"%(iter,rho)
954      iter+=1      iter+=1
955    
956  # At this point either iter > iter_max or rho < tol.  # At this point either iter > iter_max or rho < tol.
957  # It's time to compute x and leave.          # It's time to compute x and leave.
958    
959     if verbose: print "GMRES: iteration stopped after %s step."%iter     if verbose: print "GMRES: iteration stopped after %s step."%iter
960     if iter > 0 :     if iter > 0 :
961       y=numarray.zeros(iter,numarray.Float64)           y=numarray.zeros(iter,numarray.Float64)
962       y[iter-1] = g[iter-1] / h[iter-1,iter-1]       y[iter-1] = g[iter-1] / h[iter-1,iter-1]
963       if iter > 1 :         if iter > 1 :
964          i=iter-2            i=iter-2
965          while i>=0 :          while i>=0 :
966            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]
967            i=i-1            i=i-1
968       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
969       for i in range(iter-1):       for i in range(iter-1):
970      xhat += v[i]*y[i]      xhat += v[i]*y[i]
971     else:     else:
972       xhat=v[0] * 0       xhat=v[0] * 0
973    
974     x += xhat     x += xhat
975     if iter<iter_restart-1:     if iter<iter_restart-1:
976        stopped=True        stopped=True
977     else:     else:
978        stopped=False        stopped=False
979    
980     return x,stopped     return x,stopped
981    
982  def MINRES(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100):  def MINRES(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
983      """      """
984      Solver for      Solver for
985    
986      M{Ax=b}      M{Ax=b}
987    
988      with a symmetric and positive definite operator A (more details required!).      with a symmetric and positive definite operator A (more details required!).
989      It uses the minimum residual method (MINRES) with preconditioner M providing an approximation of A.      It uses the minimum residual method (MINRES) with preconditioner M
990        providing an approximation of A.
991      The iteration is terminated if    
992        The iteration is terminated if
993    
994      M{|r| <= atol+rtol*|r0|}      M{|r| <= atol+rtol*|r0|}
995    
996      where M{r0} is the initial residual and M{|.|} is the energy norm. In fact      where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
997    
998      M{|r| =sqrt( bilinearform(Msolve(r),r))}      M{|r| = sqrt( bilinearform(Msolve(r),r))}
999    
1000      For details on the preconditioned conjugate gradient method see the book:      For details on the preconditioned conjugate gradient method see the book:
1001    
1002      Templates for the Solution of Linear Systems by R. Barrett, M. Berry,      I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
1003      T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,      T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
1004      C. Romine, and H. van der Vorst.      C. Romine, and H. van der Vorst}.
1005    
1006      @param r: initial residual M{r=b-Ax}. C{r} is altered.      @param r: initial residual M{r=b-Ax}. C{r} is altered.
1007      @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)      @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1008      @param x: an initial guess for the solution.      @param x: an initial guess for the solution
1009      @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)
1010      @param Aprod: returns the value Ax      @param Aprod: returns the value Ax
1011      @type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}.      @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
1012          The returned object needs to be of the same type like argument C{r}.                   argument C{x}. The returned object needs to be of the same
1013                     type like argument C{r}.
1014      @param Msolve: solves Mx=r      @param Msolve: solves Mx=r
1015      @type Msolve: function C{Msolve(r)} where C{r} is of the same type like argument C{r}.      @type Msolve: function C{Msolve(r)} where C{r} is of the same type like
1016          The returned object needs to be of the same type like argument C{x}.                    argument C{r}. The returned object needs to be of the same
1017                      type like argument C{x}.
1018      @param bilinearform: inner product C{<x,r>}      @param bilinearform: inner product C{<x,r>}
1019      @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is.      @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
1020         The returned value is a C{float}.                          type like argument C{x} and C{r} is. The returned value
1021                            is a C{float}.
1022      @param atol: absolute tolerance      @param atol: absolute tolerance
1023      @type atol: non-negative C{float}      @type atol: non-negative C{float}
1024      @param rtol: relative tolerance      @param rtol: relative tolerance
1025      @type rtol: non-negative C{float}      @type rtol: non-negative C{float}
1026      @param iter_max: maximum number of iteration steps.      @param iter_max: maximum number of iteration steps
1027      @type iter_max: C{int}      @type iter_max: C{int}
1028      @return: the solution approximation and the corresponding residual      @return: the solution approximation and the corresponding residual
1029      @rtype: C{tuple}      @rtype: C{tuple}
1030      @warning: C{r} and C{x} are altered.      @warning: C{r} and C{x} are altered.
1031      """      """
1032      #------------------------------------------------------------------      #------------------------------------------------------------------
# Line 1009  def MINRES(r, Aprod, x, Msolve, bilinear Line 1037  def MINRES(r, Aprod, x, Msolve, bilinear
1037      r1    = r      r1    = r
1038      y = Msolve(r)      y = Msolve(r)
1039      beta1 = bilinearform(y,r)      beta1 = bilinearform(y,r)
1040    
1041      if beta1< 0: raise NegativeNorm,"negative norm."      if beta1< 0: raise NegativeNorm,"negative norm."
1042    
1043      #  If r = 0 exactly, stop with x      #  If r = 0 exactly, stop with x
1044      if beta1==0: return x      if beta1==0: return x
1045    
1046      if beta1> 0: beta1  = math.sqrt(beta1)            if beta1> 0: beta1  = math.sqrt(beta1)
1047    
1048      #------------------------------------------------------------------      #------------------------------------------------------------------
1049      # Initialize quantities.      # Initialize quantities.
# Line 1062  def MINRES(r, Aprod, x, Msolve, bilinear Line 1090  def MINRES(r, Aprod, x, Msolve, bilinear
1090          #-----------------------------------------------------------------          #-----------------------------------------------------------------
1091          s = 1/beta                 # Normalize previous vector (in y).          s = 1/beta                 # Normalize previous vector (in y).
1092          v = s*y                    # v = vk if P = I          v = s*y                    # v = vk if P = I
1093        
1094          y      = Aprod(v)          y      = Aprod(v)
1095        
1096          if iter >= 2:          if iter >= 2:
1097            y = y - (beta/oldb)*r1            y = y - (beta/oldb)*r1
1098    
1099          alfa   = bilinearform(v,y)              # alphak          alfa   = bilinearform(v,y)              # alphak
1100          y      += (- alfa/beta)*r2          y      += (- alfa/beta)*r2
1101          r1     = r2          r1     = r2
1102          r2     = y          r2     = y
1103          y = Msolve(r2)          y = Msolve(r2)
# Line 1079  def MINRES(r, Aprod, x, Msolve, bilinear Line 1107  def MINRES(r, Aprod, x, Msolve, bilinear
1107    
1108          beta   = math.sqrt( beta )          beta   = math.sqrt( beta )
1109          tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta          tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
1110            
1111          if iter==1:                 # Initialize a few things.          if iter==1:                 # Initialize a few things.
1112            gmax   = abs( alfa )      # alpha1            gmax   = abs( alfa )      # alpha1
1113            gmin   = gmax             # alpha1            gmin   = gmax             # alpha1
# Line 1087  def MINRES(r, Aprod, x, Msolve, bilinear Line 1115  def MINRES(r, Aprod, x, Msolve, bilinear
1115          # Apply previous rotation Qk-1 to get          # Apply previous rotation Qk-1 to get
1116          #   [deltak epslnk+1] = [cs  sn][dbark    0   ]          #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
1117          #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].          #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
1118        
1119          oldeps = epsln          oldeps = epsln
1120          delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak          delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
1121          gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k          gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
# Line 1097  def MINRES(r, Aprod, x, Msolve, bilinear Line 1125  def MINRES(r, Aprod, x, Msolve, bilinear
1125          # Compute the next plane rotation Qk          # Compute the next plane rotation Qk
1126    
1127          gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak          gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
1128          gamma  = max(gamma,eps)          gamma  = max(gamma,eps)
1129          cs     = gbar / gamma             # ck          cs     = gbar / gamma             # ck
1130          sn     = beta / gamma             # sk          sn     = beta / gamma             # sk
1131          phi    = cs * phibar              # phik          phi    = cs * phibar              # phik
# Line 1105  def MINRES(r, Aprod, x, Msolve, bilinear Line 1133  def MINRES(r, Aprod, x, Msolve, bilinear
1133    
1134          # Update  x.          # Update  x.
1135    
1136          denom = 1/gamma          denom = 1/gamma
1137          w1    = w2          w1    = w2
1138          w2    = w          w2    = w
1139          w     = (v - oldeps*w1 - delta*w2) * denom          w     = (v - oldeps*w1 - delta*w2) * denom
1140          x     +=  phi*w          x     +=  phi*w
1141    
# Line 1122  def MINRES(r, Aprod, x, Msolve, bilinear Line 1150  def MINRES(r, Aprod, x, Msolve, bilinear
1150    
1151          # Estimate various norms and test for convergence.          # Estimate various norms and test for convergence.
1152    
1153          Anorm  = math.sqrt( tnorm2 )          Anorm  = math.sqrt( tnorm2 )
1154          ynorm  = math.sqrt( ynorm2 )          ynorm  = math.sqrt( ynorm2 )
1155    
1156          rnorm  = phibar          rnorm  = phibar
1157    
# Line 1131  def MINRES(r, Aprod, x, Msolve, bilinear Line 1159  def MINRES(r, Aprod, x, Msolve, bilinear
1159    
1160  def TFQMR(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100):  def TFQMR(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1161    """    """
1162    Solver for    Solver for
1163    
1164    M{Ax=b}    M{Ax=b}
1165    
1166    with a general operator A (more details required!).    with a general operator A (more details required!).
1167    It uses the generalized minimum residual method (GMRES).    It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).
1168    
1169    The iteration is terminated if      The iteration is terminated if
1170    
1171    M{|r| <= atol+rtol*|r0|}    M{|r| <= atol+rtol*|r0|}
1172    
1173    where M{r0} is the initial residual and M{|.|} is the energy norm. In fact    where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
1174    
1175    M{|r| =sqrt( bilinearform(r,r))}    M{|r| = sqrt( bilinearform(r,r))}
1176    
1177    @param r: initial residual M{r=b-Ax}. C{r} is altered.    @param r: initial residual M{r=b-Ax}. C{r} is altered.
1178    @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)    @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1179    @param x: an initial guess for the solution.    @param x: an initial guess for the solution
1180    @type x: same like C{r}    @type x: same like C{r}
1181    @param Aprod: returns the value Ax    @param Aprod: returns the value Ax
1182    @type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}.    @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
1183        The returned object needs to be of the same type like argument C{r}.                 argument C{x}. The returned object needs to be of the same type
1184                   like argument C{r}.
1185    @param bilinearform: inner product C{<x,r>}    @param bilinearform: inner product C{<x,r>}
1186    @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is.    @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
1187        The returned value is a C{float}.                        type like argument C{x} and C{r}. The returned value is
1188                          a C{float}.
1189    @param atol: absolute tolerance    @param atol: absolute tolerance
1190    @type atol: non-negative C{float}    @type atol: non-negative C{float}
1191    @param rtol: relative tolerance    @param rtol: relative tolerance
1192    @type rtol: non-negative C{float}    @type rtol: non-negative C{float}
1193    @param iter_max: maximum number of iteration steps.    @param iter_max: maximum number of iteration steps
1194    @type iter_max: C{int}    @type iter_max: C{int}
1195    @rtype: C{tuple}    @rtype: C{tuple}
1196    @warning: C{r} and C{x} are altered.    @warning: C{r} and C{x} are altered.
1197    """    """
1198    u1=0    u1=0
# Line 1171  def TFQMR(r, Aprod, x, bilinearform, ato Line 1201  def TFQMR(r, Aprod, x, bilinearform, ato
1201    y2=0    y2=0
1202    
1203    w = r    w = r
1204    y1 = r    y1 = r
1205    iter = 0    iter = 0
1206    d = 0    d = 0
1207    v = Aprod(y1)    v = Aprod(y1)
1208    u1 = v    u1 = v
1209      
1210    theta = 0.0;    theta = 0.0;
1211    eta = 0.0;    eta = 0.0;
1212    rho=bilinearform(r,r)    rho=bilinearform(r,r)
# Line 1188  def TFQMR(r, Aprod, x, bilinearform, ato Line 1218  def TFQMR(r, Aprod, x, bilinearform, ato
1218    
1219      sigma = bilinearform(r,v)      sigma = bilinearform(r,v)
1220      if sigma == 0.0: raise IterationBreakDown,'TFQMR breakdown, sigma=0'      if sigma == 0.0: raise IterationBreakDown,'TFQMR breakdown, sigma=0'
1221    
1222      alpha = rho / sigma      alpha = rho / sigma
1223    
1224      for j in range(2):      for j in range(2):
# Line 1197  def TFQMR(r, Aprod, x, bilinearform, ato Line 1228  def TFQMR(r, Aprod, x, bilinearform, ato
1228        if ( j == 1 ):        if ( j == 1 ):
1229          y2 = y1 - alpha * v          y2 = y1 - alpha * v
1230          u2 = Aprod(y2)          u2 = Aprod(y2)
1231    
1232        m = 2 * (iter+1) - 2 + (j+1)        m = 2 * (iter+1) - 2 + (j+1)
1233        if j==0:        if j==0:
1234           w = w - alpha * u1           w = w - alpha * u1
1235           d = y1 + ( theta * theta * eta / alpha ) * d           d = y1 + ( theta * theta * eta / alpha ) * d
1236        if j==1:        if j==1:
# Line 1222  def TFQMR(r, Aprod, x, bilinearform, ato Line 1253  def TFQMR(r, Aprod, x, bilinearform, ato
1253      y1 = w + beta * y2;      y1 = w + beta * y2;
1254      u1 = Aprod(y1)      u1 = Aprod(y1)
1255      v = u1 + beta * ( u2 + beta * v )      v = u1 + beta * ( u2 + beta * v )
1256        
1257      iter += 1      iter += 1
1258    
1259    return x    return x
# Line 1232  def TFQMR(r, Aprod, x, bilinearform, ato Line 1263  def TFQMR(r, Aprod, x, bilinearform, ato
1263    
1264  class ArithmeticTuple(object):  class ArithmeticTuple(object):
1265     """     """
1266     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
1267       ArithmeticTuple and C{a} is a float.
1268    
1269     example of usage:     Example of usage::
1270    
1271     from esys.escript import Data         from esys.escript import Data
1272     from numarray import array         from numarray import array
1273     a=Data(...)         a=Data(...)
1274     b=array([1.,4.])         b=array([1.,4.])
1275     x=ArithmeticTuple(a,b)         x=ArithmeticTuple(a,b)
1276     y=5.*x         y=5.*x
1277    
1278     """     """
1279     def __init__(self,*args):     def __init__(self,*args):
1280         """         """
1281         initialize object with elements args.         Initializes object with elements C{args}.
1282    
1283         @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
1284                        scaling (x=a*y)
1285         """         """
1286         self.__items=list(args)         self.__items=list(args)
1287    
1288     def __len__(self):     def __len__(self):
1289         """         """
1290         number of items         Returns the number of items.
1291    
1292         @return: number of items         @return: number of items
1293         @rtype: C{int}         @rtype: C{int}
# Line 1263  class ArithmeticTuple(object): Line 1296  class ArithmeticTuple(object):
1296    
1297     def __getitem__(self,index):     def __getitem__(self,index):
1298         """         """
1299         get an item         Returns item at specified position.
1300    
1301         @param index: item to be returned         @param index: index of item to be returned
1302         @type index: C{int}         @type index: C{int}
1303         @return: item with index C{index}         @return: item with index C{index}
1304         """         """
# Line 1273  class ArithmeticTuple(object): Line 1306  class ArithmeticTuple(object):
1306    
1307     def __mul__(self,other):     def __mul__(self,other):
1308         """         """
1309         scaling from the right         Scales by C{other} from the right.
1310    
1311         @param other: scaling factor         @param other: scaling factor
1312         @type other: C{float}         @type other: C{float}
# Line 1281  class ArithmeticTuple(object): Line 1314  class ArithmeticTuple(object):
1314         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1315         """         """
1316         out=[]         out=[]
1317         try:           try:
1318             l=len(other)             l=len(other)
1319             if l!=len(self):             if l!=len(self):
1320                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1321             for i in range(l): out.append(self[i]*other[i])             for i in range(l): out.append(self[i]*other[i])
1322         except TypeError:         except TypeError:
1323         for i in range(len(self)): out.append(self[i]*other)             for i in range(len(self)): out.append(self[i]*other)
1324         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1325    
1326     def __rmul__(self,other):     def __rmul__(self,other):
1327         """         """
1328         scaling from the left         Scales by C{other} from the left.
1329    
1330         @param other: scaling factor         @param other: scaling factor
1331         @type other: C{float}         @type other: C{float}
# Line 1300  class ArithmeticTuple(object): Line 1333  class ArithmeticTuple(object):
1333         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1334         """         """
1335         out=[]         out=[]
1336         try:           try:
1337             l=len(other)             l=len(other)
1338             if l!=len(self):             if l!=len(self):
1339                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1340             for i in range(l): out.append(other[i]*self[i])             for i in range(l): out.append(other[i]*self[i])
1341         except TypeError:         except TypeError:
1342         for i in range(len(self)): out.append(other*self[i])             for i in range(len(self)): out.append(other*self[i])
1343         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1344    
1345     def __div__(self,other):     def __div__(self,other):
1346         """         """
1347         dividing from the right         Scales by (1/C{other}) from the right.
1348    
1349         @param other: scaling factor         @param other: scaling factor
1350         @type other: C{float}         @type other: C{float}
# Line 1322  class ArithmeticTuple(object): Line 1355  class ArithmeticTuple(object):
1355    
1356     def __rdiv__(self,other):     def __rdiv__(self,other):
1357         """         """
1358         dividing from the left         Scales by (1/C{other}) from the left.
1359    
1360         @param other: scaling factor         @param other: scaling factor
1361         @type other: C{float}         @type other: C{float}
# Line 1330  class ArithmeticTuple(object): Line 1363  class ArithmeticTuple(object):
1363         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1364         """         """
1365         out=[]         out=[]
1366         try:           try:
1367             l=len(other)             l=len(other)
1368             if l!=len(self):             if l!=len(self):
1369                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1370             for i in range(l): out.append(other[i]/self[i])             for i in range(l): out.append(other[i]/self[i])
1371         except TypeError:         except TypeError:
1372         for i in range(len(self)): out.append(other/self[i])             for i in range(len(self)): out.append(other/self[i])
1373         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1374      
1375     def __iadd__(self,other):     def __iadd__(self,other):
1376         """         """
1377         in-place add of other to self         Inplace addition of C{other} to self.
1378    
1379         @param other: increment         @param other: increment
1380         @type other: C{ArithmeticTuple}         @type other: C{ArithmeticTuple}
1381         """         """
1382         if len(self) != len(other):         if len(self) != len(other):
1383             raise ValueError,"tuple length must match."             raise ValueError,"tuple lengths must match."
1384         for i in range(len(self)):         for i in range(len(self)):
1385             self.__items[i]+=other[i]             self.__items[i]+=other[i]
1386         return self         return self
1387    
1388     def __add__(self,other):     def __add__(self,other):
1389         """         """
1390         add other to self         Adds C{other} to self.
1391    
1392         @param other: increment         @param other: increment
1393         @type other: C{ArithmeticTuple}         @type other: C{ArithmeticTuple}
1394         """         """
1395         out=[]         out=[]
1396         try:           try:
1397             l=len(other)             l=len(other)
1398             if l!=len(self):             if l!=len(self):
1399                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1400             for i in range(l): out.append(self[i]+other[i])             for i in range(l): out.append(self[i]+other[i])
1401         except TypeError:         except TypeError:
1402         for i in range(len(self)): out.append(self[i]+other)             for i in range(len(self)): out.append(self[i]+other)
1403         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1404    
1405     def __sub__(self,other):     def __sub__(self,other):
1406         """         """
1407         subtract other from self         Subtracts C{other} from self.
1408    
1409         @param other: increment         @param other: decrement
1410         @type other: C{ArithmeticTuple}         @type other: C{ArithmeticTuple}
1411         """         """
1412         out=[]         out=[]
1413         try:           try:
1414             l=len(other)             l=len(other)
1415             if l!=len(self):             if l!=len(self):
1416                 raise ValueError,"length of of arguments don't match."                 raise ValueError,"length of arguments don't match."
1417             for i in range(l): out.append(self[i]-other[i])             for i in range(l): out.append(self[i]-other[i])
1418         except TypeError:         except TypeError:
1419         for i in range(len(self)): out.append(self[i]-other)             for i in range(len(self)): out.append(self[i]-other)
1420         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1421      
1422     def __isub__(self,other):     def __isub__(self,other):
1423         """         """
1424         subtract other from self         Inplace subtraction of C{other} from self.
1425    
1426         @param other: increment         @param other: decrement
1427         @type other: C{ArithmeticTuple}         @type other: C{ArithmeticTuple}
1428         """         """
1429         if len(self) != len(other):         if len(self) != len(other):
# Line 1401  class ArithmeticTuple(object): Line 1434  class ArithmeticTuple(object):
1434    
1435     def __neg__(self):     def __neg__(self):
1436         """         """
1437         negate         Negates values.
   
1438         """         """
1439         out=[]         out=[]
1440         for i in range(len(self)):         for i in range(len(self)):
# Line 1412  class ArithmeticTuple(object): Line 1444  class ArithmeticTuple(object):
1444    
1445  class HomogeneousSaddlePointProblem(object):  class HomogeneousSaddlePointProblem(object):
1446        """        """
1447        This provides a framwork for solving linear homogeneous saddle point problem of the form        This class provides a framework for solving linear homogeneous saddle
1448          point problems of the form::
1449    
1450               Av+B^*p=f            M{Av+B^*p=f}
1451               Bv    =0            M{Bv     =0}
1452    
1453        for the unknowns v and p and given operators A and B and given right hand side f.        for the unknowns M{v} and M{p} and given operators M{A} and M{B} and
1454        B^* is the adjoint operator of B.        given right hand side M{f}. M{B^*} is the adjoint operator of M{B}.
1455        """        """
1456        def __init__(self,**kwargs):        def __init__(self,**kwargs):
1457          self.setTolerance()          self.setTolerance()
# Line 1428  class HomogeneousSaddlePointProblem(obje Line 1461  class HomogeneousSaddlePointProblem(obje
1461        #=============================================================        #=============================================================
1462        def initialize(self):        def initialize(self):
1463          """          """
1464          initialize the problem (overwrite)          Initializes the problem (overwrite).
1465          """          """
1466          pass          pass
1467    
1468        def B(self,v):        def B(self,v):
1469          """          """
1470          returns Bv (overwrite)          Returns Bv (overwrite).
         @rtype: equal to the type of p  
1471    
1472            @rtype: equal to the type of p
1473          @note: boundary conditions on p should be zero!          @note: boundary conditions on p should be zero!
1474          """          """
1475          raise NotImplementedError,"no operator B implemented."          raise NotImplementedError, "no operator B implemented."
1476    
1477        def inner_pBv(self,p,Bv):        def inner_pBv(self,p,Bv):
1478           """           """
1479           returns inner product of element p and Bv  (overwrite)           Returns inner product of element p and Bv (overwrite).
1480            
1481           @type p: equal to the type of p           @type p: equal to the type of p
1482           @type Bv: equal to the type of result of operator B           @type Bv: equal to the type of result of operator B
          @rtype: C{float}  
   
1483           @rtype: equal to the type of p           @rtype: equal to the type of p
1484           """           """
1485           raise NotImplementedError,"no inner product for p implemented."           raise NotImplementedError,"no inner product for p implemented."
1486    
1487        def inner_p(self,p0,p1):        def inner_p(self,p0,p1):
1488           """           """
1489           returns inner product of element p0 and p1  (overwrite)           Returns inner product of element p0 and p1 (overwrite).
1490            
1491           @type p0: equal to the type of p           @type p0: equal to the type of p
1492           @type p1: equal to the type of p           @type p1: equal to the type of p
1493           @rtype: equal to the type of p           @rtype: equal to the type of p
# Line 1465  class HomogeneousSaddlePointProblem(obje Line 1496  class HomogeneousSaddlePointProblem(obje
1496    
1497        def inner_v(self,v0,v1):        def inner_v(self,v0,v1):
1498           """           """
1499           returns inner product of two element v0 and v1  (overwrite)           Returns inner product of two elements v0 and v1 (overwrite).
1500            
1501           @type v0: equal to the type of v           @type v0: equal to the type of v
1502           @type v1: equal to the type of v           @type v1: equal to the type of v
          @rtype: C{float}  
   
1503           @rtype: equal to the type of v           @rtype: equal to the type of v
1504           """           """
1505           raise NotImplementedError,"no inner product for v implemented."           raise NotImplementedError,"no inner product for v implemented."
# Line 1478  class HomogeneousSaddlePointProblem(obje Line 1507  class HomogeneousSaddlePointProblem(obje
1507    
1508        def solve_A(self,u,p):        def solve_A(self,u,p):
1509           """           """
1510           solves Av=f-Au-B^*p with accuracy self.getSubProblemTolerance() (overwrite)           Solves M{Av=f-Au-B^*p} with accuracy L{self.getSubProblemTolerance()}
1511             (overwrite).
1512    
1513           @rtype: equal to the type of v           @rtype: equal to the type of v
1514           @note: boundary conditions on v should be zero!           @note: boundary conditions on v should be zero!
# Line 1487  class HomogeneousSaddlePointProblem(obje Line 1517  class HomogeneousSaddlePointProblem(obje
1517    
1518        def solve_prec(self,p):        def solve_prec(self,p):
1519           """           """
1520           provides a preconditioner for BA^{-1}B^* with accuracy self.getSubProblemTolerance() (overwrite)           Provides a preconditioner for M{BA^{-1}B^*} with accuracy
1521             L{self.getSubProblemTolerance()} (overwrite).
1522    
1523           @rtype: equal to the type of p           @rtype: equal to the type of p
1524           @note: boundary conditions on p should be zero!           @note: boundary conditions on p should be zero!
# Line 1495  class HomogeneousSaddlePointProblem(obje Line 1526  class HomogeneousSaddlePointProblem(obje
1526           raise NotImplementedError,"no preconditioner for Schur complement implemented."           raise NotImplementedError,"no preconditioner for Schur complement implemented."
1527        #=============================================================        #=============================================================
1528        def __Aprod_PCG(self,p):        def __Aprod_PCG(self,p):
1529            # return (v,Bv) with v=A^-1B*p            # return (v,Bv) with v=A^-1B*p
1530            #solve Av =B^*p as Av =f-Az-B^*(-p)            #solve Av =B^*p as Av =f-Az-B^*(-p)
1531            v=self.solve_A(self.__z,-p)            v=self.solve_A(self.__z,-p)
1532            return ArithmeticTuple(v, self.B(v))            return ArithmeticTuple(v, self.B(v))
# Line 1514  class HomogeneousSaddlePointProblem(obje Line 1545  class HomogeneousSaddlePointProblem(obje
1545            p=x[1]            p=x[1]
1546            w=self.solve_A(self.__z-v,-p)            w=self.solve_A(self.__z-v,-p)
1547            Bw=self.B(w-v)            Bw=self.B(w-v)
1548        q=self.solve_prec(Bw)            q=self.solve_prec(Bw)
1549            return ArithmeticTuple(w,q)            return ArithmeticTuple(w,q)
1550    
1551        def __inner_GMRES(self,a1,a2):        def __inner_GMRES(self,a1,a2):
# Line 1529  class HomogeneousSaddlePointProblem(obje Line 1560  class HomogeneousSaddlePointProblem(obje
1560    
1561        def solve(self,v,p,max_iter=20, verbose=False, show_details=False, useUzawa=True, iter_restart=20,max_correction_steps=3):        def solve(self,v,p,max_iter=20, verbose=False, show_details=False, useUzawa=True, iter_restart=20,max_correction_steps=3):
1562           """           """
1563           solves the saddle point problem using initial guesses v and p.           Solves the saddle point problem using initial guesses v and p.
1564    
1565           @param v: initial guess for velocity           @param v: initial guess for velocity
1566           @param p: initial guess for pressure           @param p: initial guess for pressure
1567           @type v: L{Data}           @type v: L{Data}
1568           @type p: L{Data}           @type p: L{Data}
1569           @param useUzawa: indicate the usage of the Uzawa scheme. Otherwise the problem is solved in coupled form.           @param useUzawa: indicates the usage of the Uzawa scheme. Otherwise
1570           @param max_iter: maximum number of iteration steps per correction attempt.                            the problem is solved in coupled form.
1571           @param verbose: show information on the progress of the saddlepoint problem solver.           @param max_iter: maximum number of iteration steps per correction
1572           @param show_details: show details of the sub problems solves                            attempt
1573           @param iter_restart: restart the iteration after C{iter_restart} steps (only used if useUzaw=False)           @param verbose: if True, shows information on the progress of the
1574           @param max_correction_steps: maximum number of iteration steps in the attempt get |Bv| to zero.                           saddlepoint problem solver.
1575             @param show_details: if True, shows details of the sub problem solver
1576             @param iter_restart: restart the iteration after C{iter_restart} steps
1577                                  (only used if useUzaw=False)
1578             @param max_correction_steps: maximum number of iteration steps in the
1579                                          attempt to get |Bv| to zero.
1580           @return: new approximations for velocity and pressure           @return: new approximations for velocity and pressure
1581           @type useUzawa: C{bool}           @type useUzawa: C{bool}
1582           @type max_iter: C{int}           @type max_iter: C{int}
# Line 1552  class HomogeneousSaddlePointProblem(obje Line 1588  class HomogeneousSaddlePointProblem(obje
1588           """           """
1589           self.verbose=verbose           self.verbose=verbose
1590           self.show_details=show_details and self.verbose           self.show_details=show_details and self.verbose
1591           #=================================================================================           #=====================================================================
1592           # Az=f is solved as A(z-v)=f-Av-B^*0 (typically with z-v=0 on boundary)           # Az=f is solved as A(z-v)=f-Av-B^*0 (typically with z-v=0 on boundary)
1593           self.__z=v+self.solve_A(v,p*0)           self.__z=v+self.solve_A(v,p*0)
1594           # tolerances:           # tolerances:
# Line 1577  class HomogeneousSaddlePointProblem(obje Line 1613  class HomogeneousSaddlePointProblem(obje
1613           while not converged :           while not converged :
1614              if useUzawa:              if useUzawa:
1615                 # assume p is known: then v=z-A^-1B^*p                 # assume p is known: then v=z-A^-1B^*p
1616                 #                       #
1617                 # which leads to BA^-1B^*p = Bz                 # which leads to BA^-1B^*p = Bz
1618                 #                 #
1619                 # note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv                 # note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
# Line 1586  class HomogeneousSaddlePointProblem(obje Line 1622  class HomogeneousSaddlePointProblem(obje
1622                 # the norm of the right hand side Bv = f0                 # the norm of the right hand side Bv = f0
1623                 #                 #
1624                 #                  and the initial residual                 #                  and the initial residual
1625                 #                           #
1626                 #     r=Bz-BA^-1B^*q = B(z-A^{-1}B^*q)=Bw with A(w-z)=Az-Az-B^*q = f -A*0 - B^{*}q                 #     r=Bz-BA^-1B^*q = B(z-A^{-1}B^*q)=Bw with A(w-z)=Az-Az-B^*q = f -A*0 - B^{*}q
1627                 #                 #
1628                 w=self.solve_A(self.__z,q)+self.__z                 w=self.solve_A(self.__z,q)+self.__z
1629                 if self.verbose: print "enter PCG method (iter_max=%s, atol=%e, subtolerance=%e)"%(max_iter,ATOL, self.getSubProblemTolerance())                 if self.verbose: print "enter PCG method (iter_max=%s, atol=%e, subtolerance=%e)"%(max_iter,ATOL, self.getSubProblemTolerance())
1630                 q,r=PCG(ArithmeticTuple(w,self.B(w)),self.__Aprod_PCG,q,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)                 q,r=PCG(ArithmeticTuple(w,self.B(w)),self.__Aprod_PCG,q,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
1631             u=r[0]                   u=r[0]
1632             Bu=r[1]                 Bu=r[1]
1633              else:              else:
1634                 #                 #
1635                 #  with v=dv+z                 #  with v=dv+z
# Line 1607  class HomogeneousSaddlePointProblem(obje Line 1643  class HomogeneousSaddlePointProblem(obje
1643                 if self.verbose: print "enter GMRES (iter_max=%s, atol=%e, subtolerance=%e)"%(max_iter,ATOL,self.getSubProblemTolerance())                 if self.verbose: print "enter GMRES (iter_max=%s, atol=%e, subtolerance=%e)"%(max_iter,ATOL,self.getSubProblemTolerance())
1644                 x=GMRES(ArithmeticTuple(w,self.solve_prec(self.B(u+w))),self.__Aprod_GMRES, ArithmeticTuple(u,q), \                 x=GMRES(ArithmeticTuple(w,self.solve_prec(self.B(u+w))),self.__Aprod_GMRES, ArithmeticTuple(u,q), \
1645                           self.__inner_GMRES,atol=ATOL, rtol=0.,iter_max=max_iter, iter_restart=iter_restart, verbose=self.verbose)                           self.__inner_GMRES,atol=ATOL, rtol=0.,iter_max=max_iter, iter_restart=iter_restart, verbose=self.verbose)
1646             u=x[0]                 u=x[0]
1647                 q=x[1]                 q=x[1]
1648             Bu=self.B(u)                 Bu=self.B(u)
1649              # now we check if |Bu| ~ 0 or more precise |Bu|_p  <= rtol * |v|_v              # now we check if |Bu| ~ 0 or more precise |Bu|_p  <= rtol * |v|_v
1650              nu=self.inner_v(u,u)              nu=self.inner_v(u,u)
1651              p2=self.solve_prec(Bu)              p2=self.solve_prec(Bu)
# Line 1631  class HomogeneousSaddlePointProblem(obje Line 1667  class HomogeneousSaddlePointProblem(obje
1667           if self.verbose: print "saddle point solver: tolerance reached."           if self.verbose: print "saddle point solver: tolerance reached."
1668       return u,q       return u,q
1669    
1670        #==========================================================================================================================        #========================================================================
1671        def setTolerance(self,tolerance=1.e-4):        def setTolerance(self,tolerance=1.e-4):
1672           """           """
1673           sets the relative tolerance for (v,p)           Sets the relative tolerance for (v,p).
1674    
1675           @param tolerance: tolerance to be used           @param tolerance: tolerance to be used
1676           @type tolerance: non-negative C{float}           @type tolerance: non-negative C{float}
1677           """           """
1678           if tolerance<0:           if tolerance<0:
# Line 1646  class HomogeneousSaddlePointProblem(obje Line 1682  class HomogeneousSaddlePointProblem(obje
1682    
1683        def getTolerance(self):        def getTolerance(self):
1684           """           """
1685           returns the relative tolerance           Returns the relative tolerance.
1686    
1687           @return: relative tolerance           @return: relative tolerance
1688           @rtype: C{float}           @rtype: C{float}
1689           """           """
1690           return self.__rtol           return self.__rtol
1691    
1692        def setAbsoluteTolerance(self,tolerance=0.):        def setAbsoluteTolerance(self,tolerance=0.):
1693           """           """
1694           sets the absolute tolerance           Sets the absolute tolerance.
1695    
1696           @param tolerance: tolerance to be used           @param tolerance: tolerance to be used
1697           @type tolerance: non-negative C{float}           @type tolerance: non-negative C{float}
1698           """           """
1699           if tolerance<0:           if tolerance<0:
1700               raise ValueError,"tolerance must be non-negative."               raise ValueError,"tolerance must be non-negative."
1701           self.__atol=tolerance           self.__atol=tolerance
1702    
1703        def getAbsoluteTolerance(self):        def getAbsoluteTolerance(self):
1704           """           """
1705           returns the absolute tolerance           Returns the absolute tolerance.
1706    
1707           @return: absolute tolerance           @return: absolute tolerance
1708           @rtype: C{float}           @rtype: C{float}
1709           """           """
1710           return self.__atol           return self.__atol
1711    
1712        def setSubProblemTolerance(self,rtol=None):        def setSubProblemTolerance(self,rtol=None):
1713           """           """
1714           sets the relative tolerance to solve the subproblem(s).           Sets the relative tolerance to solve the subproblem(s).
1715    
1716           @param rtol: relative tolerence           @param rtol: relative tolerence
1717           @type rtol: positive C{float}           @type rtol: positive C{float}
1718           """           """
1719           if rtol == None:           if rtol == None:
1720                rtol=max(200.*util.EPSILON,self.getTolerance()**2)                rtol=max(200.*util.EPSILON,self.getTolerance()**2)
1721           if rtol<=0:           if rtol<=0:
1722               raise ValueError,"tolerance must be positive."               raise ValueError,"tolerance must be positive."
1723           self.__sub_tol=rtol           self.__sub_tol=rtol
1724    
1725        def getSubProblemTolerance(self):        def getSubProblemTolerance(self):
1726           """           """
1727           returns the subproblem reduction factor           Returns the subproblem reduction factor.
1728    
1729           @return: subproblem reduction factor           @return: subproblem reduction factor
1730           @rtype: C{float}           @rtype: C{float}
1731           """           """
1732           return self.__sub_tol           return self.__sub_tol
1733    
1734  def MaskFromBoundaryTag(domain,*tags):  def MaskFromBoundaryTag(domain,*tags):
1735     """     """
1736     create a mask on the Solution(domain) function space which one for samples     Creates a mask on the Solution(domain) function space where the value is
1737     that touch the boundary tagged by tags.     one for samples that touch the boundary tagged by tags.
1738    
1739     usage: m=MaskFromBoundaryTag(domain,"left", "right")     Usage: m=MaskFromBoundaryTag(domain, "left", "right")
1740    
1741     @param domain: a given domain     @param domain: domain to be used
1742     @type domain: L{escript.Domain}     @type domain: L{escript.Domain}
1743     @param tags: boundray tags     @param tags: boundary tags
1744     @type tags: C{str}     @type tags: C{str}
1745     @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
1746                by any of the given tags
1747     @rtype: L{escript.Data} of rank 0     @rtype: L{escript.Data} of rank 0
1748     """     """
1749     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
# Line 1711  def MaskFromBoundaryTag(domain,*tags): Line 1751  def MaskFromBoundaryTag(domain,*tags):
1751     for t in tags: d.setTaggedValue(t,1.)     for t in tags: d.setTaggedValue(t,1.)
1752     pde.setValue(y=d)     pde.setValue(y=d)
1753     return util.whereNonZero(pde.getRightHandSide())     return util.whereNonZero(pde.getRightHandSide())
1754  #============================================================================================================================================  
1755    #==============================================================================
1756  class SaddlePointProblem(object):  class SaddlePointProblem(object):
1757     """     """
1758     This implements a solver for a saddlepoint problem     This implements a solver for a saddle point problem
1759    
1760     M{f(u,p)=0}     M{f(u,p)=0}
1761     M{g(u)=0}     M{g(u)=0}
# Line 1724  class SaddlePointProblem(object): Line 1765  class SaddlePointProblem(object):
1765     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})}
1766     M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}     M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}
1767    
1768     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
1769     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
1770     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,
1771       non-linear conjugate gradient method is applied to solve for p, so Q_g plays
1772     in fact the role of a preconditioner.     in fact the role of a preconditioner.
1773     """     """
1774     def __init__(self,verbose=False,*args):     def __init__(self,verbose=False,*args):
1775         """         """
1776         initializes the problem         Initializes the problem.
1777    
1778         @param verbose: switches on the printing out some information         @param verbose: if True, some information is printed in the process
1779         @type verbose: C{bool}         @type verbose: C{bool}
1780         @note: this method may be overwritten by a particular saddle point problem         @note: this method may be overwritten by a particular saddle point
1781                  problem
1782         """         """
1783         print "SaddlePointProblem should not be used anymore!"         print "SaddlePointProblem should not be used anymore!"
1784         if not isinstance(verbose,bool):         if not isinstance(verbose,bool):
# Line 1746  class SaddlePointProblem(object): Line 1789  class SaddlePointProblem(object):
1789    
1790     def trace(self,text):     def trace(self,text):
1791         """         """
1792         prints text if verbose has been set         Prints C{text} if verbose has been set.
1793    
1794         @param text: a text message         @param text: a text message
1795         @type text: C{str}         @type text: C{str}
# Line 1755  class SaddlePointProblem(object): Line 1798  class SaddlePointProblem(object):
1798    
1799     def solve_f(self,u,p,tol=1.e-8):     def solve_f(self,u,p,tol=1.e-8):
1800         """         """
1801         solves         Solves
1802    
1803         A_f du = f(u,p)         A_f du = f(u,p)
1804    
1805         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
1806           to u.
1807    
1808         @param u: current approximation of u         @param u: current approximation of u
1809         @type u: L{escript.Data}         @type u: L{escript.Data}
# Line 1769  class SaddlePointProblem(object): Line 1813  class SaddlePointProblem(object):
1813         @type tol: C{float}         @type tol: C{float}
1814         @return: increment du         @return: increment du
1815         @rtype: L{escript.Data}         @rtype: L{escript.Data}
1816         @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
1817                  problem
1818         """         """
1819         pass         pass
1820    
1821     def solve_g(self,u,tol=1.e-8):     def solve_g(self,u,tol=1.e-8):
1822         """         """
1823         solves         Solves
1824    
1825         Q_g dp = g(u)         Q_g dp = g(u)
1826    
1827         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
1828           Jacobian of g with respect to p.
1829    
1830         @param u: current approximation of u         @param u: current approximation of u
1831         @type u: L{escript.Data}         @type u: L{escript.Data}
# Line 1787  class SaddlePointProblem(object): Line 1833  class SaddlePointProblem(object):
1833         @type tol: C{float}         @type tol: C{float}
1834         @return: increment dp         @return: increment dp
1835         @rtype: L{escript.Data}         @rtype: L{escript.Data}
1836         @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
1837                  problem
1838         """         """
1839         pass         pass
1840    
1841     def inner(self,p0,p1):     def inner(self,p0,p1):
1842         """         """
1843         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
1844           C{integrate(p0*p1)}.
1845         @return: inner product of p0 and p1         @return: inner product of p0 and p1
1846         @rtype: C{float}         @rtype: C{float}
1847         """         """
# Line 1802  class SaddlePointProblem(object): Line 1850  class SaddlePointProblem(object):
1850     subiter_max=3     subiter_max=3
1851     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):
1852          """          """
1853          runs the solver          Runs the solver.
1854    
1855          @param u0: initial guess for C{u}          @param u0: initial guess for C{u}
1856          @type u0: L{esys.escript.Data}          @type u0: L{esys.escript.Data}
# Line 1810  class SaddlePointProblem(object): Line 1858  class SaddlePointProblem(object):
1858          @type p0: L{esys.escript.Data}          @type p0: L{esys.escript.Data}
1859          @param tolerance: tolerance for relative error in C{u} and C{p}          @param tolerance: tolerance for relative error in C{u} and C{p}
1860          @type tolerance: positive C{float}          @type tolerance: positive C{float}
1861          @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
1862                                from C{tolerance}
1863          @type tolerance_u: positive C{float}          @type tolerance_u: positive C{float}
1864          @param iter_max: maximum number of iteration steps.          @param iter_max: maximum number of iteration steps
1865          @type iter_max: C{int}          @type iter_max: C{int}
1866          @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
1867                                     relaxation factor. If C{accepted_reduction=None} no backtracking is used.                                     C{accepted_reduction} backtracking to adapt
1868                                       the relaxation factor. If
1869                                       C{accepted_reduction=None} no backtracking
1870                                       is used.
1871          @type accepted_reduction: positive C{float} or C{None}          @type accepted_reduction: positive C{float} or C{None}
1872          @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.          @param relaxation: initial relaxation factor. If C{relaxation==None},
1873                               the last relaxation factor is used.
1874          @type relaxation: C{float} or C{None}          @type relaxation: C{float} or C{None}
1875          """          """
1876          tol=1.e-2          tol=1.e-2
# Line 1916  class SaddlePointProblem(object): Line 1969  class SaddlePointProblem(object):
1969              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
1970          self.trace("convergence after %s steps."%self.iter)          self.trace("convergence after %s steps."%self.iter)
1971          return u,p          return u,p
1972    

Legend:
Removed from v.2168  
changed lines
  Added in v.2169

  ViewVC Help
Powered by ViewVC 1.1.26