/[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 782 by bcumming, Tue Jul 18 00:47:47 2006 UTC revision 2498 by gross, Mon Jun 29 04:25:23 2009 UTC
# Line 1  Line 1 
1  # $Id$  
2    ########################################################
3    #
4    # Copyright (c) 2003-2008 by University of Queensland
5    # Earth Systems Science Computational Center (ESSCC)
6    # http://www.uq.edu.au/esscc
7    #
8    # Primary Business: Queensland, Australia
9    # Licensed under the Open Software License version 3.0
10    # http://www.opensource.org/licenses/osl-3.0.php
11    #
12    ########################################################
13    
14    __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15    Earth Systems Science Computational Center (ESSCC)
16    http://www.uq.edu.au/esscc
17    Primary Business: Queensland, Australia"""
18    __license__="""Licensed under the Open Software License version 3.0
19    http://www.opensource.org/licenses/osl-3.0.php"""
20    __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
30    
31  @var __author__: name of author  @var __author__: name of author
32  @var __copyright__: copyrights  @var __copyright__: copyrights
# Line 17  Currently includes: Line 37  Currently includes:
37  """  """
38    
39  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
 __copyright__="""  Copyright (c) 2006 by ACcESS MNRF  
                     http://www.access.edu.au  
                 Primary Business: Queensland, Australia"""  
 __license__="""Licensed under the Open Software License version 3.0  
              http://www.opensource.org/licenses/osl-3.0.php"""  
 __url__="http://www.iservo.edu.au/esys"  
 __version__="$Revision$"  
 __date__="$Date$"  
40    
41    
42  import escript  import escript
43  import linearPDEs  import linearPDEs
44  import numarray  import numpy
45  import util  import util
46    import math
47    
48    ##### Added by Artak
49    # from Numeric import zeros,Int,Float64
50    ###################################
51    
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 50  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 67  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 76  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 93  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 108  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    def __del__(self):      """
155      return      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())
168      if input_data.getRank()==0:      if input_data.getRank()==0:
169          self.__pde.setValue(Y = input_data)          self.__pde.setValue(Y = input_data)
170          out=self.__pde.getSolution()          out=self.__pde.getSolution()
# Line 170  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 218  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 260  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 280  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
328           @type where: L{escript.FunctionSpace}
329           @param x: location(s) of the Locator
330           @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         self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).mindp()         iterative=False
337           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=[]
347               for p in x:
348                  self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())
349           else:
350               self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint()
351    
352       def __str__(self):       def __str__(self):
353         """         """
354         Returns the coordinates of the Locator as a string.         Returns the coordinates of the Locator as a string.
355         """         """
356         return "<Locator %s>"%str(self.getX())         x=self.getX()
357           if instance(x,list):
358              out="["
359              first=True
360              for xx in x:
361                if not first:
362                    out+=","
363                else:
364                    first=False
365                out+=str(xx)
366              out+="]>"
367           else:
368              out=str(x)
369           return out
370    
371         def getX(self):
372            """
373            Returns the exact coordinates of the Locator.
374            """
375            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):       def getId(self,item=None):
384          """          """
385      Returns the identifier of the location.          Returns the identifier of the location.
     """  
         return self.__id  
   
      def getX(self):  
386          """          """
387      Returns the exact coordinates of the Locator.          if item == None:
388      """             return self.__id
389          return self(self.getFunctionSpace().getX())          else:
390               if isinstance(self.__id,list):
391                  return self.__id[item]
392               else:
393                  return self.__id
394    
395    
396       def __call__(self,data):       def __call__(self,data):
397          """          """
398      Returns the value of data at the Locator of a Data object otherwise          Returns the value of data at the Locator of a Data object.
399      the object is returned.          """
     """  
400          return self.getValue(data)          return self.getValue(data)
401    
402       def getValue(self,data):       def getValue(self,data):
403          """          """
404      Returns the value of data at the Locator if data is a Data object          Returns the value of C{data} at the Locator if C{data} is a L{Data}
405      otherwise the object is returned.          object otherwise the object is returned.
406      """          """
407          if isinstance(data,escript.Data):          if isinstance(data,escript.Data):
408             if data.getFunctionSpace()==self.getFunctionSpace():             if data.getFunctionSpace()==self.getFunctionSpace():
409               out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])               dat=data
              #out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])  
410             else:             else:
411               out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])               dat=data.interpolate(self.getFunctionSpace())
412               #out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])             id=self.getId()
413             if data.getRank()==0:             r=data.getRank()
414                return out[0]             if isinstance(id,list):
415                   out=[]
416                   for i in id:
417                      o=numpy.array(data.getTupleForGlobalDataPoint(*i))
418                      if data.getRank()==0:
419                         out.append(o[0])
420                      else:
421                         out.append(o)
422                   return out
423             else:             else:
424                return out               out=numpy.array(data.getTupleForGlobalDataPoint(*id))
425                 if data.getRank()==0:
426                    return out[0]
427                 else:
428                    return out
429          else:          else:
430             return data             return data
431    
432  # vim: expandtab shiftwidth=4:  class SolverSchemeException(Exception):
433       """
434       This is a generic exception thrown by solvers.
435       """
436       pass
437    
438    class IndefinitePreconditioner(SolverSchemeException):
439       """
440       Exception thrown if the preconditioner is not positive definite.
441       """
442       pass
443    
444    class MaxIterReached(SolverSchemeException):
445       """
446       Exception thrown if the maximum number of iteration steps is reached.
447       """
448       pass
449    
450    class CorrectionFailed(SolverSchemeException):
451       """
452       Exception thrown if no convergence has been achieved in the solution
453       correction scheme.
454       """
455       pass
456    
457    class IterationBreakDown(SolverSchemeException):
458       """
459       Exception thrown if the iteration scheme encountered an incurable breakdown.
460       """
461       pass
462    
463    class NegativeNorm(SolverSchemeException):
464       """
465       Exception thrown if a norm calculation returns a negative norm.
466       """
467       pass
468    
469    def PCG(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100, initial_guess=True, verbose=False):
470       """
471       Solver for
472    
473       M{Ax=b}
474    
475       with a symmetric and positive definite operator A (more details required!).
476       It uses the conjugate gradient method with preconditioner M providing an
477       approximation of A.
478    
479       The iteration is terminated if
480    
481       M{|r| <= atol+rtol*|r0|}
482    
483       where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
484    
485       M{|r| = sqrt( bilinearform(Msolve(r),r))}
486    
487       For details on the preconditioned conjugate gradient method see the book:
488    
489       I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
490       T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
491       C. Romine, and H. van der Vorst}.
492    
493       @param r: initial residual M{r=b-Ax}. C{r} is altered.
494       @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
495       @param x: an initial guess for the solution
496       @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
497       @param Aprod: returns the value Ax
498       @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
499                    argument C{x}. The returned object needs to be of the same type
500                    like argument C{r}.
501       @param Msolve: solves Mx=r
502       @type Msolve: function C{Msolve(r)} where C{r} is of the same type like
503                     argument C{r}. The returned object needs to be of the same
504                     type like argument C{x}.
505       @param bilinearform: inner product C{<x,r>}
506       @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
507                           type like argument C{x} and C{r} is. The returned value
508                           is a C{float}.
509       @param atol: absolute tolerance
510       @type atol: non-negative C{float}
511       @param rtol: relative tolerance
512       @type rtol: non-negative C{float}
513       @param iter_max: maximum number of iteration steps
514       @type iter_max: C{int}
515       @return: the solution approximation and the corresponding residual
516       @rtype: C{tuple}
517       @warning: C{r} and C{x} are altered.
518       """
519       iter=0
520       rhat=Msolve(r)
521       d = rhat
522       rhat_dot_r = bilinearform(rhat, r)
523       if rhat_dot_r<0: raise NegativeNorm,"negative norm."
524       norm_r0=math.sqrt(rhat_dot_r)
525       atol2=atol+rtol*norm_r0
526       if atol2<=0:
527          raise ValueError,"Non-positive tolarance."
528       atol2=max(atol2, 100. * util.EPSILON * norm_r0)
529    
530       if verbose: print "PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_r0, atol2)
531    
532    
533       while not math.sqrt(rhat_dot_r) <= atol2:
534           iter+=1
535           if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
536    
537           q=Aprod(d)
538           alpha = rhat_dot_r / bilinearform(d, q)
539           x += alpha * d
540           if isinstance(q,ArithmeticTuple):
541           r += q * (-alpha)      # Doing it the other way calls the float64.__mul__ not AT.__rmul__
542           else:
543               r += (-alpha) * q
544           rhat=Msolve(r)
545           rhat_dot_r_new = bilinearform(rhat, r)
546           beta = rhat_dot_r_new / rhat_dot_r
547           rhat+=beta * d
548           d=rhat
549    
550           rhat_dot_r = rhat_dot_r_new
551           if rhat_dot_r<0: raise NegativeNorm,"negative norm."
552           if verbose: print "PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))
553       if verbose: print "PCG: tolerance reached after %s steps."%iter
554       return x,r,math.sqrt(rhat_dot_r)
555    
556    class Defect(object):
557        """
558        Defines a non-linear defect F(x) of a variable x.
559        """
560        def __init__(self):
561            """
562            Initializes defect.
563            """
564            self.setDerivativeIncrementLength()
565    
566        def bilinearform(self, x0, x1):
567            """
568            Returns the inner product of x0 and x1
569    
570            @param x0: value for x0
571            @param x1: value for x1
572            @return: the inner product of x0 and x1
573            @rtype: C{float}
574            """
575            return 0
576    
577        def norm(self,x):
578            """
579            Returns the norm of argument C{x}.
580    
581            @param x: a value
582            @return: norm of argument x
583            @rtype: C{float}
584            @note: by default C{sqrt(self.bilinearform(x,x)} is returned.
585            """
586            s=self.bilinearform(x,x)
587            if s<0: raise NegativeNorm,"negative norm."
588            return math.sqrt(s)
589    
590        def eval(self,x):
591            """
592            Returns the value F of a given C{x}.
593    
594            @param x: value for which the defect C{F} is evaluated
595            @return: value of the defect at C{x}
596            """
597            return 0
598    
599        def __call__(self,x):
600            return self.eval(x)
601    
602        def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):
603            """
604            Sets the relative length of the increment used to approximate the
605            derivative of the defect. The increment is inc*norm(x)/norm(v)*v in the
606            direction of v with x as a starting point.
607    
608            @param inc: relative increment length
609            @type inc: positive C{float}
610            """
611            if inc<=0: raise ValueError,"positive increment required."
612            self.__inc=inc
613    
614        def getDerivativeIncrementLength(self):
615            """
616            Returns the relative increment length used to approximate the
617            derivative of the defect.
618            @return: value of the defect at C{x}
619            @rtype: positive C{float}
620            """
621            return self.__inc
622    
623        def derivative(self, F0, x0, v, v_is_normalised=True):
624            """
625            Returns the directional derivative at C{x0} in the direction of C{v}.
626    
627            @param F0: value of this defect at x0
628            @param x0: value at which derivative is calculated
629            @param v: direction
630            @param v_is_normalised: True to indicate that C{v} is nomalized
631                                    (self.norm(v)=0)
632            @return: derivative of this defect at x0 in the direction of C{v}
633            @note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is
634                   used but this method maybe overwritten to use exact evaluation.
635            """
636            normx=self.norm(x0)
637            if normx>0:
638                 epsnew = self.getDerivativeIncrementLength() * normx
639            else:
640                 epsnew = self.getDerivativeIncrementLength()
641            if not v_is_normalised:
642                normv=self.norm(v)
643                if normv<=0:
644                   return F0*0
645                else:
646                   epsnew /= normv
647            F1=self.eval(x0 + epsnew * v)
648            return (F1-F0)/epsnew
649    
650    ######################################
651    def NewtonGMRES(defect, x, iter_max=100, sub_iter_max=20, atol=0,rtol=1.e-4, sub_tol_max=0.5, gamma=0.9, verbose=False):
652       """
653       Solves a non-linear problem M{F(x)=0} for unknown M{x} using the stopping
654       criterion:
655    
656       M{norm(F(x) <= atol + rtol * norm(F(x0)}
657    
658       where M{x0} is the initial guess.
659    
660       @param defect: object defining the function M{F}. C{defect.norm} defines the
661                      M{norm} used in the stopping criterion.
662       @type defect: L{Defect}
663       @param x: initial guess for the solution, C{x} is altered.
664       @type x: any object type allowing basic operations such as
665                C{numpy.ndarray}, L{Data}
666       @param iter_max: maximum number of iteration steps
667       @type iter_max: positive C{int}
668       @param sub_iter_max: maximum number of inner iteration steps
669       @type sub_iter_max: positive C{int}
670       @param atol: absolute tolerance for the solution
671       @type atol: positive C{float}
672       @param rtol: relative tolerance for the solution
673       @type rtol: positive C{float}
674       @param gamma: tolerance safety factor for inner iteration
675       @type gamma: positive C{float}, less than 1
676       @param sub_tol_max: upper bound for inner tolerance
677       @type sub_tol_max: positive C{float}, less than 1
678       @return: an approximation of the solution with the desired accuracy
679       @rtype: same type as the initial guess
680       """
681       lmaxit=iter_max
682       if atol<0: raise ValueError,"atol needs to be non-negative."
683       if rtol<0: raise ValueError,"rtol needs to be non-negative."
684       if rtol+atol<=0: raise ValueError,"rtol or atol needs to be non-negative."
685       if gamma<=0 or gamma>=1: raise ValueError,"tolerance safety factor for inner iteration (gamma =%s) needs to be positive and less than 1."%gamma
686       if sub_tol_max<=0 or sub_tol_max>=1: raise ValueError,"upper bound for inner tolerance for inner iteration (sub_tol_max =%s) needs to be positive and less than 1."%sub_tol_max
687    
688       F=defect(x)
689       fnrm=defect.norm(F)
690       stop_tol=atol + rtol*fnrm
691       sub_tol=sub_tol_max
692       if verbose: print "NewtonGMRES: initial residual = %e."%fnrm
693       if verbose: print "             tolerance = %e."%sub_tol
694       iter=1
695       #
696       # main iteration loop
697       #
698       while not fnrm<=stop_tol:
699                if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
700                #
701            #   adjust sub_tol_
702            #
703                if iter > 1:
704               rat=fnrm/fnrmo
705                   sub_tol_old=sub_tol
706               sub_tol=gamma*rat**2
707               if gamma*sub_tol_old**2 > .1: sub_tol=max(sub_tol,gamma*sub_tol_old**2)
708               sub_tol=max(min(sub_tol,sub_tol_max), .5*stop_tol/fnrm)
709            #
710            # calculate newton increment xc
711                #     if iter_max in __FDGMRES is reached MaxIterReached is thrown
712                #     if iter_restart -1 is returned as sub_iter
713                #     if  atol is reached sub_iter returns the numer of steps performed to get there
714                #
715                #
716                if verbose: print "             subiteration (GMRES) is called with relative tolerance %e."%sub_tol
717                try:
718                   xc, sub_iter=__FDGMRES(F, defect, x, sub_tol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)
719                except MaxIterReached:
720                   raise MaxIterReached,"maximum number of %s steps reached."%iter_max
721                if sub_iter<0:
722                   iter+=sub_iter_max
723                else:
724                   iter+=sub_iter
725                # ====
726            x+=xc
727                F=defect(x)
728            iter+=1
729                fnrmo, fnrm=fnrm, defect.norm(F)
730                if verbose: print "             step %s: residual %e."%(iter,fnrm)
731       if verbose: print "NewtonGMRES: completed after %s steps."%iter
732       return x
733    
734    def __givapp(c,s,vin):
735        """
736        Applies a sequence of Givens rotations (c,s) recursively to the vector
737        C{vin}
738    
739        @warning: C{vin} is altered.
740        """
741        vrot=vin
742        if isinstance(c,float):
743            vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
744        else:
745            for i in range(len(c)):
746                w1=c[i]*vrot[i]-s[i]*vrot[i+1]
747            w2=s[i]*vrot[i]+c[i]*vrot[i+1]
748                vrot[i]=w1
749                vrot[i+1]=w2
750        return vrot
751    
752    def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
753       h=numpy.zeros((iter_restart,iter_restart),numpy.float64)
754       c=numpy.zeros(iter_restart,numpy.float64)
755       s=numpy.zeros(iter_restart,numpy.float64)
756       g=numpy.zeros(iter_restart,numpy.float64)
757       v=[]
758    
759       rho=defect.norm(F0)
760       if rho<=0.: return x0*0
761    
762       v.append(-F0/rho)
763       g[0]=rho
764       iter=0
765       while rho > atol and iter<iter_restart-1:
766            if iter  >= iter_max:
767                raise MaxIterReached,"maximum number of %s steps reached."%iter_max
768    
769            p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
770            v.append(p)
771    
772            v_norm1=defect.norm(v[iter+1])
773    
774            # Modified Gram-Schmidt
775            for j in range(iter+1):
776                h[j,iter]=defect.bilinearform(v[j],v[iter+1])
777                v[iter+1]-=h[j,iter]*v[j]
778    
779            h[iter+1,iter]=defect.norm(v[iter+1])
780            v_norm2=h[iter+1,iter]
781    
782            # Reorthogonalize if needed
783            if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
784                for j in range(iter+1):
785                    hr=defect.bilinearform(v[j],v[iter+1])
786                    h[j,iter]=h[j,iter]+hr
787                    v[iter+1] -= hr*v[j]
788    
789                v_norm2=defect.norm(v[iter+1])
790                h[iter+1,iter]=v_norm2
791            #   watch out for happy breakdown
792            if not v_norm2 == 0:
793                v[iter+1]=v[iter+1]/h[iter+1,iter]
794    
795            #   Form and store the information for the new Givens rotation
796            if iter > 0 :
797                hhat=numpy.zeros(iter+1,numpy.float64)
798                for i in range(iter+1) : hhat[i]=h[i,iter]
799                hhat=__givapp(c[0:iter],s[0:iter],hhat);
800                for i in range(iter+1) : h[i,iter]=hhat[i]
801    
802            mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
803    
804            if mu!=0 :
805                c[iter]=h[iter,iter]/mu
806                s[iter]=-h[iter+1,iter]/mu
807                h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
808                h[iter+1,iter]=0.0
809                gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
810                g[iter]=gg[0]
811                g[iter+1]=gg[1]
812    
813            # Update the residual norm
814            rho=abs(g[iter+1])
815            iter+=1
816    
817       # At this point either iter > iter_max or rho < tol.
818       # It's time to compute x and leave.
819       if iter > 0 :
820         y=numpy.zeros(iter,numpy.float64)
821         y[iter-1] = g[iter-1] / h[iter-1,iter-1]
822         if iter > 1 :
823            i=iter-2
824            while i>=0 :
825              y[i] = ( g[i] - numpy.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
826              i=i-1
827         xhat=v[iter-1]*y[iter-1]
828         for i in range(iter-1):
829        xhat += v[i]*y[i]
830       else :
831          xhat=v[0] * 0
832    
833       if iter<iter_restart-1:
834          stopped=iter
835       else:
836          stopped=-1
837    
838       return xhat,stopped
839    
840    def GMRES(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100, iter_restart=20, verbose=False,P_R=None):
841       """
842       Solver for
843    
844       M{Ax=b}
845    
846       with a general operator A (more details required!).
847       It uses the generalized minimum residual method (GMRES).
848    
849       The iteration is terminated if
850    
851       M{|r| <= atol+rtol*|r0|}
852    
853       where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
854    
855       M{|r| = sqrt( bilinearform(r,r))}
856    
857       @param r: initial residual M{r=b-Ax}. C{r} is altered.
858       @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
859       @param x: an initial guess for the solution
860       @type x: same like C{r}
861       @param Aprod: returns the value Ax
862       @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
863                    argument C{x}. The returned object needs to be of the same
864                    type like argument C{r}.
865       @param bilinearform: inner product C{<x,r>}
866       @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
867                           type like argument C{x} and C{r}. The returned value is
868                           a C{float}.
869       @param atol: absolute tolerance
870       @type atol: non-negative C{float}
871       @param rtol: relative tolerance
872       @type rtol: non-negative C{float}
873       @param iter_max: maximum number of iteration steps
874       @type iter_max: C{int}
875       @param iter_restart: in order to save memory the orthogonalization process
876                            is terminated after C{iter_restart} steps and the
877                            iteration is restarted.
878       @type iter_restart: C{int}
879       @return: the solution approximation and the corresponding residual
880       @rtype: C{tuple}
881       @warning: C{r} and C{x} are altered.
882       """
883       m=iter_restart
884       restarted=False
885       iter=0
886       if rtol>0:
887          r_dot_r = bilinearform(r, r)
888          if r_dot_r<0: raise NegativeNorm,"negative norm."
889          atol2=atol+rtol*math.sqrt(r_dot_r)
890          if verbose: print "GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)
891       else:
892          atol2=atol
893          if verbose: print "GMRES: absolute tolerance = %e"%atol2
894       if atol2<=0:
895          raise ValueError,"Non-positive tolarance."
896    
897       while True:
898          if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
899          if restarted:
900             r2 = r-Aprod(x-x2)
901          else:
902             r2=1*r
903          x2=x*1.
904          x,stopped=_GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose,P_R=P_R)
905          iter+=iter_restart
906          if stopped: break
907          if verbose: print "GMRES: restart."
908          restarted=True
909       if verbose: print "GMRES: tolerance has been reached."
910       return x
911    
912    def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False, P_R=None):
913       iter=0
914    
915       h=numpy.zeros((iter_restart+1,iter_restart),numpy.float64)
916       c=numpy.zeros(iter_restart,numpy.float64)
917       s=numpy.zeros(iter_restart,numpy.float64)
918       g=numpy.zeros(iter_restart+1,numpy.float64)
919       v=[]
920    
921       r_dot_r = bilinearform(r, r)
922       if r_dot_r<0: raise NegativeNorm,"negative norm."
923       rho=math.sqrt(r_dot_r)
924    
925       v.append(r/rho)
926       g[0]=rho
927    
928       if verbose: print "GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)
929       while not (rho<=atol or iter==iter_restart):
930    
931        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
932    
933            if P_R!=None:
934                p=Aprod(P_R(v[iter]))
935            else:
936            p=Aprod(v[iter])
937        v.append(p)
938    
939        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
940    
941    # Modified Gram-Schmidt
942        for j in range(iter+1):
943          h[j,iter]=bilinearform(v[j],v[iter+1])
944          v[iter+1]-=h[j,iter]*v[j]
945    
946        h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
947        v_norm2=h[iter+1,iter]
948    
949    # Reorthogonalize if needed
950        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
951         for j in range(iter+1):
952            hr=bilinearform(v[j],v[iter+1])
953                h[j,iter]=h[j,iter]+hr
954                v[iter+1] -= hr*v[j]
955    
956         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
957         h[iter+1,iter]=v_norm2
958    
959    #   watch out for happy breakdown
960            if not v_norm2 == 0:
961             v[iter+1]=v[iter+1]/h[iter+1,iter]
962    
963    #   Form and store the information for the new Givens rotation
964        if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])
965        mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
966    
967        if mu!=0 :
968            c[iter]=h[iter,iter]/mu
969            s[iter]=-h[iter+1,iter]/mu
970            h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
971            h[iter+1,iter]=0.0
972                    gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
973                    g[iter]=gg[0]
974                    g[iter+1]=gg[1]
975    # Update the residual norm
976    
977            rho=abs(g[iter+1])
978            if verbose: print "GMRES: iteration step %s: residual %e"%(iter,rho)
979        iter+=1
980    
981    # At this point either iter > iter_max or rho < tol.
982    # It's time to compute x and leave.
983    
984       if verbose: print "GMRES: iteration stopped after %s step."%iter
985       if iter > 0 :
986         y=numpy.zeros(iter,numpy.float64)
987         y[iter-1] = g[iter-1] / h[iter-1,iter-1]
988         if iter > 1 :
989            i=iter-2
990            while i>=0 :
991              y[i] = ( g[i] - numpy.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
992              i=i-1
993         xhat=v[iter-1]*y[iter-1]
994         for i in range(iter-1):
995        xhat += v[i]*y[i]
996       else:
997         xhat=v[0] * 0
998       if P_R!=None:
999          x += P_R(xhat)
1000       else:
1001          x += xhat
1002       if iter<iter_restart-1:
1003          stopped=True
1004       else:
1005          stopped=False
1006    
1007       return x,stopped
1008    
1009    def MINRES(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1010        """
1011        Solver for
1012    
1013        M{Ax=b}
1014    
1015        with a symmetric and positive definite operator A (more details required!).
1016        It uses the minimum residual method (MINRES) with preconditioner M
1017        providing an approximation of A.
1018    
1019        The iteration is terminated if
1020    
1021        M{|r| <= atol+rtol*|r0|}
1022    
1023        where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
1024    
1025        M{|r| = sqrt( bilinearform(Msolve(r),r))}
1026    
1027        For details on the preconditioned conjugate gradient method see the book:
1028    
1029        I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
1030        T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
1031        C. Romine, and H. van der Vorst}.
1032    
1033        @param r: initial residual M{r=b-Ax}. C{r} is altered.
1034        @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1035        @param x: an initial guess for the solution
1036        @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1037        @param Aprod: returns the value Ax
1038        @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
1039                     argument C{x}. The returned object needs to be of the same
1040                     type like argument C{r}.
1041        @param Msolve: solves Mx=r
1042        @type Msolve: function C{Msolve(r)} where C{r} is of the same type like
1043                      argument C{r}. The returned object needs to be of the same
1044                      type like argument C{x}.
1045        @param bilinearform: inner product C{<x,r>}
1046        @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
1047                            type like argument C{x} and C{r} is. The returned value
1048                            is a C{float}.
1049        @param atol: absolute tolerance
1050        @type atol: non-negative C{float}
1051        @param rtol: relative tolerance
1052        @type rtol: non-negative C{float}
1053        @param iter_max: maximum number of iteration steps
1054        @type iter_max: C{int}
1055        @return: the solution approximation and the corresponding residual
1056        @rtype: C{tuple}
1057        @warning: C{r} and C{x} are altered.
1058        """
1059        #------------------------------------------------------------------
1060        # Set up y and v for the first Lanczos vector v1.
1061        # y  =  beta1 P' v1,  where  P = C**(-1).
1062        # v is really P' v1.
1063        #------------------------------------------------------------------
1064        r1    = r
1065        y = Msolve(r)
1066        beta1 = bilinearform(y,r)
1067    
1068        if beta1< 0: raise NegativeNorm,"negative norm."
1069    
1070        #  If r = 0 exactly, stop with x
1071        if beta1==0: return x
1072    
1073        if beta1> 0: beta1  = math.sqrt(beta1)
1074    
1075        #------------------------------------------------------------------
1076        # Initialize quantities.
1077        # ------------------------------------------------------------------
1078        iter   = 0
1079        Anorm = 0
1080        ynorm = 0
1081        oldb   = 0
1082        beta   = beta1
1083        dbar   = 0
1084        epsln  = 0
1085        phibar = beta1
1086        rhs1   = beta1
1087        rhs2   = 0
1088        rnorm  = phibar
1089        tnorm2 = 0
1090        ynorm2 = 0
1091        cs     = -1
1092        sn     = 0
1093        w      = r*0.
1094        w2     = r*0.
1095        r2     = r1
1096        eps    = 0.0001
1097    
1098        #---------------------------------------------------------------------
1099        # Main iteration loop.
1100        # --------------------------------------------------------------------
1101        while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL
1102    
1103        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1104            iter    = iter  +  1
1105    
1106            #-----------------------------------------------------------------
1107            # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
1108            # The general iteration is similar to the case k = 1 with v0 = 0:
1109            #
1110            #   p1      = Operator * v1  -  beta1 * v0,
1111            #   alpha1  = v1'p1,
1112            #   q2      = p2  -  alpha1 * v1,
1113            #   beta2^2 = q2'q2,
1114            #   v2      = (1/beta2) q2.
1115            #
1116            # Again, y = betak P vk,  where  P = C**(-1).
1117            #-----------------------------------------------------------------
1118            s = 1/beta                 # Normalize previous vector (in y).
1119            v = s*y                    # v = vk if P = I
1120    
1121            y      = Aprod(v)
1122    
1123            if iter >= 2:
1124              y = y - (beta/oldb)*r1
1125    
1126            alfa   = bilinearform(v,y)              # alphak
1127            y      += (- alfa/beta)*r2
1128            r1     = r2
1129            r2     = y
1130            y = Msolve(r2)
1131            oldb   = beta                           # oldb = betak
1132            beta   = bilinearform(y,r2)             # beta = betak+1^2
1133            if beta < 0: raise NegativeNorm,"negative norm."
1134    
1135            beta   = math.sqrt( beta )
1136            tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
1137    
1138            if iter==1:                 # Initialize a few things.
1139              gmax   = abs( alfa )      # alpha1
1140              gmin   = gmax             # alpha1
1141    
1142            # Apply previous rotation Qk-1 to get
1143            #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
1144            #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
1145    
1146            oldeps = epsln
1147            delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
1148            gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
1149            epsln  =               sn * beta  # epsln2 = 0         epslnk+1
1150            dbar   =            -  cs * beta  # dbar 2 = beta2     dbar k+1
1151    
1152            # Compute the next plane rotation Qk
1153    
1154            gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
1155            gamma  = max(gamma,eps)
1156            cs     = gbar / gamma             # ck
1157            sn     = beta / gamma             # sk
1158            phi    = cs * phibar              # phik
1159            phibar = sn * phibar              # phibark+1
1160    
1161            # Update  x.
1162    
1163            denom = 1/gamma
1164            w1    = w2
1165            w2    = w
1166            w     = (v - oldeps*w1 - delta*w2) * denom
1167            x     +=  phi*w
1168    
1169            # Go round again.
1170    
1171            gmax   = max(gmax,gamma)
1172            gmin   = min(gmin,gamma)
1173            z      = rhs1 / gamma
1174            ynorm2 = z*z  +  ynorm2
1175            rhs1   = rhs2 -  delta*z
1176            rhs2   =      -  epsln*z
1177    
1178            # Estimate various norms and test for convergence.
1179    
1180            Anorm  = math.sqrt( tnorm2 )
1181            ynorm  = math.sqrt( ynorm2 )
1182    
1183            rnorm  = phibar
1184    
1185        return x
1186    
1187    def TFQMR(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1188      """
1189      Solver for
1190    
1191      M{Ax=b}
1192    
1193      with a general operator A (more details required!).
1194      It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).
1195    
1196      The iteration is terminated if
1197    
1198      M{|r| <= atol+rtol*|r0|}
1199    
1200      where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
1201    
1202      M{|r| = sqrt( bilinearform(r,r))}
1203    
1204      @param r: initial residual M{r=b-Ax}. C{r} is altered.
1205      @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1206      @param x: an initial guess for the solution
1207      @type x: same like C{r}
1208      @param Aprod: returns the value Ax
1209      @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
1210                   argument C{x}. The returned object needs to be of the same type
1211                   like argument C{r}.
1212      @param bilinearform: inner product C{<x,r>}
1213      @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
1214                          type like argument C{x} and C{r}. The returned value is
1215                          a C{float}.
1216      @param atol: absolute tolerance
1217      @type atol: non-negative C{float}
1218      @param rtol: relative tolerance
1219      @type rtol: non-negative C{float}
1220      @param iter_max: maximum number of iteration steps
1221      @type iter_max: C{int}
1222      @rtype: C{tuple}
1223      @warning: C{r} and C{x} are altered.
1224      """
1225      u1=0
1226      u2=0
1227      y1=0
1228      y2=0
1229    
1230      w = r
1231      y1 = r
1232      iter = 0
1233      d = 0
1234      v = Aprod(y1)
1235      u1 = v
1236    
1237      theta = 0.0;
1238      eta = 0.0;
1239      rho=bilinearform(r,r)
1240      if rho < 0: raise NegativeNorm,"negative norm."
1241      tau = math.sqrt(rho)
1242      norm_r0=tau
1243      while tau>atol+rtol*norm_r0:
1244        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1245    
1246        sigma = bilinearform(r,v)
1247        if sigma == 0.0: raise IterationBreakDown,'TFQMR breakdown, sigma=0'
1248    
1249        alpha = rho / sigma
1250    
1251        for j in range(2):
1252    #
1253    #   Compute y2 and u2 only if you have to
1254    #
1255          if ( j == 1 ):
1256            y2 = y1 - alpha * v
1257            u2 = Aprod(y2)
1258    
1259          m = 2 * (iter+1) - 2 + (j+1)
1260          if j==0:
1261             w = w - alpha * u1
1262             d = y1 + ( theta * theta * eta / alpha ) * d
1263          if j==1:
1264             w = w - alpha * u2
1265             d = y2 + ( theta * theta * eta / alpha ) * d
1266    
1267          theta = math.sqrt(bilinearform(w,w))/ tau
1268          c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1269          tau = tau * theta * c
1270          eta = c * c * alpha
1271          x = x + eta * d
1272    #
1273    #  Try to terminate the iteration at each pass through the loop
1274    #
1275        if rho == 0.0: raise IterationBreakDown,'TFQMR breakdown, rho=0'
1276    
1277        rhon = bilinearform(r,w)
1278        beta = rhon / rho;
1279        rho = rhon;
1280        y1 = w + beta * y2;
1281        u1 = Aprod(y1)
1282        v = u1 + beta * ( u2 + beta * v )
1283    
1284        iter += 1
1285    
1286      return x
1287    
1288    
1289    #############################################
1290    
1291    class ArithmeticTuple(object):
1292       """
1293       Tuple supporting inplace update x+=y and scaling x=a*y where C{x,y} is an
1294       ArithmeticTuple and C{a} is a float.
1295    
1296       Example of usage::
1297    
1298           from esys.escript import Data
1299           from numpy import array
1300           a=Data(...)
1301           b=array([1.,4.])
1302           x=ArithmeticTuple(a,b)
1303           y=5.*x
1304    
1305       """
1306       def __init__(self,*args):
1307           """
1308           Initializes object with elements C{args}.
1309    
1310           @param args: tuple of objects that support inplace add (x+=y) and
1311                        scaling (x=a*y)
1312           """
1313           self.__items=list(args)
1314    
1315       def __len__(self):
1316           """
1317           Returns the number of items.
1318    
1319           @return: number of items
1320           @rtype: C{int}
1321           """
1322           return len(self.__items)
1323    
1324       def __getitem__(self,index):
1325           """
1326           Returns item at specified position.
1327    
1328           @param index: index of item to be returned
1329           @type index: C{int}
1330           @return: item with index C{index}
1331           """
1332           return self.__items.__getitem__(index)
1333    
1334       def __mul__(self,other):
1335           """
1336           Scales by C{other} from the right.
1337    
1338           @param other: scaling factor
1339           @type other: C{float}
1340           @return: itemwise self*other
1341           @rtype: L{ArithmeticTuple}
1342           """
1343           out=[]
1344           try:
1345               l=len(other)
1346               if l!=len(self):
1347                   raise ValueError,"length of arguments don't match."
1348               for i in range(l): out.append(self[i]*other[i])
1349           except TypeError:
1350               for i in range(len(self)): out.append(self[i]*other)
1351           return ArithmeticTuple(*tuple(out))
1352    
1353       def __rmul__(self,other):
1354           """
1355           Scales by C{other} from the left.
1356    
1357           @param other: scaling factor
1358           @type other: C{float}
1359           @return: itemwise other*self
1360           @rtype: L{ArithmeticTuple}
1361           """
1362           out=[]
1363           try:
1364               l=len(other)
1365               if l!=len(self):
1366                   raise ValueError,"length of arguments don't match."
1367               for i in range(l): out.append(other[i]*self[i])
1368           except TypeError:
1369               for i in range(len(self)): out.append(other*self[i])
1370           return ArithmeticTuple(*tuple(out))
1371    
1372       def __div__(self,other):
1373           """
1374           Scales by (1/C{other}) from the right.
1375    
1376           @param other: scaling factor
1377           @type other: C{float}
1378           @return: itemwise self/other
1379           @rtype: L{ArithmeticTuple}
1380           """
1381           return self*(1/other)
1382    
1383       def __rdiv__(self,other):
1384           """
1385           Scales by (1/C{other}) from the left.
1386    
1387           @param other: scaling factor
1388           @type other: C{float}
1389           @return: itemwise other/self
1390           @rtype: L{ArithmeticTuple}
1391           """
1392           out=[]
1393           try:
1394               l=len(other)
1395               if l!=len(self):
1396                   raise ValueError,"length of arguments don't match."
1397               for i in range(l): out.append(other[i]/self[i])
1398           except TypeError:
1399               for i in range(len(self)): out.append(other/self[i])
1400           return ArithmeticTuple(*tuple(out))
1401    
1402       def __iadd__(self,other):
1403           """
1404           Inplace addition of C{other} to self.
1405    
1406           @param other: increment
1407           @type other: C{ArithmeticTuple}
1408           """
1409           if len(self) != len(other):
1410               raise ValueError,"tuple lengths must match."
1411           for i in range(len(self)):
1412               self.__items[i]+=other[i]
1413           return self
1414    
1415       def __add__(self,other):
1416           """
1417           Adds C{other} to self.
1418    
1419           @param other: increment
1420           @type other: C{ArithmeticTuple}
1421           """
1422           out=[]
1423           try:
1424               l=len(other)
1425               if l!=len(self):
1426                   raise ValueError,"length of arguments don't match."
1427               for i in range(l): out.append(self[i]+other[i])
1428           except TypeError:
1429               for i in range(len(self)): out.append(self[i]+other)
1430           return ArithmeticTuple(*tuple(out))
1431    
1432       def __sub__(self,other):
1433           """
1434           Subtracts C{other} from self.
1435    
1436           @param other: decrement
1437           @type other: C{ArithmeticTuple}
1438           """
1439           out=[]
1440           try:
1441               l=len(other)
1442               if l!=len(self):
1443                   raise ValueError,"length of arguments don't match."
1444               for i in range(l): out.append(self[i]-other[i])
1445           except TypeError:
1446               for i in range(len(self)): out.append(self[i]-other)
1447           return ArithmeticTuple(*tuple(out))
1448    
1449       def __isub__(self,other):
1450           """
1451           Inplace subtraction of C{other} from self.
1452    
1453           @param other: decrement
1454           @type other: C{ArithmeticTuple}
1455           """
1456           if len(self) != len(other):
1457               raise ValueError,"tuple length must match."
1458           for i in range(len(self)):
1459               self.__items[i]-=other[i]
1460           return self
1461    
1462       def __neg__(self):
1463           """
1464           Negates values.
1465           """
1466           out=[]
1467           for i in range(len(self)):
1468               out.append(-self[i])
1469           return ArithmeticTuple(*tuple(out))
1470    
1471    
1472    class HomogeneousSaddlePointProblem(object):
1473          """
1474          This class provides a framework for solving linear homogeneous saddle
1475          point problems of the form::
1476    
1477              M{Av+B^*p=f}
1478              M{Bv     =0}
1479    
1480          for the unknowns M{v} and M{p} and given operators M{A} and M{B} and
1481          given right hand side M{f}. M{B^*} is the adjoint operator of M{B}.
1482          """
1483          def __init__(self, adaptSubTolerance=True, **kwargs):
1484        """
1485        initializes the saddle point problem
1486        
1487        @param adaptSubTolerance: If True the tolerance for subproblem is set automatically.
1488        @type adaptSubTolerance: C{bool}
1489        """
1490            self.setTolerance()
1491            self.setAbsoluteTolerance()
1492        self.__adaptSubTolerance=adaptSubTolerance
1493          #=============================================================
1494          def initialize(self):
1495            """
1496            Initializes the problem (overwrite).
1497            """
1498            pass
1499    
1500          def inner_pBv(self,p,Bv):
1501             """
1502             Returns inner product of element p and Bv (overwrite).
1503    
1504             @param p: a pressure increment
1505             @param v: a residual
1506             @return: inner product of element p and Bv
1507             @rtype: C{float}
1508             @note: used if PCG is applied.
1509             """
1510             raise NotImplementedError,"no inner product for p and Bv implemented."
1511    
1512          def inner_p(self,p0,p1):
1513             """
1514             Returns inner product of p0 and p1 (overwrite).
1515    
1516             @param p0: a pressure
1517             @param p1: a pressure
1518             @return: inner product of p0 and p1
1519             @rtype: C{float}
1520             """
1521             raise NotImplementedError,"no inner product for p implemented."
1522      
1523          def norm_v(self,v):
1524             """
1525             Returns the norm of v (overwrite).
1526    
1527             @param v: a velovity
1528             @return: norm of v
1529             @rtype: non-negative C{float}
1530             """
1531             raise NotImplementedError,"no norm of v implemented."
1532          def getV(self, p, v0):
1533             """
1534             return the value for v for a given p (overwrite)
1535    
1536             @param p: a pressure
1537             @param v0: a initial guess for the value v to return.
1538             @return: v given as M{v= A^{-1} (f-B^*p)}
1539             """
1540             raise NotImplementedError,"no v calculation implemented."
1541    
1542            
1543          def Bv(self,v):
1544            """
1545            Returns Bv (overwrite).
1546    
1547            @rtype: equal to the type of p
1548            @note: boundary conditions on p should be zero!
1549            """
1550            raise NotImplementedError, "no operator B implemented."
1551    
1552          def norm_Bv(self,Bv):
1553            """
1554            Returns the norm of Bv (overwrite).
1555    
1556            @rtype: equal to the type of p
1557            @note: boundary conditions on p should be zero!
1558            """
1559            raise NotImplementedError, "no norm of Bv implemented."
1560    
1561          def solve_AinvBt(self,p):
1562             """
1563             Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
1564             (overwrite).
1565    
1566             @param p: a pressure increment
1567             @return: the solution of M{Av=B^*p}
1568             @note: boundary conditions on v should be zero!
1569             """
1570             raise NotImplementedError,"no operator A implemented."
1571    
1572          def solve_prec(self,Bv):
1573             """
1574             Provides a preconditioner for M{BA^{-1}B^*} applied to Bv with accuracy
1575             L{self.getSubProblemTolerance()} (overwrite).
1576    
1577             @rtype: equal to the type of p
1578             @note: boundary conditions on p should be zero!
1579             """
1580             raise NotImplementedError,"no preconditioner for Schur complement implemented."
1581          def setSubProblemTolerance(self):
1582             """
1583         Updates the tolerance for subproblems
1584         @note: method is typically the method is overwritten.
1585             """
1586             pass
1587          #=============================================================
1588          def __Aprod_PCG(self,p):
1589              dv=self.solve_AinvBt(p)
1590              return ArithmeticTuple(dv,self.Bv(dv))
1591    
1592          def __inner_PCG(self,p,r):
1593             return self.inner_pBv(p,r[1])
1594    
1595          def __Msolve_PCG(self,r):
1596              return self.solve_prec(r[1])
1597          #=============================================================
1598          def __Aprod_GMRES(self,p):
1599              return self.solve_prec(self.Bv(self.solve_AinvBt(p)))
1600          def __inner_GMRES(self,p0,p1):
1601             return self.inner_p(p0,p1)
1602    
1603          #=============================================================
1604          def norm_p(self,p):
1605              """
1606              calculates the norm of C{p}
1607              
1608              @param p: a pressure
1609              @return: the norm of C{p} using the inner product for pressure
1610              @rtype: C{float}
1611              """
1612              f=self.inner_p(p,p)
1613              if f<0: raise ValueError,"negative pressure norm."
1614              return math.sqrt(f)
1615          def adaptSubTolerance(self):
1616          """
1617          Returns True if tolerance adaption for subproblem is choosen.
1618          """
1619              self.__adaptSubTolerance
1620          
1621          def solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1622             """
1623             Solves the saddle point problem using initial guesses v and p.
1624    
1625             @param v: initial guess for velocity
1626             @param p: initial guess for pressure
1627             @type v: L{Data}
1628             @type p: L{Data}
1629             @param usePCG: indicates the usage of the PCG rather than GMRES scheme.
1630             @param max_iter: maximum number of iteration steps per correction
1631                              attempt
1632             @param verbose: if True, shows information on the progress of the
1633                             saddlepoint problem solver.
1634             @param iter_restart: restart the iteration after C{iter_restart} steps
1635                                  (only used if useUzaw=False)
1636             @type usePCG: C{bool}
1637             @type max_iter: C{int}
1638             @type verbose: C{bool}
1639             @type iter_restart: C{int}
1640             @rtype: C{tuple} of L{Data} objects
1641             """
1642             self.verbose=verbose
1643             rtol=self.getTolerance()
1644             atol=self.getAbsoluteTolerance()
1645         if self.adaptSubTolerance(): self.setSubProblemTolerance()
1646             correction_step=0
1647             converged=False
1648             while not converged:
1649                  # calculate velocity for current pressure:
1650                  v=self.getV(p,v)
1651                  Bv=self.Bv(v)
1652                  norm_v=self.norm_v(v)
1653                  norm_Bv=self.norm_Bv(Bv)
1654                  ATOL=norm_v*rtol+atol
1655                  if self.verbose: print "HomogeneousSaddlePointProblem: norm v= %e, norm_Bv= %e, tolerance = %e."%(norm_v, norm_Bv,ATOL)
1656                  if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."
1657                  if norm_Bv <= ATOL:
1658                     converged=True
1659                  else:
1660                     correction_step+=1
1661                     if correction_step>max_correction_steps:
1662                          raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1663                     dp=self.solve_prec(Bv)
1664                     if usePCG:
1665                       norm2=self.inner_pBv(dp,Bv)
1666                       if norm2<0: raise ValueError,"negative PCG norm."
1667                       norm2=math.sqrt(norm2)
1668                     else:
1669                       norm2=self.norm_p(dp)
1670                     ATOL_ITER=ATOL/norm_Bv*norm2*0.5
1671                     if self.verbose: print "HomogeneousSaddlePointProblem: tolerance for solver: %e"%ATOL_ITER
1672                     if usePCG:
1673                           p,v0,a_norm=PCG(ArithmeticTuple(v,Bv),self.__Aprod_PCG,p,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL_ITER, rtol=0.,iter_max=max_iter, verbose=self.verbose)
1674                     else:
1675                           p=GMRES(dp,self.__Aprod_GMRES, p, self.__inner_GMRES,atol=ATOL_ITER, rtol=0.,iter_max=max_iter, iter_restart=iter_restart, verbose=self.verbose)
1676             if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached."
1677         return v,p
1678    
1679          #========================================================================
1680          def setTolerance(self,tolerance=1.e-4):
1681             """
1682             Sets the relative tolerance for (v,p).
1683    
1684             @param tolerance: tolerance to be used
1685             @type tolerance: non-negative C{float}
1686             """
1687             if tolerance<0:
1688                 raise ValueError,"tolerance must be positive."
1689             self.__rtol=tolerance
1690    
1691          def getTolerance(self):
1692             """
1693             Returns the relative tolerance.
1694    
1695             @return: relative tolerance
1696             @rtype: C{float}
1697             """
1698             return self.__rtol
1699    
1700          def setAbsoluteTolerance(self,tolerance=0.):
1701             """
1702             Sets the absolute tolerance.
1703    
1704             @param tolerance: tolerance to be used
1705             @type tolerance: non-negative C{float}
1706             """
1707             if tolerance<0:
1708                 raise ValueError,"tolerance must be non-negative."
1709             self.__atol=tolerance
1710    
1711          def getAbsoluteTolerance(self):
1712             """
1713             Returns the absolute tolerance.
1714    
1715             @return: absolute tolerance
1716             @rtype: C{float}
1717             """
1718             return self.__atol
1719    
1720          def getSubProblemTolerance(self):
1721             """
1722             Sets the relative tolerance to solve the subproblem(s).
1723    
1724             @param rtol: relative tolerence
1725             @type rtol: positive C{float}
1726             """
1727             return max(200.*util.EPSILON,self.getTolerance()**2)
1728    
1729    def MaskFromBoundaryTag(domain,*tags):
1730       """
1731       Creates a mask on the Solution(domain) function space where the value is
1732       one for samples that touch the boundary tagged by tags.
1733    
1734       Usage: m=MaskFromBoundaryTag(domain, "left", "right")
1735    
1736       @param domain: domain to be used
1737       @type domain: L{escript.Domain}
1738       @param tags: boundary tags
1739       @type tags: C{str}
1740       @return: a mask which marks samples that are touching the boundary tagged
1741                by any of the given tags
1742       @rtype: L{escript.Data} of rank 0
1743       """
1744       pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1745       d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
1746       for t in tags: d.setTaggedValue(t,1.)
1747       pde.setValue(y=d)
1748       return util.whereNonZero(pde.getRightHandSide())
1749    
1750    def MaskFromTag(domain,*tags):
1751       """
1752       Creates a mask on the Solution(domain) function space where the value is
1753       one for samples that touch regions tagged by tags.
1754    
1755       Usage: m=MaskFromTag(domain, "ham")
1756    
1757       @param domain: domain to be used
1758       @type domain: L{escript.Domain}
1759       @param tags: boundary tags
1760       @type tags: C{str}
1761       @return: a mask which marks samples that are touching the boundary tagged
1762                by any of the given tags
1763       @rtype: L{escript.Data} of rank 0
1764       """
1765       pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1766       d=escript.Scalar(0.,escript.Function(domain))
1767       for t in tags: d.setTaggedValue(t,1.)
1768       pde.setValue(Y=d)
1769       return util.whereNonZero(pde.getRightHandSide())
1770    
1771    

Legend:
Removed from v.782  
changed lines
  Added in v.2498

  ViewVC Help
Powered by ViewVC 1.1.26