/[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 2440 by gross, Wed May 27 08:45:55 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 numarray
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.setSolverMethod(linearPDEs.LinearPDE.LUMPING)
149      self.__pde.setSymmetryOn()      self.__pde.setSymmetryOn()
150      self.__pde.setReducedOrderTo(reduce)      self.__pde.setReducedOrderTo(reduce)
151      self.__pde.setValue(D = 1.)      self.__pde.setValue(D = 1.)
152      return      return
153    
   def __del__(self):  
     return  
   
154    def __call__(self, input_data):    def __call__(self, input_data):
155      """      """
156      Projects input_data onto a continuous function      Projects C{input_data} onto a continuous function.
157    
158      @param input_data: The input_data to be projected.      @param input_data: the data to be projected
159      """      """
160      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
161        self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())
162      if input_data.getRank()==0:      if input_data.getRank()==0:
163          self.__pde.setValue(Y = input_data)          self.__pde.setValue(Y = input_data)
164          out=self.__pde.getSolution()          out=self.__pde.getSolution()
# Line 170  class Projector: Line 188  class Projector:
188    
189  class NoPDE:  class NoPDE:
190       """       """
191       solves the following problem for u:       Solves the following problem for u:
192    
193       M{kronecker[i,j]*D[j]*u[j]=Y[i]}       M{kronecker[i,j]*D[j]*u[j]=Y[i]}
194    
195       with constraint       with constraint
196    
197       M{u[j]=r[j]}  where M{q[j]>0}       M{u[j]=r[j]}  where M{q[j]>0}
198    
199       where D, Y, r and q are given functions of rank 1.       where M{D}, M{Y}, M{r} and M{q} are given functions of rank 1.
200    
201       In the case of scalars this takes the form       In the case of scalars this takes the form
202    
203       M{D*u=Y}       M{D*u=Y}
204    
205       with constraint       with constraint
206    
207       M{u=r}  where M{q>0}       M{u=r} where M{q>0}
208    
209       where D, Y, r and q are given scalar functions.       where M{D}, M{Y}, M{r} and M{q} are given scalar functions.
210    
211       The constraint is overwriting any other condition.       The constraint overwrites any other condition.
212    
213       @note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention       @note: This class is similar to the L{linearPDEs.LinearPDE} class with
214              that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole              A=B=C=X=0 but has the intention that all input parameters are given
215              thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.              in L{Solution} or L{ReducedSolution}.
216       """       """
217         # The whole thing is a bit strange and I blame Rob Woodcock (CSIRO) for
218         # this.
219       def __init__(self,domain,D=None,Y=None,q=None,r=None):       def __init__(self,domain,D=None,Y=None,q=None,r=None):
220           """           """
221           initialize the problem           Initializes the problem.
222    
223           @param domain: domain of the PDE.           @param domain: domain of the PDE
224           @type domain: L{Domain}           @type domain: L{Domain}
225           @param D: coefficient of the solution.           @param D: coefficient of the solution
226           @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}           @type D: C{float}, C{int}, C{numarray.NumArray}, L{Data}
227           @param Y: right hand side           @param Y: right hand side
228           @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}           @type Y: C{float}, C{int}, C{numarray.NumArray}, L{Data}
229           @param q: location of constraints           @param q: location of constraints
230           @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}           @type q: C{float}, C{int}, C{numarray.NumArray}, L{Data}
231           @param r: value of solution at locations of constraints           @param r: value of solution at locations of constraints
232           @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}           @type r: C{float}, C{int}, C{numarray.NumArray}, L{Data}
233           """           """
234           self.__domain=domain           self.__domain=domain
235           self.__D=D           self.__D=D
# Line 218  class NoPDE: Line 238  class NoPDE:
238           self.__r=r           self.__r=r
239           self.__u=None           self.__u=None
240           self.__function_space=escript.Solution(self.__domain)           self.__function_space=escript.Solution(self.__domain)
241    
242       def setReducedOn(self):       def setReducedOn(self):
243           """           """
244           sets the L{FunctionSpace} of the solution to L{ReducedSolution}           Sets the L{FunctionSpace} of the solution to L{ReducedSolution}.
245           """           """
246           self.__function_space=escript.ReducedSolution(self.__domain)           self.__function_space=escript.ReducedSolution(self.__domain)
247           self.__u=None           self.__u=None
248    
249       def setReducedOff(self):       def setReducedOff(self):
250           """           """
251           sets the L{FunctionSpace} of the solution to L{Solution}           Sets the L{FunctionSpace} of the solution to L{Solution}.
252           """           """
253           self.__function_space=escript.Solution(self.__domain)           self.__function_space=escript.Solution(self.__domain)
254           self.__u=None           self.__u=None
255            
256       def setValue(self,D=None,Y=None,q=None,r=None):       def setValue(self,D=None,Y=None,q=None,r=None):
257           """           """
258           assigns values to the parameters.           Assigns values to the parameters.
259    
260           @param D: coefficient of the solution.           @param D: coefficient of the solution
261           @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}           @type D: C{float}, C{int}, C{numarray.NumArray}, L{Data}
262           @param Y: right hand side           @param Y: right hand side
263           @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}           @type Y: C{float}, C{int}, C{numarray.NumArray}, L{Data}
264           @param q: location of constraints           @param q: location of constraints
265           @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}           @type q: C{float}, C{int}, C{numarray.NumArray}, L{Data}
266           @param r: value of solution at locations of constraints           @param r: value of solution at locations of constraints
267           @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}           @type r: C{float}, C{int}, C{numarray.NumArray}, L{Data}
268           """           """
269           if not D==None:           if not D==None:
270              self.__D=D              self.__D=D
# Line 260  class NoPDE: Line 281  class NoPDE:
281    
282       def getSolution(self):       def getSolution(self):
283           """           """
284           returns the solution           Returns the solution.
285            
286           @return: the solution of the problem           @return: the solution of the problem
287           @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.           @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or
288                     L{ReducedSolution}
289           """           """
290           if self.__u==None:           if self.__u==None:
291              if self.__D==None:              if self.__D==None:
# Line 280  class NoPDE: Line 302  class NoPDE:
302                  self.__u*=(1.-q)                  self.__u*=(1.-q)
303                  if not self.__r==None: self.__u+=q*self.__r                  if not self.__r==None: self.__u+=q*self.__r
304           return self.__u           return self.__u
305                
306  class Locator:  class Locator:
307       """       """
308       Locator provides access to the values of data objects at a given       Locator provides access to the values of data objects at a given spatial
309       spatial coordinate x.         coordinate x.
310        
311       In fact, a Locator object finds the sample in the set of samples of a       In fact, a Locator object finds the sample in the set of samples of a
312       given function space or domain where which is closest to the given       given function space or domain which is closest to the given point x.
      point x.  
