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

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

  ViewVC Help
Powered by ViewVC 1.1.26