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

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

  ViewVC Help
Powered by ViewVC 1.1.26