313       """       """
314    
315       def __init__(self,where,x=numarray.zeros((3,))):       def __init__(self,where,x=numarray.zeros((3,))):
316         """         """
317         Initializes a Locator to access values in Data objects on the Doamin         Initializes a Locator to access values in Data objects on the Doamin
318         or FunctionSpace where for the sample point which         or FunctionSpace for the sample point which is closest to the given
319         closest to the given point x.         point x.
320    
321           @param where: function space
322           @type where: L{escript.FunctionSpace}
323           @param x: coefficient of the solution
324           @type x: C{numarray.NumArray} or C{list} of C{numarray.NumArray}
325         """         """
326         if isinstance(where,escript.FunctionSpace):         if isinstance(where,escript.FunctionSpace):
327            self.__function_space=where            self.__function_space=where
328         else:         else:
329            self.__function_space=escript.ContinuousFunction(where)            self.__function_space=escript.ContinuousFunction(where)
330         self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).mindp()         if isinstance(x, list):
331               self.__id=[]
332               for p in x:
333                  self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())
334           else:
335               self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint()
336    
337       def __str__(self):       def __str__(self):
338         """         """
339         Returns the coordinates of the Locator as a string.         Returns the coordinates of the Locator as a string.
340         """         """
341         return "<Locator %s>"%str(self.getX())         x=self.getX()
342           if instance(x,list):
343              out="["
344              first=True
345              for xx in x:
346                if not first:
347                    out+=","
348                else:
349                    first=False
350                out+=str(xx)
351              out+="]>"
352           else:
353              out=str(x)
354           return out
355    
356         def getX(self):
357            """
358            Returns the exact coordinates of the Locator.
359            """
360            return self(self.getFunctionSpace().getX())
361    
362       def getFunctionSpace(self):       def getFunctionSpace(self):
363          """          """
364      Returns the function space of the Locator.          Returns the function space of the Locator.
365      """          """
366          return self.__function_space          return self.__function_space
367    
368       def getId(self):       def getId(self,item=None):
369          """          """
370      Returns the identifier of the location.          Returns the identifier of the location.
     """  
         return self.__id  
   
      def getX(self):  
371          """          """
372      Returns the exact coordinates of the Locator.          if item == None:
373      """             return self.__id
374          return self(self.getFunctionSpace().getX())          else:
375               if isinstance(self.__id,list):
376                  return self.__id[item]
377               else:
378                  return self.__id
379    
380    
381       def __call__(self,data):       def __call__(self,data):
382          """          """
383      Returns the value of data at the Locator of a Data object otherwise          Returns the value of data at the Locator of a Data object.
384      the object is returned.          """
     """  
385          return self.getValue(data)          return self.getValue(data)
386    
387       def getValue(self,data):       def getValue(self,data):
388          """          """
389      Returns the value of data at the Locator if data is a Data object          Returns the value of C{data} at the Locator if C{data} is a L{Data}
390      otherwise the object is returned.          object otherwise the object is returned.
391      """          """
392          if isinstance(data,escript.Data):          if isinstance(data,escript.Data):
393             if data.getFunctionSpace()==self.getFunctionSpace():             if data.getFunctionSpace()==self.getFunctionSpace():
394               out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])               dat=data
              #out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])  
395             else:             else:
396               out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])               dat=data.interpolate(self.getFunctionSpace())
397               #out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1],self.getId()[2])             id=self.getId()
398             if data.getRank()==0:             r=data.getRank()
399                return out[0]             if isinstance(id,list):
400                   out=[]
401                   for i in id:
402                      o=data.getValueOfGlobalDataPoint(*i)
403                      if data.getRank()==0:
404                         out.append(o[0])
405                      else:
406                         out.append(o)
407                   return out
408             else:             else:
409                return out               out=data.getValueOfGlobalDataPoint(*id)
410                 if data.getRank()==0:
411                    return out[0]
412                 else:
413                    return out
414          else:          else:
415             return data             return data
416    
417  # vim: expandtab shiftwidth=4:  class SolverSchemeException(Exception):
418       """
419       This is a generic exception thrown by solvers.
420       """
421       pass
422    
423    class IndefinitePreconditioner(SolverSchemeException):
424       """
425       Exception thrown if the preconditioner is not positive definite.
426       """
427       pass
428    
429    class MaxIterReached(SolverSchemeException):
430       """
431       Exception thrown if the maximum number of iteration steps is reached.
432       """
433       pass
434    
435    class CorrectionFailed(SolverSchemeException):
436       """
437       Exception thrown if no convergence has been achieved in the solution
438       correction scheme.
439       """
440       pass
441    
442    class IterationBreakDown(SolverSchemeException):
443       """
444       Exception thrown if the iteration scheme encountered an incurable breakdown.
445       """
446       pass
447    
448    class NegativeNorm(SolverSchemeException):
449       """
450       Exception thrown if a norm calculation returns a negative norm.
451       """
452       pass
453    
454    def PCG(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100, initial_guess=True, verbose=False):
455       """
456       Solver for
457    
458       M{Ax=b}
459    
460       with a symmetric and positive definite operator A (more details required!).
461       It uses the conjugate gradient method with preconditioner M providing an
462       approximation of A.
463    
464       The iteration is terminated if
465    
466       M{|r| <= atol+rtol*|r0|}
467    
468       where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
469    
470       M{|r| = sqrt( bilinearform(Msolve(r),r))}
471    
472       For details on the preconditioned conjugate gradient method see the book:
473    
474       I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
475       T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
476       C. Romine, and H. van der Vorst}.
477    
478       @param r: initial residual M{r=b-Ax}. C{r} is altered.
479       @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
480       @param x: an initial guess for the solution
481       @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
482       @param Aprod: returns the value Ax
483       @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
484                    argument C{x}. The returned object needs to be of the same type
485                    like argument C{r}.
486       @param Msolve: solves Mx=r
487       @type Msolve: function C{Msolve(r)} where C{r} is of the same type like
488                     argument C{r}. The returned object needs to be of the same
489                     type like argument C{x}.
490       @param bilinearform: inner product C{<x,r>}
491       @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
492                           type like argument C{x} and C{r} is. The returned value
493                           is a C{float}.
494       @param atol: absolute tolerance
495       @type atol: non-negative C{float}
496       @param rtol: relative tolerance
497       @type rtol: non-negative C{float}
498       @param iter_max: maximum number of iteration steps
499       @type iter_max: C{int}
500       @return: the solution approximation and the corresponding residual
501       @rtype: C{tuple}
502       @warning: C{r} and C{x} are altered.
503       """
504       iter=0
505       rhat=Msolve(r)
506       d = rhat
507       rhat_dot_r = bilinearform(rhat, r)
508       if rhat_dot_r<0: raise NegativeNorm,"negative norm."
509       norm_r0=math.sqrt(rhat_dot_r)
510       atol2=atol+rtol*norm_r0
511       if atol2<=0:
512          raise ValueError,"Non-positive tolarance."
513       atol2=max(atol2, 100. * util.EPSILON * norm_r0)
514    
515       if verbose: print "PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_r0, atol2)
516    
517    
518       while not math.sqrt(rhat_dot_r) <= atol2:
519           iter+=1
520           if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
521    
522           q=Aprod(d)
523           alpha = rhat_dot_r / bilinearform(d, q)
524           x += alpha * d
525           r += (-alpha) * q
526    
527           rhat=Msolve(r)
528           rhat_dot_r_new = bilinearform(rhat, r)
529           beta = rhat_dot_r_new / rhat_dot_r
530           rhat+=beta * d
531           d=rhat
532    
533           rhat_dot_r = rhat_dot_r_new
534           if rhat_dot_r<0: raise NegativeNorm,"negative norm."
535           if verbose: print "PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))
536       if verbose: print "PCG: tolerance reached after %s steps."%iter
537       return x,r,math.sqrt(rhat_dot_r)
538    
539    class Defect(object):
540        """
541        Defines a non-linear defect F(x) of a variable x.
542        """
543        def __init__(self):
544            """
545            Initializes defect.
546            """
547            self.setDerivativeIncrementLength()
548    
549        def bilinearform(self, x0, x1):
550            """
551            Returns the inner product of x0 and x1
552    
553            @param x0: value for x0
554            @param x1: value for x1
555            @return: the inner product of x0 and x1
556            @rtype: C{float}
557            """
558            return 0
559    
560        def norm(self,x):
561            """
562            Returns the norm of argument C{x}.
563    
564            @param x: a value
565            @return: norm of argument x
566            @rtype: C{float}
567            @note: by default C{sqrt(self.bilinearform(x,x)} is returned.
568            """
569            s=self.bilinearform(x,x)
570            if s<0: raise NegativeNorm,"negative norm."
571            return math.sqrt(s)
572    
573        def eval(self,x):
574            """
575            Returns the value F of a given C{x}.
576    
577            @param x: value for which the defect C{F} is evaluated
578            @return: value of the defect at C{x}
579            """
580            return 0
581    
582        def __call__(self,x):
583            return self.eval(x)
584    
585        def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):
586            """
587            Sets the relative length of the increment used to approximate the
588            derivative of the defect. The increment is inc*norm(x)/norm(v)*v in the
589            direction of v with x as a starting point.
590    
591            @param inc: relative increment length
592            @type inc: positive C{float}
593            """
594            if inc<=0: raise ValueError,"positive increment required."
595            self.__inc=inc
596    
597        def getDerivativeIncrementLength(self):
598            """
599            Returns the relative increment length used to approximate the
600            derivative of the defect.
601            @return: value of the defect at C{x}
602            @rtype: positive C{float}
603            """
604            return self.__inc
605    
606        def derivative(self, F0, x0, v, v_is_normalised=True):
607            """
608            Returns the directional derivative at C{x0} in the direction of C{v}.
609    
610            @param F0: value of this defect at x0
611            @param x0: value at which derivative is calculated
612            @param v: direction
613            @param v_is_normalised: True to indicate that C{v} is nomalized
614                                    (self.norm(v)=0)
615            @return: derivative of this defect at x0 in the direction of C{v}
616            @note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is
617                   used but this method maybe overwritten to use exact evaluation.
618            """
619            normx=self.norm(x0)
620            if normx>0:
621                 epsnew = self.getDerivativeIncrementLength() * normx
622            else:
623                 epsnew = self.getDerivativeIncrementLength()
624            if not v_is_normalised:
625                normv=self.norm(v)
626                if normv<=0:
627                   return F0*0
628                else:
629                   epsnew /= normv
630            F1=self.eval(x0 + epsnew * v)
631            return (F1-F0)/epsnew
632    
633    ######################################
634    def NewtonGMRES(defect, x, iter_max=100, sub_iter_max=20, atol=0,rtol=1.e-4, sub_tol_max=0.5, gamma=0.9, verbose=False):
635       """
636       Solves a non-linear problem M{F(x)=0} for unknown M{x} using the stopping
637       criterion:
638    
639       M{norm(F(x) <= atol + rtol * norm(F(x0)}
640    
641       where M{x0} is the initial guess.
642    
643       @param defect: object defining the function M{F}. C{defect.norm} defines the
644                      M{norm} used in the stopping criterion.
645       @type defect: L{Defect}
646       @param x: initial guess for the solution, C{x} is altered.
647       @type x: any object type allowing basic operations such as
648                C{numarray.NumArray}, L{Data}
649       @param iter_max: maximum number of iteration steps
650       @type iter_max: positive C{int}
651       @param sub_iter_max: maximum number of inner iteration steps
652       @type sub_iter_max: positive C{int}
653       @param atol: absolute tolerance for the solution
654       @type atol: positive C{float}
655       @param rtol: relative tolerance for the solution
656       @type rtol: positive C{float}
657       @param gamma: tolerance safety factor for inner iteration
658       @type gamma: positive C{float}, less than 1
659       @param sub_tol_max: upper bound for inner tolerance
660       @type sub_tol_max: positive C{float}, less than 1
661       @return: an approximation of the solution with the desired accuracy
662       @rtype: same type as the initial guess
663       """
664       lmaxit=iter_max
665       if atol<0: raise ValueError,"atol needs to be non-negative."
666       if rtol<0: raise ValueError,"rtol needs to be non-negative."
667       if rtol+atol<=0: raise ValueError,"rtol or atol needs to be non-negative."
668       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
669       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
670    
671       F=defect(x)
672       fnrm=defect.norm(F)
673       stop_tol=atol + rtol*fnrm
674       sub_tol=sub_tol_max
675       if verbose: print "NewtonGMRES: initial residual = %e."%fnrm
676       if verbose: print "             tolerance = %e."%sub_tol
677       iter=1
678       #
679       # main iteration loop
680       #
681       while not fnrm<=stop_tol:
682                if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
683                #
684            #   adjust sub_tol_
685            #
686                if iter > 1:
687               rat=fnrm/fnrmo
688                   sub_tol_old=sub_tol
689               sub_tol=gamma*rat**2
690               if gamma*sub_tol_old**2 > .1: sub_tol=max(sub_tol,gamma*sub_tol_old**2)
691               sub_tol=max(min(sub_tol,sub_tol_max), .5*stop_tol/fnrm)
692            #
693            # calculate newton increment xc
694                #     if iter_max in __FDGMRES is reached MaxIterReached is thrown
695                #     if iter_restart -1 is returned as sub_iter
696                #     if  atol is reached sub_iter returns the numer of steps performed to get there
697                #
698                #
699                if verbose: print "             subiteration (GMRES) is called with relative tolerance %e."%sub_tol
700                try:
701                   xc, sub_iter=__FDGMRES(F, defect, x, sub_tol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)
702                except MaxIterReached:
703                   raise MaxIterReached,"maximum number of %s steps reached."%iter_max
704                if sub_iter<0:
705                   iter+=sub_iter_max
706                else:
707                   iter+=sub_iter
708                # ====
709            x+=xc
710                F=defect(x)
711            iter+=1
712                fnrmo, fnrm=fnrm, defect.norm(F)
713                if verbose: print "             step %s: residual %e."%(iter,fnrm)
714       if verbose: print "NewtonGMRES: completed after %s steps."%iter
715       return x
716    
717    def __givapp(c,s,vin):
718        """
719        Applies a sequence of Givens rotations (c,s) recursively to the vector
720        C{vin}
721    
722        @warning: C{vin} is altered.
723        """
724        vrot=vin
725        if isinstance(c,float):
726            vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
727        else:
728            for i in range(len(c)):
729                w1=c[i]*vrot[i]-s[i]*vrot[i+1]
730            w2=s[i]*vrot[i]+c[i]*vrot[i+1]
731                vrot[i]=w1
732                vrot[i+1]=w2
733        return vrot
734    
735    def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
736       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
737       c=numarray.zeros(iter_restart,numarray.Float64)
738       s=numarray.zeros(iter_restart,numarray.Float64)
739       g=numarray.zeros(iter_restart,numarray.Float64)
740       v=[]
741    
742       rho=defect.norm(F0)
743       if rho<=0.: return x0*0
744    
745       v.append(-F0/rho)
746       g[0]=rho
747       iter=0
748       while rho > atol and iter<iter_restart-1:
749            if iter  >= iter_max:
750                raise MaxIterReached,"maximum number of %s steps reached."%iter_max
751    
752            p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
753            v.append(p)
754    
755            v_norm1=defect.norm(v[iter+1])
756    
757            # Modified Gram-Schmidt
758            for j in range(iter+1):
759                h[j,iter]=defect.bilinearform(v[j],v[iter+1])
760                v[iter+1]-=h[j,iter]*v[j]
761    
762            h[iter+1,iter]=defect.norm(v[iter+1])
763            v_norm2=h[iter+1,iter]
764    
765            # Reorthogonalize if needed
766            if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
767                for j in range(iter+1):
768                    hr=defect.bilinearform(v[j],v[iter+1])
769                    h[j,iter]=h[j,iter]+hr
770                    v[iter+1] -= hr*v[j]
771    
772                v_norm2=defect.norm(v[iter+1])
773                h[iter+1,iter]=v_norm2
774            #   watch out for happy breakdown
775            if not v_norm2 == 0:
776                v[iter+1]=v[iter+1]/h[iter+1,iter]
777    
778            #   Form and store the information for the new Givens rotation
779            if iter > 0 :
780                hhat=numarray.zeros(iter+1,numarray.Float64)
781                for i in range(iter+1) : hhat[i]=h[i,iter]
782                hhat=__givapp(c[0:iter],s[0:iter],hhat);
783                for i in range(iter+1) : h[i,iter]=hhat[i]
784    
785            mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
786    
787            if mu!=0 :
788                c[iter]=h[iter,iter]/mu
789                s[iter]=-h[iter+1,iter]/mu
790                h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
791                h[iter+1,iter]=0.0
792                gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
793                g[iter]=gg[0]
794                g[iter+1]=gg[1]
795    
796            # Update the residual norm
797            rho=abs(g[iter+1])
798            iter+=1
799    
800       # At this point either iter > iter_max or rho < tol.
801       # It's time to compute x and leave.
802       if iter > 0 :
803         y=numarray.zeros(iter,numarray.Float64)
804         y[iter-1] = g[iter-1] / h[iter-1,iter-1]
805         if iter > 1 :
806            i=iter-2
807            while i>=0 :
808              y[i] = ( g[i] - numarray.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
809              i=i-1
810         xhat=v[iter-1]*y[iter-1]
811         for i in range(iter-1):
812        xhat += v[i]*y[i]
813       else :
814          xhat=v[0] * 0
815    
816       if iter<iter_restart-1:
817          stopped=iter
818       else:
819          stopped=-1
820    
821       return xhat,stopped
822    
823    def GMRES(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100, iter_restart=20, verbose=False,P_R=None):
824       """
825       Solver for
826    
827       M{Ax=b}
828    
829       with a general operator A (more details required!).
830       It uses the generalized minimum residual method (GMRES).
831    
832       The iteration is terminated if
833    
834       M{|r| <= atol+rtol*|r0|}
835    
836       where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
837    
838       M{|r| = sqrt( bilinearform(r,r))}
839    
840       @param r: initial residual M{r=b-Ax}. C{r} is altered.
841       @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
842       @param x: an initial guess for the solution
843       @type x: same like C{r}
844       @param Aprod: returns the value Ax
845       @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
846                    argument C{x}. The returned object needs to be of the same
847                    type like argument C{r}.
848       @param bilinearform: inner product C{<x,r>}
849       @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
850                           type like argument C{x} and C{r}. The returned value is
851                           a C{float}.
852       @param atol: absolute tolerance
853       @type atol: non-negative C{float}
854       @param rtol: relative tolerance
855       @type rtol: non-negative C{float}
856       @param iter_max: maximum number of iteration steps
857       @type iter_max: C{int}
858       @param iter_restart: in order to save memory the orthogonalization process
859                            is terminated after C{iter_restart} steps and the
860                            iteration is restarted.
861       @type iter_restart: C{int}
862       @return: the solution approximation and the corresponding residual
863       @rtype: C{tuple}
864       @warning: C{r} and C{x} are altered.
865       """
866       m=iter_restart
867       restarted=False
868       iter=0
869       if rtol>0:
870          r_dot_r = bilinearform(r, r)
871          if r_dot_r<0: raise NegativeNorm,"negative norm."
872          atol2=atol+rtol*math.sqrt(r_dot_r)
873          if verbose: print "GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)
874       else:
875          atol2=atol
876          if verbose: print "GMRES: absolute tolerance = %e"%atol2
877       if atol2<=0:
878          raise ValueError,"Non-positive tolarance."
879    
880       while True:
881          if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
882          if restarted:
883             r2 = r-Aprod(x-x2)
884          else:
885             r2=1*r
886          x2=x*1.
887          x,stopped=_GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose,P_R=P_R)
888          iter+=iter_restart
889          if stopped: break
890          if verbose: print "GMRES: restart."
891          restarted=True
892       if verbose: print "GMRES: tolerance has been reached."
893       return x
894    
895    def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False, P_R=None):
896       iter=0
897    
898       h=numarray.zeros((iter_restart+1,iter_restart),numarray.Float64)
899       c=numarray.zeros(iter_restart,numarray.Float64)
900       s=numarray.zeros(iter_restart,numarray.Float64)
901       g=numarray.zeros(iter_restart+1,numarray.Float64)
902       v=[]
903    
904       r_dot_r = bilinearform(r, r)
905       if r_dot_r<0: raise NegativeNorm,"negative norm."
906       rho=math.sqrt(r_dot_r)
907    
908       v.append(r/rho)
909       g[0]=rho
910    
911       if verbose: print "GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)
912       while not (rho<=atol or iter==iter_restart):
913    
914        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
915    
916            if P_R!=None:
917                p=Aprod(P_R(v[iter]))
918            else:
919            p=Aprod(v[iter])
920        v.append(p)
921    
922        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
923    
924    # Modified Gram-Schmidt
925        for j in range(iter+1):
926          h[j,iter]=bilinearform(v[j],v[iter+1])
927          v[iter+1]-=h[j,iter]*v[j]
928    
929        h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
930        v_norm2=h[iter+1,iter]
931    
932    # Reorthogonalize if needed
933        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
934         for j in range(iter+1):
935            hr=bilinearform(v[j],v[iter+1])
936                h[j,iter]=h[j,iter]+hr
937                v[iter+1] -= hr*v[j]
938    
939         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
940         h[iter+1,iter]=v_norm2
941    
942    #   watch out for happy breakdown
943            if not v_norm2 == 0:
944             v[iter+1]=v[iter+1]/h[iter+1,iter]
945    
946    #   Form and store the information for the new Givens rotation
947        if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])
948        mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
949    
950        if mu!=0 :
951            c[iter]=h[iter,iter]/mu
952            s[iter]=-h[iter+1,iter]/mu
953            h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
954            h[iter+1,iter]=0.0
955                    gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
956                    g[iter]=gg[0]
957                    g[iter+1]=gg[1]
958    # Update the residual norm
959    
960            rho=abs(g[iter+1])
961            if verbose: print "GMRES: iteration step %s: residual %e"%(iter,rho)
962        iter+=1
963    
964    # At this point either iter > iter_max or rho < tol.
965    # It's time to compute x and leave.
966    
967       if verbose: print "GMRES: iteration stopped after %s step."%iter
968       if iter > 0 :
969         y=numarray.zeros(iter,numarray.Float64)
970         y[iter-1] = g[iter-1] / h[iter-1,iter-1]
971         if iter > 1 :
972            i=iter-2
973            while i>=0 :
974              y[i] = ( g[i] - numarray.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
975              i=i-1
976         xhat=v[iter-1]*y[iter-1]
977         for i in range(iter-1):
978        xhat += v[i]*y[i]
979       else:
980         xhat=v[0] * 0
981       if P_R!=None:
982          x += P_R(xhat)
983       else:
984          x += xhat
985       if iter<iter_restart-1:
986          stopped=True
987       else:
988          stopped=False
989    
990       return x,stopped
991    
992    def MINRES(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
993        """
994        Solver for
995    
996        M{Ax=b}
997    
998        with a symmetric and positive definite operator A (more details required!).
999        It uses the minimum residual method (MINRES) with preconditioner M
1000        providing an approximation of A.
1001    
1002        The iteration is terminated if
1003    
1004        M{|r| <= atol+rtol*|r0|}
1005    
1006        where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
1007    
1008        M{|r| = sqrt( bilinearform(Msolve(r),r))}
1009    
1010        For details on the preconditioned conjugate gradient method see the book:
1011    
1012        I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
1013        T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
1014        C. Romine, and H. van der Vorst}.
1015    
1016        @param r: initial residual M{r=b-Ax}. C{r} is altered.
1017        @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1018        @param x: an initial guess for the solution
1019        @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1020        @param Aprod: returns the value Ax
1021        @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
1022                     argument C{x}. The returned object needs to be of the same
1023                     type like argument C{r}.
1024        @param Msolve: solves Mx=r
1025        @type Msolve: function C{Msolve(r)} where C{r} is of the same type like
1026                      argument C{r}. The returned object needs to be of the same
1027                      type like argument C{x}.
1028        @param bilinearform: inner product C{<x,r>}
1029        @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
1030                            type like argument C{x} and C{r} is. The returned value
1031                            is a C{float}.
1032        @param atol: absolute tolerance
1033        @type atol: non-negative C{float}
1034        @param rtol: relative tolerance
1035        @type rtol: non-negative C{float}
1036        @param iter_max: maximum number of iteration steps
1037        @type iter_max: C{int}
1038        @return: the solution approximation and the corresponding residual
1039        @rtype: C{tuple}
1040        @warning: C{r} and C{x} are altered.
1041        """
1042        #------------------------------------------------------------------
1043        # Set up y and v for the first Lanczos vector v1.
1044        # y  =  beta1 P' v1,  where  P = C**(-1).
1045        # v is really P' v1.
1046        #------------------------------------------------------------------
1047        r1    = r
1048        y = Msolve(r)
1049        beta1 = bilinearform(y,r)
1050    
1051        if beta1< 0: raise NegativeNorm,"negative norm."
1052    
1053        #  If r = 0 exactly, stop with x
1054        if beta1==0: return x
1055    
1056        if beta1> 0: beta1  = math.sqrt(beta1)
1057    
1058        #------------------------------------------------------------------
1059        # Initialize quantities.
1060        # ------------------------------------------------------------------
1061        iter   = 0
1062        Anorm = 0
1063        ynorm = 0
1064        oldb   = 0
1065        beta   = beta1
1066        dbar   = 0
1067        epsln  = 0
1068        phibar = beta1
1069        rhs1   = beta1
1070        rhs2   = 0
1071        rnorm  = phibar
1072        tnorm2 = 0
1073        ynorm2 = 0
1074        cs     = -1
1075        sn     = 0
1076        w      = r*0.
1077        w2     = r*0.
1078        r2     = r1
1079        eps    = 0.0001
1080    
1081        #---------------------------------------------------------------------
1082        # Main iteration loop.
1083        # --------------------------------------------------------------------
1084        while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL
1085    
1086        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1087            iter    = iter  +  1
1088    
1089            #-----------------------------------------------------------------
1090            # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
1091            # The general iteration is similar to the case k = 1 with v0 = 0:
1092            #
1093            #   p1      = Operator * v1  -  beta1 * v0,
1094            #   alpha1  = v1'p1,
1095            #   q2      = p2  -  alpha1 * v1,
1096            #   beta2^2 = q2'q2,
1097            #   v2      = (1/beta2) q2.
1098            #
1099            # Again, y = betak P vk,  where  P = C**(-1).
1100            #-----------------------------------------------------------------
1101            s = 1/beta                 # Normalize previous vector (in y).
1102            v = s*y                    # v = vk if P = I
1103    
1104            y      = Aprod(v)
1105    
1106            if iter >= 2:
1107              y = y - (beta/oldb)*r1
1108    
1109            alfa   = bilinearform(v,y)              # alphak
1110            y      += (- alfa/beta)*r2
1111            r1     = r2
1112            r2     = y
1113            y = Msolve(r2)
1114            oldb   = beta                           # oldb = betak
1115            beta   = bilinearform(y,r2)             # beta = betak+1^2
1116            if beta < 0: raise NegativeNorm,"negative norm."
1117    
1118            beta   = math.sqrt( beta )
1119            tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
1120    
1121            if iter==1:                 # Initialize a few things.
1122              gmax   = abs( alfa )      # alpha1
1123              gmin   = gmax             # alpha1
1124    
1125            # Apply previous rotation Qk-1 to get
1126            #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
1127            #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
1128    
1129            oldeps = epsln
1130            delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
1131            gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
1132            epsln  =               sn * beta  # epsln2 = 0         epslnk+1
1133            dbar   =            -  cs * beta  # dbar 2 = beta2     dbar k+1
1134    
1135            # Compute the next plane rotation Qk
1136    
1137            gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
1138            gamma  = max(gamma,eps)
1139            cs     = gbar / gamma             # ck
1140            sn     = beta / gamma             # sk
1141            phi    = cs * phibar              # phik
1142            phibar = sn * phibar              # phibark+1
1143    
1144            # Update  x.
1145    
1146            denom = 1/gamma
1147            w1    = w2
1148            w2    = w
1149            w     = (v - oldeps*w1 - delta*w2) * denom
1150            x     +=  phi*w
1151    
1152            # Go round again.
1153    
1154            gmax   = max(gmax,gamma)
1155            gmin   = min(gmin,gamma)
1156            z      = rhs1 / gamma
1157            ynorm2 = z*z  +  ynorm2
1158            rhs1   = rhs2 -  delta*z
1159            rhs2   =      -  epsln*z
1160    
1161            # Estimate various norms and test for convergence.
1162    
1163            Anorm  = math.sqrt( tnorm2 )
1164            ynorm  = math.sqrt( ynorm2 )
1165    
1166            rnorm  = phibar
1167    
1168        return x
1169    
1170    def TFQMR(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1171      """
1172      Solver for
1173    
1174      M{Ax=b}
1175    
1176      with a general operator A (more details required!).
1177      It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).
1178    
1179      The iteration is terminated if
1180    
1181      M{|r| <= atol+rtol*|r0|}
1182    
1183      where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
1184    
1185      M{|r| = sqrt( bilinearform(r,r))}
1186    
1187      @param r: initial residual M{r=b-Ax}. C{r} is altered.
1188      @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1189      @param x: an initial guess for the solution
1190      @type x: same like C{r}
1191      @param Aprod: returns the value Ax
1192      @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
1193                   argument C{x}. The returned object needs to be of the same type
1194                   like argument C{r}.
1195      @param bilinearform: inner product C{<x,r>}
1196      @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
1197                          type like argument C{x} and C{r}. The returned value is
1198                          a C{float}.
1199      @param atol: absolute tolerance
1200      @type atol: non-negative C{float}
1201      @param rtol: relative tolerance
1202      @type rtol: non-negative C{float}
1203      @param iter_max: maximum number of iteration steps
1204      @type iter_max: C{int}
1205      @rtype: C{tuple}
1206      @warning: C{r} and C{x} are altered.
1207      """
1208      u1=0
1209      u2=0
1210      y1=0
1211      y2=0
1212    
1213      w = r
1214      y1 = r
1215      iter = 0
1216      d = 0
1217      v = Aprod(y1)
1218      u1 = v
1219    
1220      theta = 0.0;
1221      eta = 0.0;
1222      rho=bilinearform(r,r)
1223      if rho < 0: raise NegativeNorm,"negative norm."
1224      tau = math.sqrt(rho)
1225      norm_r0=tau
1226      while tau>atol+rtol*norm_r0:
1227        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1228    
1229        sigma = bilinearform(r,v)
1230        if sigma == 0.0: raise IterationBreakDown,'TFQMR breakdown, sigma=0'
1231    
1232        alpha = rho / sigma
1233    
1234        for j in range(2):
1235    #
1236    #   Compute y2 and u2 only if you have to
1237    #
1238          if ( j == 1 ):
1239            y2 = y1 - alpha * v
1240            u2 = Aprod(y2)
1241    
1242          m = 2 * (iter+1) - 2 + (j+1)
1243          if j==0:
1244             w = w - alpha * u1
1245             d = y1 + ( theta * theta * eta / alpha ) * d
1246          if j==1:
1247             w = w - alpha * u2
1248             d = y2 + ( theta * theta * eta / alpha ) * d
1249    
1250          theta = math.sqrt(bilinearform(w,w))/ tau
1251          c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1252          tau = tau * theta * c
1253          eta = c * c * alpha
1254          x = x + eta * d
1255    #
1256    #  Try to terminate the iteration at each pass through the loop
1257    #
1258        if rho == 0.0: raise IterationBreakDown,'TFQMR breakdown, rho=0'
1259    
1260        rhon = bilinearform(r,w)
1261        beta = rhon / rho;
1262        rho = rhon;
1263        y1 = w + beta * y2;
1264        u1 = Aprod(y1)
1265        v = u1 + beta * ( u2 + beta * v )
1266    
1267        iter += 1
1268    
1269      return x
1270    
1271    
1272    #############################################
1273    
1274    class ArithmeticTuple(object):
1275       """
1276       Tuple supporting inplace update x+=y and scaling x=a*y where C{x,y} is an
1277       ArithmeticTuple and C{a} is a float.
1278    
1279       Example of usage::
1280    
1281           from esys.escript import Data
1282           from numarray import array
1283           a=Data(...)
1284           b=array([1.,4.])
1285           x=ArithmeticTuple(a,b)
1286           y=5.*x
1287    
1288       """
1289       def __init__(self,*args):
1290           """
1291           Initializes object with elements C{args}.
1292    
1293           @param args: tuple of objects that support inplace add (x+=y) and
1294                        scaling (x=a*y)
1295           """
1296           self.__items=list(args)
1297    
1298       def __len__(self):
1299           """
1300           Returns the number of items.
1301    
1302           @return: number of items
1303           @rtype: C{int}
1304           """
1305           return len(self.__items)
1306    
1307       def __getitem__(self,index):
1308           """
1309           Returns item at specified position.
1310    
1311           @param index: index of item to be returned
1312           @type index: C{int}
1313           @return: item with index C{index}
1314           """
1315           return self.__items.__getitem__(index)
1316    
1317       def __mul__(self,other):
1318           """
1319           Scales by C{other} from the right.
1320    
1321           @param other: scaling factor
1322           @type other: C{float}
1323           @return: itemwise self*other
1324           @rtype: L{ArithmeticTuple}
1325           """
1326           out=[]
1327           try:
1328               l=len(other)
1329               if l!=len(self):
1330                   raise ValueError,"length of arguments don't match."
1331               for i in range(l): out.append(self[i]*other[i])
1332           except TypeError:
1333               for i in range(len(self)): out.append(self[i]*other)
1334           return ArithmeticTuple(*tuple(out))
1335    
1336       def __rmul__(self,other):
1337           """
1338           Scales by C{other} from the left.
1339    
1340           @param other: scaling factor
1341           @type other: C{float}
1342           @return: itemwise other*self
1343           @rtype: L{ArithmeticTuple}
1344           """
1345           out=[]
1346           try:
1347               l=len(other)
1348               if l!=len(self):
1349                   raise ValueError,"length of arguments don't match."
1350               for i in range(l): out.append(other[i]*self[i])
1351           except TypeError:
1352               for i in range(len(self)): out.append(other*self[i])
1353           return ArithmeticTuple(*tuple(out))
1354    
1355       def __div__(self,other):
1356           """
1357           Scales by (1/C{other}) from the right.
1358    
1359           @param other: scaling factor
1360           @type other: C{float}
1361           @return: itemwise self/other
1362           @rtype: L{ArithmeticTuple}
1363           """
1364           return self*(1/other)
1365    
1366       def __rdiv__(self,other):
1367           """
1368           Scales by (1/C{other}) from the left.
1369    
1370           @param other: scaling factor
1371           @type other: C{float}
1372           @return: itemwise other/self
1373           @rtype: L{ArithmeticTuple}
1374           """
1375           out=[]
1376           try:
1377               l=len(other)
1378               if l!=len(self):
1379                   raise ValueError,"length of arguments don't match."
1380               for i in range(l): out.append(other[i]/self[i])
1381           except TypeError:
1382               for i in range(len(self)): out.append(other/self[i])
1383           return ArithmeticTuple(*tuple(out))
1384    
1385       def __iadd__(self,other):
1386           """
1387           Inplace addition of C{other} to self.
1388    
1389           @param other: increment
1390           @type other: C{ArithmeticTuple}
1391           """
1392           if len(self) != len(other):
1393               raise ValueError,"tuple lengths must match."
1394           for i in range(len(self)):
1395               self.__items[i]+=other[i]
1396           return self
1397    
1398       def __add__(self,other):
1399           """
1400           Adds C{other} to self.
1401    
1402           @param other: increment
1403           @type other: C{ArithmeticTuple}
1404           """
1405           out=[]
1406           try:
1407               l=len(other)
1408               if l!=len(self):
1409                   raise ValueError,"length of arguments don't match."
1410               for i in range(l): out.append(self[i]+other[i])
1411           except TypeError:
1412               for i in range(len(self)): out.append(self[i]+other)
1413           return ArithmeticTuple(*tuple(out))
1414    
1415       def __sub__(self,other):
1416           """
1417           Subtracts C{other} from self.
1418    
1419           @param other: decrement
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 __isub__(self,other):
1433           """
1434           Inplace subtraction of C{other} from self.
1435    
1436           @param other: decrement
1437           @type other: C{ArithmeticTuple}
1438           """
1439           if len(self) != len(other):
1440               raise ValueError,"tuple length must match."
1441           for i in range(len(self)):
1442               self.__items[i]-=other[i]
1443           return self
1444    
1445       def __neg__(self):
1446           """
1447           Negates values.
1448           """
1449           out=[]
1450           for i in range(len(self)):
1451               out.append(-self[i])
1452           return ArithmeticTuple(*tuple(out))
1453    
1454    
1455    class HomogeneousSaddlePointProblem(object):
1456          """
1457          This class provides a framework for solving linear homogeneous saddle
1458          point problems of the form::
1459    
1460              M{Av+B^*p=f}
1461              M{Bv     =0}
1462    
1463          for the unknowns M{v} and M{p} and given operators M{A} and M{B} and
1464          given right hand side M{f}. M{B^*} is the adjoint operator of M{B}.
1465          """
1466          def __init__(self,**kwargs):
1467            self.setTolerance()
1468            self.setAbsoluteTolerance()
1469            self.setSubProblemTolerance()
1470    
1471          #=============================================================
1472          def initialize(self):
1473            """
1474            Initializes the problem (overwrite).
1475            """
1476            pass
1477    
1478          def inner_pBv(self,p,v):
1479             """
1480             Returns inner product of element p and Bv (overwrite).
1481    
1482             @param p: a pressure increment
1483             @param v: a residual
1484             @return: inner product of element p and Bv
1485             @rtype: C{float}
1486             @note: used if PCG is applied.
1487             """
1488             raise NotImplementedError,"no inner product for p implemented."
1489    
1490          def inner_p(self,p0,p1):
1491             """
1492             Returns inner product of p0 and p1 (overwrite).
1493    
1494             @param p0: a pressure
1495             @param p1: a pressure
1496             @return: inner product of p0 and p1
1497             @rtype: C{float}
1498             """
1499             raise NotImplementedError,"no inner product for p implemented."
1500      
1501          def norm_v(self,v):
1502             """
1503             Returns the norm of v (overwrite).
1504    
1505             @param v: a velovity
1506             @return: norm of v
1507             @rtype: non-negative C{float}
1508             """
1509             raise NotImplementedError,"no norm of v implemented."
1510    
1511    
1512          def getV(self, p, v0):
1513             """
1514             return the value for v for a given p (overwrite)
1515    
1516             @param p: a pressure
1517             @param v0: a initial guess for the value v to return.
1518             @return: v given as M{v= A^{-1} (f-B^*p)}
1519             """
1520             raise NotImplementedError,"no v calculation implemented."
1521            
1522          def norm_Bv(self,v):
1523            """
1524            Returns Bv (overwrite).
1525    
1526            @rtype: equal to the type of p
1527            @note: boundary conditions on p should be zero!
1528            """
1529            raise NotImplementedError, "no operator B implemented."
1530    
1531          def solve_AinvBt(self,p):
1532             """
1533             Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
1534             (overwrite).
1535    
1536             @param p: a pressure increment
1537             @return: the solution of M{Av=B^*p}
1538             @note: boundary conditions on v should be zero!
1539             """
1540             raise NotImplementedError,"no operator A implemented."
1541    
1542          def solve_precB(self,v):
1543             """
1544             Provides a preconditioner for M{BA^{-1}B^*} with accuracy
1545             L{self.getSubProblemTolerance()} (overwrite).
1546    
1547             @rtype: equal to the type of p
1548             @note: boundary conditions on p should be zero!
1549             """
1550             raise NotImplementedError,"no preconditioner for Schur complement implemented."
1551          #=============================================================
1552          def __Aprod_PCG(self,p):
1553              return self.solve_AinvBt(p)
1554    
1555          def __inner_PCG(self,p,v):
1556             return self.inner_pBv(p,v)
1557    
1558          def __Msolve_PCG(self,v):
1559              return self.solve_precB(v)
1560          #=============================================================
1561          def __Aprod_GMRES(self,p):
1562              return self.solve_precB(self.solve_AinvBt(p))
1563          def __inner_GMRES(self,p0,p1):
1564             return self.inner_p(p0,p1)
1565          #=============================================================
1566          def norm_p(self,p):
1567              """
1568              calculates the norm of C{p}
1569              
1570              @param p: a pressure
1571              @return: the norm of C{p} using the inner product for pressure
1572              @rtype: C{float}
1573              """
1574              f=self.inner_p(p,p)
1575              if f<0: raise ValueError,"negative pressure norm."
1576              return math.sqrt(f)
1577              
1578    
1579          def solve(self,v,p,max_iter=20, verbose=False, show_details=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1580             """
1581             Solves the saddle point problem using initial guesses v and p.
1582    
1583             @param v: initial guess for velocity
1584             @param p: initial guess for pressure
1585             @type v: L{Data}
1586             @type p: L{Data}
1587             @param usePCG: indicates the usage of the PCG rather than GMRES scheme.
1588             @param max_iter: maximum number of iteration steps per correction
1589                              attempt
1590             @param verbose: if True, shows information on the progress of the
1591                             saddlepoint problem solver.
1592             @param show_details: if True, shows details of the sub problem solver
1593             @param iter_restart: restart the iteration after C{iter_restart} steps
1594                                  (only used if useUzaw=False)
1595             @type usePCG: C{bool}
1596             @type max_iter: C{int}
1597             @type verbose: C{bool}
1598             @type show_details: C{bool}
1599             @type iter_restart: C{int}
1600             @rtype: C{tuple} of L{Data} objects
1601             """
1602             self.verbose=verbose
1603             self.show_details=show_details and self.verbose
1604             rtol=self.getTolerance()
1605             atol=self.getAbsoluteTolerance()
1606             correction_step=0
1607             converged=False
1608             while not converged:
1609                  # calculate velocity for current pressure:
1610                  v=self.getV(p,v)
1611                  norm_v=self.norm_v(v)
1612                  norm_Bv=self.norm_Bv(v)
1613                  ATOL=norm_v*rtol+atol
1614                  if self.verbose: print "HomogeneousSaddlePointProblem: norm v= %e, norm_Bv= %e, tolerance = %e."%(norm_v, norm_Bv,ATOL)
1615                  if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."
1616                  if norm_Bv <= ATOL:
1617                     converged=True
1618                  else:
1619                     correction_step+=1
1620                     if correction_step>max_correction_steps:
1621                          raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1622                     dp=self.solve_precB(v)
1623                     if usePCG:
1624                       norm2=self.inner_pBv(dp,v)
1625                       if norm2<0: raise ValueError,"negative PCG norm."
1626                       norm2=math.sqrt(norm2)
1627                     else:
1628                       norm2=self.norm_p(dp)
1629                     ATOL_ITER=ATOL/norm_Bv*norm2*0.5
1630                     if self.verbose: print "HomogeneousSaddlePointProblem: tolerance for solver: %e"%ATOL_ITER
1631                     if usePCG:
1632                           p,v0,a_norm=PCG(v,self.__Aprod_PCG,p,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL_ITER, rtol=0.,iter_max=max_iter, verbose=self.verbose)
1633                     else:
1634                           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)
1635             if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached."
1636         return v,p
1637    
1638          #========================================================================
1639          def setTolerance(self,tolerance=1.e-4):
1640             """
1641             Sets the relative tolerance for (v,p).
1642    
1643             @param tolerance: tolerance to be used
1644             @type tolerance: non-negative C{float}
1645             """
1646             if tolerance<0:
1647                 raise ValueError,"tolerance must be positive."
1648             self.__rtol=tolerance
1649             self.setSubProblemTolerance()
1650    
1651          def getTolerance(self):
1652             """
1653             Returns the relative tolerance.
1654    
1655             @return: relative tolerance
1656             @rtype: C{float}
1657             """
1658             return self.__rtol
1659    
1660          def setAbsoluteTolerance(self,tolerance=0.):
1661             """
1662             Sets the absolute tolerance.
1663    
1664             @param tolerance: tolerance to be used
1665             @type tolerance: non-negative C{float}
1666             """
1667             if tolerance<0:
1668                 raise ValueError,"tolerance must be non-negative."
1669             self.__atol=tolerance
1670    
1671          def getAbsoluteTolerance(self):
1672             """
1673             Returns the absolute tolerance.
1674    
1675             @return: absolute tolerance
1676             @rtype: C{float}
1677             """
1678             return self.__atol
1679    
1680          def setSubProblemTolerance(self,rtol=None):
1681             """
1682             Sets the relative tolerance to solve the subproblem(s).
1683    
1684             @param rtol: relative tolerence
1685             @type rtol: positive C{float}
1686             """
1687             if rtol == None:
1688                  rtol=max(200.*util.EPSILON,self.getTolerance()**2)
1689             if rtol<=0:
1690                 raise ValueError,"tolerance must be positive."
1691             self.__sub_tol=rtol
1692    
1693          def getSubProblemTolerance(self):
1694             """
1695             Returns the subproblem reduction factor.
1696    
1697             @return: subproblem reduction factor
1698             @rtype: C{float}
1699             """
1700             return self.__sub_tol
1701    
1702    def MaskFromBoundaryTag(domain,*tags):
1703       """
1704       Creates a mask on the Solution(domain) function space where the value is
1705       one for samples that touch the boundary tagged by tags.
1706    
1707       Usage: m=MaskFromBoundaryTag(domain, "left", "right")
1708    
1709       @param domain: domain to be used
1710       @type domain: L{escript.Domain}
1711       @param tags: boundary tags
1712       @type tags: C{str}
1713       @return: a mask which marks samples that are touching the boundary tagged
1714                by any of the given tags
1715       @rtype: L{escript.Data} of rank 0
1716       """
1717       pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1718       d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
1719       for t in tags: d.setTaggedValue(t,1.)
1720       pde.setValue(y=d)
1721       return util.whereNonZero(pde.getRightHandSide())
1722    
1723    #==============================================================================
1724    class SaddlePointProblem(object):
1725       """
1726       This implements a solver for a saddle point problem
1727    
1728       M{f(u,p)=0}
1729       M{g(u)=0}
1730    
1731       for u and p. The problem is solved with an inexact Uszawa scheme for p:
1732    
1733       M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
1734       M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}
1735    
1736       where Q_f is an approximation of the Jacobian A_f of f with respect to u and
1737       Q_f is an approximation of A_g A_f^{-1} A_g with A_g is the Jacobian of g
1738       with respect to p. As a the construction of a 'proper' Q_g can be difficult,
1739       non-linear conjugate gradient method is applied to solve for p, so Q_g plays
1740       in fact the role of a preconditioner.
1741       """
1742       def __init__(self,verbose=False,*args):
1743           """
1744           Initializes the problem.
1745    
1746           @param verbose: if True, some information is printed in the process
1747           @type verbose: C{bool}
1748           @note: this method may be overwritten by a particular saddle point
1749                  problem
1750           """
1751           print "SaddlePointProblem should not be used anymore!"
1752           if not isinstance(verbose,bool):
1753                raise TypeError("verbose needs to be of type bool.")
1754           self.__verbose=verbose
1755           self.relaxation=1.
1756           DeprecationWarning("SaddlePointProblem should not be used anymore.")
1757    
1758       def trace(self,text):
1759           """
1760           Prints C{text} if verbose has been set.
1761    
1762           @param text: a text message
1763           @type text: C{str}
1764           """
1765           if self.__verbose: print "%s: %s"%(str(self),text)
1766    
1767       def solve_f(self,u,p,tol=1.e-8):
1768           """
1769           Solves
1770    
1771           A_f du = f(u,p)
1772    
1773           with tolerance C{tol} and returns du. A_f is Jacobian of f with respect
1774           to u.
1775    
1776           @param u: current approximation of u
1777           @type u: L{escript.Data}
1778           @param p: current approximation of p
1779           @type p: L{escript.Data}
1780           @param tol: tolerance expected for du
1781           @type tol: C{float}
1782           @return: increment du
1783           @rtype: L{escript.Data}
1784           @note: this method has to be overwritten by a particular saddle point
1785                  problem
1786           """
1787           pass
1788    
1789       def solve_g(self,u,tol=1.e-8):
1790           """
1791           Solves
1792    
1793           Q_g dp = g(u)
1794    
1795           where Q_g is a preconditioner for A_g A_f^{-1} A_g with A_g is the
1796           Jacobian of g with respect to p.
1797    
1798           @param u: current approximation of u
1799           @type u: L{escript.Data}
1800           @param tol: tolerance expected for dp
1801           @type tol: C{float}
1802           @return: increment dp
1803           @rtype: L{escript.Data}
1804           @note: this method has to be overwritten by a particular saddle point
1805                  problem
1806           """
1807           pass
1808    
1809       def inner(self,p0,p1):
1810           """
1811           Inner product of p0 and p1 approximating p. Typically this returns
1812           C{integrate(p0*p1)}.
1813           @return: inner product of p0 and p1
1814           @rtype: C{float}
1815           """
1816           pass
1817    
1818       subiter_max=3
1819       def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
1820            """
1821            Runs the solver.
1822    
1823            @param u0: initial guess for C{u}
1824            @type u0: L{esys.escript.Data}
1825            @param p0: initial guess for C{p}
1826            @type p0: L{esys.escript.Data}
1827            @param tolerance: tolerance for relative error in C{u} and C{p}
1828            @type tolerance: positive C{float}
1829            @param tolerance_u: tolerance for relative error in C{u} if different
1830                                from C{tolerance}
1831            @type tolerance_u: positive C{float}
1832            @param iter_max: maximum number of iteration steps
1833            @type iter_max: C{int}
1834            @param accepted_reduction: if the norm g cannot be reduced by
1835                                       C{accepted_reduction} backtracking to adapt
1836                                       the relaxation factor. If
1837                                       C{accepted_reduction=None} no backtracking
1838                                       is used.
1839            @type accepted_reduction: positive C{float} or C{None}
1840            @param relaxation: initial relaxation factor. If C{relaxation==None},
1841                               the last relaxation factor is used.
1842            @type relaxation: C{float} or C{None}
1843            """
1844            tol=1.e-2
1845            if tolerance_u==None: tolerance_u=tolerance
1846            if not relaxation==None: self.relaxation=relaxation
1847            if accepted_reduction ==None:
1848                  angle_limit=0.
1849            elif accepted_reduction>=1.:
1850                  angle_limit=0.
1851            else:
1852                  angle_limit=util.sqrt(1-accepted_reduction**2)
1853            self.iter=0
1854            u=u0
1855            p=p0
1856            #
1857            #   initialize things:
1858            #
1859            converged=False
1860            #
1861            #  start loop:
1862            #
1863            #  initial search direction is g
1864            #
1865            while not converged :
1866                if self.iter>iter_max:
1867                    raise ArithmeticError("no convergence after %s steps."%self.iter)
1868                f_new=self.solve_f(u,p,tol)
1869                norm_f_new = util.Lsup(f_new)
1870                u_new=u-f_new
1871                g_new=self.solve_g(u_new,tol)
1872                self.iter+=1
1873                norm_g_new = util.sqrt(self.inner(g_new,g_new))
1874                if norm_f_new==0. and norm_g_new==0.: return u, p
1875                if self.iter>1 and not accepted_reduction==None:
1876                   #
1877                   #   did we manage to reduce the norm of G? I
1878                   #   if not we start a backtracking procedure
1879                   #
1880                   # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
1881                   if norm_g_new > accepted_reduction * norm_g:
1882                      sub_iter=0
1883                      s=self.relaxation
1884                      d=g
1885                      g_last=g
1886                      self.trace("    start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
1887                      while sub_iter < self.subiter_max and  norm_g_new > accepted_reduction * norm_g:
1888                         dg= g_new-g_last
1889                         norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
1890                         rad=self.inner(g_new,dg)/self.relaxation
1891                         # print "   ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
1892                         # print "   ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
1893                         if abs(rad) < norm_dg*norm_g_new * angle_limit:
1894                             if sub_iter>0: self.trace("    no further improvements expected from backtracking.")
1895                             break
1896                         r=self.relaxation
1897                         self.relaxation= - rad/norm_dg**2
1898                         s+=self.relaxation
1899                         #####
1900                         # a=g_new+self.relaxation*dg/r
1901                         # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
1902                         #####
1903                         g_last=g_new
1904                         p+=self.relaxation*d
1905                         f_new=self.solve_f(u,p,tol)
1906                         u_new=u-f_new
1907                         g_new=self.solve_g(u_new,tol)
1908                         self.iter+=1
1909                         norm_f_new = util.Lsup(f_new)
1910                         norm_g_new = util.sqrt(self.inner(g_new,g_new))
1911                         # print "   ",sub_iter," new g norm",norm_g_new
1912                         self.trace("    %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
1913                         #
1914                         #   can we expect reduction of g?
1915                         #
1916                         # u_last=u_new
1917                         sub_iter+=1
1918                      self.relaxation=s
1919                #
1920                #  check for convergence:
1921                #
1922                norm_u_new = util.Lsup(u_new)
1923                p_new=p+self.relaxation*g_new
1924                norm_p_new = util.sqrt(self.inner(p_new,p_new))
1925                self.trace("%s th step: f = %s, f/u = %s, g = %s, g/p = %s, relaxation = %s."%(self.iter, norm_f_new ,norm_f_new/norm_u_new, norm_g_new, norm_g_new/norm_p_new, self.relaxation))
1926    
1927                if self.iter>1:
1928                   dg2=g_new-g
1929                   df2=f_new-f
1930                   norm_dg2=util.sqrt(self.inner(dg2,dg2))
1931                   norm_df2=util.Lsup(df2)
1932                   # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
1933                   tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
1934                   tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1935                   if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
1936                       converged=True
1937                f, norm_f, u, norm_u, g, norm_g, p, norm_p = f_new, norm_f_new, u_new, norm_u_new, g_new, norm_g_new, p_new, norm_p_new
1938            self.trace("convergence after %s steps."%self.iter)
1939            return u,p
1940    

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

  ViewVC Help
Powered by ViewVC 1.1.26