/[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 1956 by gross, Mon Nov 3 05:08:42 2008 UTC revision 3981 by jfenwick, Fri Sep 21 02:47:54 2012 UTC
# Line 1  Line 1 
1    
2  ########################################################  ##############################################################################
3  #  #
4  # Copyright (c) 2003-2008 by University of Queensland  # Copyright (c) 2003-2012 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # http://www.uq.edu.au
 # http://www.uq.edu.au/esscc  
6  #  #
7  # Primary Business: Queensland, Australia  # Primary Business: Queensland, Australia
8  # Licensed under the Open Software License version 3.0  # Licensed under the Open Software License version 3.0
9  # http://www.opensource.org/licenses/osl-3.0.php  # http://www.opensource.org/licenses/osl-3.0.php
10  #  #
11  ########################################################  # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    # Development since 2012 by School of Earth Sciences
13    #
14    ##############################################################################
15    
16  __copyright__="""Copyright (c) 2003-2008 by University of Queensland  __copyright__="""Copyright (c) 2003-2012 by University of Queensland
17  Earth Systems Science Computational Center (ESSCC)  http://www.uq.edu.au
 http://www.uq.edu.au/esscc  
18  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
19  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
20  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
21  __url__="http://www.uq.edu.au/esscc/escript-finley"  __url__="https://launchpad.net/escript-finley"
22    
23  """  """
24  Provides some tools related to PDEs.  Provides some tools related to PDEs.
25    
26  Currently includes:  Currently includes:
27      - Projector - to project a discontinuous      - Projector - to project a discontinuous function onto a continuous function
28      - Locator - to trace values in data objects at a certain location      - Locator - to trace values in data objects at a certain location
29      - TimeIntegrationManager - to handel extraplotion in time      - TimeIntegrationManager - to handle extrapolation in time
30      - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme      - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme
31    
32  @var __author__: name of author  :var __author__: name of author
33  @var __copyright__: copyrights  :var __copyright__: copyrights
34  @var __license__: licence agreement  :var __license__: licence agreement
35  @var __url__: url entry point on documentation  :var __url__: url entry point on documentation
36  @var __version__: version  :var __version__: version
37  @var __date__: date of the version  :var __date__: date of the version
38  """  """
39    
40  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
41    
42    
43  import escript  from . import escript
44  import linearPDEs  from . import linearPDEs
45  import numarray  import numpy
46  import util  from . import util
47  import math  import math
48    
 ##### Added by Artak  
 # from Numeric import zeros,Int,Float64  
 ###################################  
   
   
49  class TimeIntegrationManager:  class TimeIntegrationManager:
50    """    """
51    a simple mechanism to manage time dependend values.    A simple mechanism to manage time dependend values.
52    
53    typical usage is::    Typical usage is::
54    
55       dt=0.1 # time increment       dt=0.1 # time increment
56       tm=TimeIntegrationManager(inital_value,p=1)       tm=TimeIntegrationManager(inital_value,p=1)
# Line 64  class TimeIntegrationManager: Line 60  class TimeIntegrationManager:
60           tm.checkin(dt,v)           tm.checkin(dt,v)
61           t+=dt           t+=dt
62    
63    @note: currently only p=1 is supported.    :note: currently only p=1 is supported.
64    """    """
65    def __init__(self,*inital_values,**kwargs):    def __init__(self,*inital_values,**kwargs):
66       """       """
67       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 ``inital_values`` are the initial values
68         and p is the order used for extrapolation.
69       """       """
70       if kwargs.has_key("p"):       if "p" in kwargs:
71              self.__p=kwargs["p"]              self.__p=kwargs["p"]
72       else:       else:
73              self.__p=1              self.__p=1
74       if kwargs.has_key("time"):       if "time" in kwargs:
75              self.__t=kwargs["time"]              self.__t=kwargs["time"]
76       else:       else:
77              self.__t=0.              self.__t=0.
# Line 85  class TimeIntegrationManager: Line 82  class TimeIntegrationManager:
82    
83    def getTime(self):    def getTime(self):
84        return self.__t        return self.__t
85    
86    def getValue(self):    def getValue(self):
87        out=self.__v_mem[0]        out=self.__v_mem[0]
88        if len(out)==1:        if len(out)==1:
# Line 94  class TimeIntegrationManager: Line 92  class TimeIntegrationManager:
92    
93    def checkin(self,dt,*values):    def checkin(self,dt,*values):
94        """        """
95        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.
96        """        """
97        o=min(self.__order+1,self.__p)        o=min(self.__order+1,self.__p)
98        self.__order=min(self.__order+1,self.__p)        self.__order=min(self.__order+1,self.__p)
# Line 111  class TimeIntegrationManager: Line 109  class TimeIntegrationManager:
109    
110    def extrapolate(self,dt):    def extrapolate(self,dt):
111        """        """
112        extrapolates to dt forward in time.        Extrapolates to ``dt`` forward in time.
113        """        """
114        if self.__order==0:        if self.__order==0:
115           out=self.__v_mem[0]           out=self.__v_mem[0]
# Line 126  class TimeIntegrationManager: Line 124  class TimeIntegrationManager:
124           return out[0]           return out[0]
125        else:        else:
126           return out           return out
127    
128    
129  class Projector:  class Projector:
130    """    """
131    The Projector is a factory which projects a discontiuous function onto a    The Projector is a factory which projects a discontinuous function onto a
132    continuous function on the a given domain.    continuous function on a given domain.
133    """    """
134    def __init__(self, domain, reduce = True, fast=True):    def __init__(self, domain, reduce=True, fast=True):
135      """      """
136      Create a continuous function space projector for a domain.      Creates a continuous function space projector for a domain.
137    
138      @param domain: Domain of the projection.      :param domain: Domain of the projection.
139      @param reduce: Flag to reduce projection order (default is True)      :param reduce: Flag to reduce projection order
140      @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
141      """      """
142      self.__pde = linearPDEs.LinearPDE(domain)      self.__pde = linearPDEs.LinearPDE(domain)
143      if fast:      if fast:
144        self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)          self.__pde.getSolverOptions().setSolverMethod(linearPDEs.SolverOptions.LUMPING)
145      self.__pde.setSymmetryOn()      self.__pde.setSymmetryOn()
146      self.__pde.setReducedOrderTo(reduce)      self.__pde.setReducedOrderTo(reduce)
147      self.__pde.setValue(D = 1.)      self.__pde.setValue(D = 1.)
148      return      return
149      def getSolverOptions(self):
150        """
151        Returns the solver options of the PDE solver.
152        
153        :rtype: `linearPDEs.SolverOptions`
154        """
155        return self.__pde.getSolverOptions()
156    
157      def getValue(self, input_data):
158        """
159        Projects ``input_data`` onto a continuous function.
160    
161        :param input_data: the data to be projected
162        """
163        return self(input_data)
164    
165    def __call__(self, input_data):    def __call__(self, input_data):
166      """      """
167      Projects input_data onto a continuous function      Projects ``input_data`` onto a continuous function.
168    
169      @param input_data: The input_data to be projected.      :param input_data: the data to be projected
170      """      """
171      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
172      self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())      self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())
# Line 186  class Projector: Line 199  class Projector:
199    
200  class NoPDE:  class NoPDE:
201       """       """
202       solves the following problem for u:       Solves the following problem for u:
203    
204       M{kronecker[i,j]*D[j]*u[j]=Y[i]}       *kronecker[i,j]*D[j]*u[j]=Y[i]*
205    
206       with constraint       with constraint
207    
208       M{u[j]=r[j]}  where M{q[j]>0}       *u[j]=r[j]*  where *q[j]>0*
209    
210       where D, Y, r and q are given functions of rank 1.       where *D*, *Y*, *r* and *q* are given functions of rank 1.
211    
212       In the case of scalars this takes the form       In the case of scalars this takes the form
213    
214       M{D*u=Y}       *D*u=Y*
215    
216       with constraint       with constraint
217    
218       M{u=r}  where M{q>0}       *u=r* where *q>0*
219    
220       where D, Y, r and q are given scalar functions.       where *D*, *Y*, *r* and *q* are given scalar functions.
221    
222       The constraint is overwriting any other condition.       The constraint overwrites any other condition.
223    
224       @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 `linearPDEs.LinearPDE` class with
225              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
226              thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.              in `Solution` or `ReducedSolution`.
227       """       """
228         # The whole thing is a bit strange and I blame Rob Woodcock (CSIRO) for
229         # this.
230       def __init__(self,domain,D=None,Y=None,q=None,r=None):       def __init__(self,domain,D=None,Y=None,q=None,r=None):
231           """           """
232           initialize the problem           Initializes the problem.
233    
234           @param domain: domain of the PDE.           :param domain: domain of the PDE
235           @type domain: L{Domain}           :type domain: `Domain`
236           @param D: coefficient of the solution.           :param D: coefficient of the solution
237           @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type D: ``float``, ``int``, ``numpy.ndarray``, `Data`
238           @param Y: right hand side           :param Y: right hand side
239           @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type Y: ``float``, ``int``, ``numpy.ndarray``, `Data`
240           @param q: location of constraints           :param q: location of constraints
241           @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type q: ``float``, ``int``, ``numpy.ndarray``, `Data`
242           @param r: value of solution at locations of constraints           :param r: value of solution at locations of constraints
243           @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type r: ``float``, ``int``, ``numpy.ndarray``, `Data`
244           """           """
245           self.__domain=domain           self.__domain=domain
246           self.__D=D           self.__D=D
# Line 234  class NoPDE: Line 249  class NoPDE:
249           self.__r=r           self.__r=r
250           self.__u=None           self.__u=None
251           self.__function_space=escript.Solution(self.__domain)           self.__function_space=escript.Solution(self.__domain)
252    
253       def setReducedOn(self):       def setReducedOn(self):
254           """           """
255           sets the L{FunctionSpace} of the solution to L{ReducedSolution}           Sets the `FunctionSpace` of the solution to `ReducedSolution`.
256           """           """
257           self.__function_space=escript.ReducedSolution(self.__domain)           self.__function_space=escript.ReducedSolution(self.__domain)
258           self.__u=None           self.__u=None
259    
260       def setReducedOff(self):       def setReducedOff(self):
261           """           """
262           sets the L{FunctionSpace} of the solution to L{Solution}           Sets the `FunctionSpace` of the solution to `Solution`.
263           """           """
264           self.__function_space=escript.Solution(self.__domain)           self.__function_space=escript.Solution(self.__domain)
265           self.__u=None           self.__u=None
266            
267       def setValue(self,D=None,Y=None,q=None,r=None):       def setValue(self,D=None,Y=None,q=None,r=None):
268           """           """
269           assigns values to the parameters.           Assigns values to the parameters.
270    
271           @param D: coefficient of the solution.           :param D: coefficient of the solution
272           @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type D: ``float``, ``int``, ``numpy.ndarray``, `Data`
273           @param Y: right hand side           :param Y: right hand side
274           @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type Y: ``float``, ``int``, ``numpy.ndarray``, `Data`
275           @param q: location of constraints           :param q: location of constraints
276           @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type q: ``float``, ``int``, ``numpy.ndarray``, `Data`
277           @param r: value of solution at locations of constraints           :param r: value of solution at locations of constraints
278           @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}           :type r: ``float``, ``int``, ``numpy.ndarray``, `Data`
279           """           """
280           if not D==None:           if not D==None:
281              self.__D=D              self.__D=D
# Line 276  class NoPDE: Line 292  class NoPDE:
292    
293       def getSolution(self):       def getSolution(self):
294           """           """
295           returns the solution           Returns the solution.
296            
297           @return: the solution of the problem           :return: the solution of the problem
298           @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.           :rtype: `Data` object in the `FunctionSpace` `Solution` or
299                     `ReducedSolution`
300           """           """
301           if self.__u==None:           if self.__u==None:
302              if self.__D==None:              if self.__D==None:
303                 raise ValueError,"coefficient D is undefined"                 raise ValueError("coefficient D is undefined")
304              D=escript.Data(self.__D,self.__function_space)              D=escript.Data(self.__D,self.__function_space)
305              if D.getRank()>1:              if D.getRank()>1:
306                 raise ValueError,"coefficient D must have rank 0 or 1"                 raise ValueError("coefficient D must have rank 0 or 1")
307              if self.__Y==None:              if self.__Y==None:
308                 self.__u=escript.Data(0.,D.getShape(),self.__function_space)                 self.__u=escript.Data(0.,D.getShape(),self.__function_space)
309              else:              else:
310                 self.__u=util.quotient(self.__Y,D)                 self.__u=1./D*self.__Y
311              if not self.__q==None:              if not self.__q==None:
312                  q=util.wherePositive(escript.Data(self.__q,self.__function_space))                  q=util.wherePositive(escript.Data(self.__q,self.__function_space))
313                  self.__u*=(1.-q)                  self.__u*=(1.-q)
314                  if not self.__r==None: self.__u+=q*self.__r                  if not self.__r==None: self.__u+=q*self.__r
315           return self.__u           return self.__u
316                
317  class Locator:  class Locator:
318       """       """
319       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
320       spatial coordinate x.         coordinate x.
321        
322       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
323       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.  
324       """       """
325    
326       def __init__(self,where,x=numarray.zeros((3,))):       def __init__(self,where,x=numpy.zeros((3,))):
327         """         """
328         Initializes a Locator to access values in Data objects on the Doamin         Initializes a Locator to access values in Data objects on the Doamin
329         or FunctionSpace where for the sample point which         or FunctionSpace for the sample point which is closest to the given
330         closest to the given point x.         point x.
331    
332         @param where: function space         :param where: function space
333         @type where: L{escript.FunctionSpace}         :type where: `escript.FunctionSpace`
334         @param x: coefficient of the solution.         :param x: location(s) of the Locator
335         @type x: L{numarray.NumArray} or C{list} of L{numarray.NumArray}         :type x: ``numpy.ndarray`` or ``list`` of ``numpy.ndarray``
336         """         """
337         if isinstance(where,escript.FunctionSpace):         if isinstance(where,escript.FunctionSpace):
338            self.__function_space=where            self.__function_space=where
339         else:         else:
340            self.__function_space=escript.ContinuousFunction(where)            self.__function_space=escript.ContinuousFunction(where)
341           iterative=False
342         if isinstance(x, list):         if isinstance(x, list):
343               if len(x)==0:
344                  raise ValueError("At least one point must be given.")
345               try:
346                 iter(x[0])
347                 iterative=True
348               except TypeError:
349                 iterative=False
350           xxx=self.__function_space.getX()
351           if iterative:
352             self.__id=[]             self.__id=[]
353             for p in x:             for p in x:
354                self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())                self.__id.append(util.length(xxx-p[:self.__function_space.getDim()]).minGlobalDataPoint())
355         else:         else:
356             self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint()             self.__id=util.length(xxx-x[:self.__function_space.getDim()]).minGlobalDataPoint()
357    
358       def __str__(self):       def __str__(self):
359         """         """
360         Returns the coordinates of the Locator as a string.         Returns the coordinates of the Locator as a string.
361         """         """
362         x=self.getX()         x=self.getX()
363         if instance(x,list):         if isinstance(x,list):
364            out="["            out="["
365            first=True            first=True
366            for xx in x:            for xx in x:
# Line 350  class Locator: Line 376  class Locator:
376    
377       def getX(self):       def getX(self):
378          """          """
379      Returns the exact coordinates of the Locator.          Returns the exact coordinates of the Locator.
380      """          """
381          return self(self.getFunctionSpace().getX())          return self(self.getFunctionSpace().getX())
382    
383       def getFunctionSpace(self):       def getFunctionSpace(self):
384          """          """
385      Returns the function space of the Locator.          Returns the function space of the Locator.
386      """          """
387          return self.__function_space          return self.__function_space
388    
389       def getId(self,item=None):       def getId(self,item=None):
390          """          """
391      Returns the identifier of the location.          Returns the identifier of the location.
392      """          """
393          if item == None:          if item == None:
394             return self.__id             return self.__id
395          else:          else:
# Line 375  class Locator: Line 401  class Locator:
401    
402       def __call__(self,data):       def __call__(self,data):
403          """          """
404      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.
405      the object is returned.          """
     """  
406          return self.getValue(data)          return self.getValue(data)
407    
408       def getValue(self,data):       def getValue(self,data):
409          """          """
410      Returns the value of data at the Locator if data is a Data object          Returns the value of ``data`` at the Locator if ``data`` is a `Data`
411      otherwise the object is returned.          object otherwise the object is returned.
412      """          """
413          if isinstance(data,escript.Data):          if isinstance(data,escript.Data):
414             if data.getFunctionSpace()==self.getFunctionSpace():             dat=util.interpolate(data,self.getFunctionSpace())
              dat=data  
            else:  
              dat=data.interpolate(self.getFunctionSpace())  
415             id=self.getId()             id=self.getId()
416             r=data.getRank()             r=data.getRank()
417             if isinstance(id,list):             if isinstance(id,list):
418                 out=[]                 out=[]
419                 for i in id:                 for i in id:
420                    o=data.getValueOfGlobalDataPoint(*i)                    o=numpy.array(dat.getTupleForGlobalDataPoint(*i))
421                    if data.getRank()==0:                    if data.getRank()==0:
422                       out.append(o[0])                       out.append(o[0])
423                    else:                    else:
424                       out.append(o)                       out.append(o)
425                 return out                 return out
426             else:             else:
427               out=data.getValueOfGlobalDataPoint(*id)               out=numpy.array(dat.getTupleForGlobalDataPoint(*id))
428               if data.getRank()==0:               if data.getRank()==0:
429                  return out[0]                  return out[0]
430               else:               else:
431                  return out                  return out
432          else:          else:
433             return data             return data
434              
435         def setValue(self, data, v):
436          """
437          Sets the value of the ``data`` at the Locator.
438          """
439          if isinstance(data, escript.Data):
440             if data.getFunctionSpace()!=self.getFunctionSpace():
441               raise TypeError, "setValue: FunctionSpace of Locator and Data object must match."
442             data.expand()  
443             id=self.getId()
444             if isinstance(id, list):
445              for i in id:
446               data._setTupleForGlobalDataPoint(i[1], i[0], v)
447             else:
448               data._setTupleForGlobalDataPoint(id[1], id[0], v)
449          else:
450               raise TypeError, "setValue: Invalid argument type."
451    
452    
453    def getInfLocator(arg):
454        """
455        Return a Locator for a point with the inf value over all arg.
456        """
457        if not isinstance(arg, escript.Data):
458           raise TypeError("getInfLocator: Unknown argument type.")
459        a_inf=util.inf(arg)
460        loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords
461        x=arg.getFunctionSpace().getX()
462        x_min=x.getTupleForGlobalDataPoint(*loc)
463        return Locator(arg.getFunctionSpace(),x_min)
464    
465    def getSupLocator(arg):
466        """
467        Return a Locator for a point with the sup value over all arg.
468        """
469        if not isinstance(arg, escript.Data):
470           raise TypeError("getInfLocator: Unknown argument type.")
471        a_inf=util.sup(arg)
472        loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords
473        x=arg.getFunctionSpace().getX()
474        x_min=x.getTupleForGlobalDataPoint(*loc)
475        return Locator(arg.getFunctionSpace(),x_min)
476        
477    
478  class SolverSchemeException(Exception):  class SolverSchemeException(Exception):
479     """     """
480     exceptions thrown by solvers     This is a generic exception thrown by solvers.
481     """     """
482     pass     pass
483    
484  class IndefinitePreconditioner(SolverSchemeException):  class IndefinitePreconditioner(SolverSchemeException):
485     """     """
486     the preconditioner is not positive definite.     Exception thrown if the preconditioner is not positive definite.
487     """     """
488     pass     pass
489    
490  class MaxIterReached(SolverSchemeException):  class MaxIterReached(SolverSchemeException):
491     """     """
492     maxium number of iteration steps is reached.     Exception thrown if the maximum number of iteration steps is reached.
493     """     """
494     pass     pass
495  class IterationBreakDown(SolverSchemeException):  
496    class CorrectionFailed(SolverSchemeException):
497     """     """
498     iteration scheme econouters an incurable breakdown.     Exception thrown if no convergence has been achieved in the solution
499       correction scheme.
500     """     """
501     pass     pass
502  class NegativeNorm(SolverSchemeException):  
503    class IterationBreakDown(SolverSchemeException):
504     """     """
505     a norm calculation returns a negative norm.     Exception thrown if the iteration scheme encountered an incurable breakdown.
506     """     """
507     pass     pass
508    
509  class IterationHistory(object):  class NegativeNorm(SolverSchemeException):
510     """     """
511     The IterationHistory class is used to define a stopping criterium. It keeps track of the     Exception thrown if a norm calculation returns a negative norm.
    residual norms. The stoppingcriterium indicates termination if the residual norm has been reduced by  
    a given tolerance.  
512     """     """
513     def __init__(self,tolerance=math.sqrt(util.EPSILON),verbose=False):     pass
       """  
       Initialization  
   
       @param tolerance: tolerance  
       @type tolerance: positive C{float}  
       @param verbose: switches on the printing out some information  
       @type verbose: C{bool}  
       """  
       if not tolerance>0.:  
           raise ValueError,"tolerance needs to be positive."  
       self.tolerance=tolerance  
       self.verbose=verbose  
       self.history=[]  
    def stoppingcriterium(self,norm_r,r,x):  
        """  
        returns True if the C{norm_r} is C{tolerance}*C{norm_r[0]} where C{norm_r[0]}  is the residual norm at the first call.  
   
         
        @param norm_r: current residual norm  
        @type norm_r: non-negative C{float}  
        @param r: current residual (not used)  
        @param x: current solution approximation (not used)  
        @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.  
        @rtype: C{bool}  
   
        """  
        self.history.append(norm_r)  
        if self.verbose: print "iter: %s:  inner(rhat,r) = %e"#(len(self.history)-1, self.history[-1])  
        return self.history[-1]<=self.tolerance * self.history[0]  
514    
515     def stoppingcriterium2(self,norm_r,norm_b,solver="GMRES",TOL=None):  def PCG(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100, initial_guess=True, verbose=False):
516         """     """
517         returns True if the C{norm_r} is C{tolerance}*C{norm_b}     Solver for
518    
519             *Ax=b*
        @param norm_r: current residual norm  
        @type norm_r: non-negative C{float}  
        @param norm_b: norm of right hand side  
        @type norm_b: non-negative C{float}  
        @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.  
        @rtype: C{bool}  
520    
521         """     with a symmetric and positive definite operator A (more details required!).
522         if TOL==None:     It uses the conjugate gradient method with preconditioner M providing an
523            TOL=self.tolerance     approximation of A.
        self.history.append(norm_r)  
        if self.verbose: print "iter: %s:  norm(r) = %e"#(len(self.history)-1, self.history[-1])  
        return self.history[-1]<=TOL * norm_b  
524    
525  def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):     The iteration is terminated if
    """  
    Solver for  
526    
527     M{Ax=b}     *|r| <= atol+rtol*|r0|*
528    
529     with a symmetric and positive definite operator A (more details required!).     where *r0* is the initial residual and *|.|* is the energy norm. In fact
    It uses the conjugate gradient method with preconditioner M providing an approximation of A.  
530    
531     The iteration is terminated if the C{stoppingcriterium} function return C{True}.     *|r| = sqrt( bilinearform(Msolve(r),r))*
532    
533     For details on the preconditioned conjugate gradient method see the book:     For details on the preconditioned conjugate gradient method see the book:
534    
535     Templates for the Solution of Linear Systems by R. Barrett, M. Berry,     I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
536     T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,     T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
537     C. Romine, and H. van der Vorst.     C. Romine, and H. van der Vorst}.
538    
539     @param b: the right hand side of the liner system. C{b} is altered.     :param r: initial residual *r=b-Ax*. ``r`` is altered.
540     @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)     :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
541     @param Aprod: returns the value Ax     :param x: an initial guess for the solution
542     @type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}. The returned object needs to be of the same type like argument C{b}.     :type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
543     @param Msolve: solves Mx=r     :param Aprod: returns the value Ax
544     @type Msolve: function C{Msolve(r)} where C{r} is of the same type like argument C{b}. The returned object needs to be of the same     :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
545  type like argument C{x}.                  argument ``x``. The returned object needs to be of the same type
546     @param bilinearform: inner product C{<x,r>}                  like argument ``r``.
547     @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is . The returned value is a C{float}.     :param Msolve: solves Mx=r
548     @param stoppingcriterium: function which returns True if a stopping criterium is meet. C{stoppingcriterium} has the arguments C{norm_r}, C{r} and C{x} giving the current norm of the residual (=C{sqrt(bilinearform(Msolve(r),r)}), the current residual and the current solution approximation. C{stoppingcriterium} is called in each iteration step.     :type Msolve: function ``Msolve(r)`` where ``r`` is of the same type like
549     @type stoppingcriterium: function that returns C{True} or C{False}                   argument ``r``. The returned object needs to be of the same
550     @param x: an initial guess for the solution. If no C{x} is given 0*b is used.                   type like argument ``x``.
551     @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)     :param bilinearform: inner product ``<x,r>``
552     @param iter_max: maximum number of iteration steps.     :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
553     @type iter_max: C{int}                         type like argument ``x`` and ``r`` is. The returned value
554     @return: the solution approximation and the corresponding residual                         is a ``float``.
555     @rtype: C{tuple}     :param atol: absolute tolerance
556     @warning: C{b} and C{x} are altered.     :type atol: non-negative ``float``
557       :param rtol: relative tolerance
558       :type rtol: non-negative ``float``
559       :param iter_max: maximum number of iteration steps
560       :type iter_max: ``int``
561       :return: the solution approximation and the corresponding residual
562       :rtype: ``tuple``
563       :warning: ``r`` and ``x`` are altered.
564     """     """
565     iter=0     iter=0
    if x==None:  
       x=0*b  
    else:  
       b += (-1)*Aprod(x)  
    r=b  
566     rhat=Msolve(r)     rhat=Msolve(r)
567     d = rhat     d = rhat
568     rhat_dot_r = bilinearform(rhat, r)     rhat_dot_r = bilinearform(rhat, r)
569     if rhat_dot_r<0: raise NegativeNorm,"negative norm."     if rhat_dot_r<0: raise NegativeNorm("negative norm.")
570       norm_r0=math.sqrt(rhat_dot_r)
571       atol2=atol+rtol*norm_r0
572       if atol2<=0:
573          raise ValueError("Non-positive tolarance.")
574       atol2=max(atol2, 100. * util.EPSILON * norm_r0)
575    
576       if verbose: print(("PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_r0, atol2)))
577    
578    
579     while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):     while not math.sqrt(rhat_dot_r) <= atol2:
580         iter+=1         iter+=1
581         if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max         if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
582    
583         q=Aprod(d)         q=Aprod(d)
584         alpha = rhat_dot_r / bilinearform(d, q)         alpha = rhat_dot_r / bilinearform(d, q)
585         x += alpha * d         x += alpha * d
586         r += (-alpha) * q         if isinstance(q,ArithmeticTuple):
587              r += q * (-alpha)      # Doing it the other way calls the float64.__mul__ not AT.__rmul__
588           else:
589               r += (-alpha) * q
590         rhat=Msolve(r)         rhat=Msolve(r)
591         rhat_dot_r_new = bilinearform(rhat, r)         rhat_dot_r_new = bilinearform(rhat, r)
592         beta = rhat_dot_r_new / rhat_dot_r         beta = rhat_dot_r_new / rhat_dot_r
# Line 556  type like argument C{x}. Line 594  type like argument C{x}.
594         d=rhat         d=rhat
595    
596         rhat_dot_r = rhat_dot_r_new         rhat_dot_r = rhat_dot_r_new
597         if rhat_dot_r<0: raise NegativeNorm,"negative norm."         if rhat_dot_r<0: raise NegativeNorm("negative norm.")
598           if verbose: print(("PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))))
599     return x,r     if verbose: print(("PCG: tolerance reached after %s steps."%iter))
600       return x,r,math.sqrt(rhat_dot_r)
601    
602  class Defect(object):  class Defect(object):
603      """      """
604      defines a non-linear defect F(x) of a variable x      Defines a non-linear defect F(x) of a variable x.
605      """      """
606      def __init__(self):      def __init__(self):
607          """          """
608          initialize defect          Initializes defect.
609          """          """
610          self.setDerivativeIncrementLength()          self.setDerivativeIncrementLength()
611    
612      def bilinearform(self, x0, x1):      def bilinearform(self, x0, x1):
613          """          """
614          returns the inner product of x0 and x1          Returns the inner product of x0 and x1
615          @param x0: a value for x  
616          @param x1: a value for x          :param x0: value for x0
617          @return: the inner product of x0 and x1          :param x1: value for x1
618          @rtype: C{float}          :return: the inner product of x0 and x1
619            :rtype: ``float``
620          """          """
621          return 0          return 0
622          
623      def norm(self,x):      def norm(self,x):
624          """          """
625          the norm of argument C{x}          Returns the norm of argument ``x``.
626    
627          @param x: a value for x          :param x: a value
628          @return: norm of argument x          :return: norm of argument x
629          @rtype: C{float}          :rtype: ``float``
630          @note: by default C{sqrt(self.bilinearform(x,x)} is retrurned.          :note: by default ``sqrt(self.bilinearform(x,x)`` is returned.
631          """          """
632          s=self.bilinearform(x,x)          s=self.bilinearform(x,x)
633          if s<0: raise NegativeNorm,"negative norm."          if s<0: raise NegativeNorm("negative norm.")
634          return math.sqrt(s)          return math.sqrt(s)
635    
   
636      def eval(self,x):      def eval(self,x):
637          """          """
638          returns the value F of a given x          Returns the value F of a given ``x``.
639    
640          @param x: value for which the defect C{F} is evalulated.          :param x: value for which the defect ``F`` is evaluated
641          @return: value of the defect at C{x}          :return: value of the defect at ``x``
642          """          """
643          return 0          return 0
644    
645      def __call__(self,x):      def __call__(self,x):
646          return self.eval(x)          return self.eval(x)
647    
648      def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):      def setDerivativeIncrementLength(self,inc=1000.*math.sqrt(util.EPSILON)):
649          """          """
650          sets the relative length of the increment used to approximate the derivative of the defect          Sets the relative length of the increment used to approximate the
651          the increment is inc*norm(x)/norm(v)*v in the direction of v with x as a staring point.          derivative of the defect. The increment is inc*norm(x)/norm(v)*v in the
652            direction of v with x as a starting point.
653    
654          @param inc: relative increment length          :param inc: relative increment length
655          @type inc: positive C{float}          :type inc: positive ``float``
656          """          """
657          if inc<=0: raise ValueError,"positive increment required."          if inc<=0: raise ValueError("positive increment required.")
658          self.__inc=inc          self.__inc=inc
659    
660      def getDerivativeIncrementLength(self):      def getDerivativeIncrementLength(self):
661          """          """
662          returns the relative increment length used to approximate the derivative of the defect          Returns the relative increment length used to approximate the
663          @return: value of the defect at C{x}          derivative of the defect.
664          @rtype: positive C{float}          :return: value of the defect at ``x``
665            :rtype: positive ``float``
666          """          """
667          return self.__inc          return self.__inc
668    
669      def derivative(self, F0, x0, v, v_is_normalised=True):      def derivative(self, F0, x0, v, v_is_normalised=True):
670          """          """
671          returns the directional derivative at x0 in the direction of v          Returns the directional derivative at ``x0`` in the direction of ``v``.
672    
673          @param F0: value of this defect at x0          :param F0: value of this defect at x0
674          @param x0: value at which derivative is calculated.          :param x0: value at which derivative is calculated
675          @param v: direction          :param v: direction
676          @param v_is_normalised: is true to indicate that C{v} is nomalized (self.norm(v)=0)          :param v_is_normalised: True to indicate that ``v`` is nomalized
677          @return: derivative of this defect at x0 in the direction of C{v}                                  (self.norm(v)=0)
678          @note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is used but this method          :return: derivative of this defect at x0 in the direction of ``v``
679          maybe oepsnew verwritten to use exact evalution.          :note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is
680                   used but this method maybe overwritten to use exact evaluation.
681          """          """
682          normx=self.norm(x0)          normx=self.norm(x0)
683          if normx>0:          if normx>0:
# Line 651  class Defect(object): Line 693  class Defect(object):
693          F1=self.eval(x0 + epsnew * v)          F1=self.eval(x0 + epsnew * v)
694          return (F1-F0)/epsnew          return (F1-F0)/epsnew
695    
696  ######################################      ######################################
697  def NewtonGMRES(defect, x, iter_max=100, sub_iter_max=20, atol=0,rtol=1.e-4, sub_tol_max=0.5, gamma=0.9, verbose=False):  def NewtonGMRES(defect, x, iter_max=100, sub_iter_max=20, atol=0,rtol=1.e-4, subtol_max=0.5, gamma=0.9, verbose=False):
698     """     """
699     solves a non-linear problem M{F(x)=0} for unknown M{x} using the stopping criterion:     Solves a non-linear problem *F(x)=0* for unknown *x* using the stopping
700       criterion:
701    
702       *norm(F(x) <= atol + rtol * norm(F(x0)*
703    
704     M{norm(F(x) <= atol + rtol * norm(F(x0)}     where *x0* is the initial guess.
705      
706     where M{x0} is the initial guess.     :param defect: object defining the function *F*. ``defect.norm`` defines the
707                      *norm* used in the stopping criterion.
708     @param defect: object defining the the function M{F}, C{defect.norm} defines the M{norm} used in the stopping criterion.     :type defect: `Defect`
709     @type defect: L{Defect}     :param x: initial guess for the solution, ``x`` is altered.
710     @param x: initial guess for the solution, C{x} is altered.     :type x: any object type allowing basic operations such as
711     @type x: any object type allowing basic operations such as  L{numarray.NumArray}, L{Data}              ``numpy.ndarray``, `Data`
712     @param iter_max: maximum number of iteration steps     :param iter_max: maximum number of iteration steps
713     @type iter_max: positive C{int}     :type iter_max: positive ``int``
714     @param sub_iter_max:     :param sub_iter_max: maximum number of inner iteration steps
715     @type sub_iter_max:     :type sub_iter_max: positive ``int``
716     @param atol: absolute tolerance for the solution     :param atol: absolute tolerance for the solution
717     @type atol: positive C{float}     :type atol: positive ``float``
718     @param rtol: relative tolerance for the solution     :param rtol: relative tolerance for the solution
719     @type rtol: positive C{float}     :type rtol: positive ``float``
720     @param gamma: tolerance safety factor for inner iteration     :param gamma: tolerance safety factor for inner iteration
721     @type gamma: positive C{float}, less than 1     :type gamma: positive ``float``, less than 1
722     @param sub_tol_max: upper bound for inner tolerance.     :param subtol_max: upper bound for inner tolerance
723     @type sub_tol_max: positive C{float}, less than 1     :type subtol_max: positive ``float``, less than 1
724     @return: an approximation of the solution with the desired accuracy     :return: an approximation of the solution with the desired accuracy
725     @rtype: same type as the initial guess.     :rtype: same type as the initial guess
726     """     """
727     lmaxit=iter_max     lmaxit=iter_max
728     if atol<0: raise ValueError,"atol needs to be non-negative."     if atol<0: raise ValueError("atol needs to be non-negative.")
729     if rtol<0: raise ValueError,"rtol needs to be non-negative."     if rtol<0: raise ValueError("rtol needs to be non-negative.")
730     if rtol+atol<=0: raise ValueError,"rtol or atol needs to be non-negative."     if rtol+atol<=0: raise ValueError("rtol or atol needs to be non-negative.")
731     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     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)
732     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     if subtol_max<=0 or subtol_max>=1: raise ValueError("upper bound for inner tolerance for inner iteration (subtol_max =%s) needs to be positive and less than 1."%subtol_max)
733    
734     F=defect(x)     F=defect(x)
735     fnrm=defect.norm(F)     fnrm=defect.norm(F)
736     stop_tol=atol + rtol*fnrm     stop_tol=atol + rtol*fnrm
737     sub_tol=sub_tol_max     subtol=subtol_max
738     if verbose: print "NewtonGMRES: initial residual = %e."%fnrm     if verbose: print(("NewtonGMRES: initial residual = %e."%fnrm))
739     if verbose: print "             tolerance = %e."%sub_tol     if verbose: print(("             tolerance = %e."%subtol))
740     iter=1     iter=1
741     #     #
742     # main iteration loop     # main iteration loop
743     #     #
744     while not fnrm<=stop_tol:     while not fnrm<=stop_tol:
745              if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max              if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
746              #              #
747          #   adjust sub_tol_          #   adjust subtol_
748          #          #
749              if iter > 1:              if iter > 1:
750             rat=fnrm/fnrmo                 rat=fnrm/fnrmo
751                 sub_tol_old=sub_tol                 subtol_old=subtol
752             sub_tol=gamma*rat**2                 subtol=gamma*rat**2
753             if gamma*sub_tol_old**2 > .1: sub_tol=max(sub_tol,gamma*sub_tol_old**2)                 if gamma*subtol_old**2 > .1: subtol=max(subtol,gamma*subtol_old**2)
754             sub_tol=max(min(sub_tol,sub_tol_max), .5*stop_tol/fnrm)                 subtol=max(min(subtol,subtol_max), .5*stop_tol/fnrm)
755          #          #
756          # calculate newton increment xc          # calculate newton increment xc
757              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown
758              #     if iter_restart -1 is returned as sub_iter              #     if iter_restart -1 is returned as sub_iter
759              #     if  atol is reached sub_iter returns the numer of steps performed to get there              #     if  atol is reached sub_iter returns the numer of steps performed to get there
760              #              #
761              #                #
762              if verbose: print "             subiteration (GMRES) is called with relative tolerance %e."%sub_tol              if verbose: print(("             subiteration (GMRES) is called with relative tolerance %e."%subtol))
763              try:              try:
764                 xc, sub_iter=__FDGMRES(F, defect, x, sub_tol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)                 xc, sub_iter=__FDGMRES(F, defect, x, subtol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)
765              except MaxIterReached:              except MaxIterReached:
766                 raise MaxIterReached,"maximum number of %s steps reached."%iter_max                 raise MaxIterReached("maximum number of %s steps reached."%iter_max)
767              if sub_iter<0:              if sub_iter<0:
768                 iter+=sub_iter_max                 iter+=sub_iter_max
769              else:              else:
770                 iter+=sub_iter                 iter+=sub_iter
771              # ====              # ====
772          x+=xc              x+=xc
773              F=defect(x)              F=defect(x)
774          iter+=1              iter+=1
775              fnrmo, fnrm=fnrm, defect.norm(F)              fnrmo, fnrm=fnrm, defect.norm(F)
776              if verbose: print "             step %s: residual %e."%(iter,fnrm)              if verbose: print(("             step %s: residual %e."%(iter,fnrm)))
777     if verbose: print "NewtonGMRES: completed after %s steps."%iter     if verbose: print(("NewtonGMRES: completed after %s steps."%iter))
778     return x     return x
779    
780  def __givapp(c,s,vin):  def __givapp(c,s,vin):
781      """      """
782      apply a sequence of Givens rotations (c,s) to the recuirsively to the vector vin      Applies a sequence of Givens rotations (c,s) recursively to the vector
783      @warning: C{vin} is altered.      ``vin``
784    
785        :warning: ``vin`` is altered.
786      """      """
787      vrot=vin      vrot=vin
788      if isinstance(c,float):      if isinstance(c,float):
789          vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]          vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
790      else:      else:
791          for i in range(len(c)):          for i in range(len(c)):
792              w1=c[i]*vrot[i]-s[i]*vrot[i+1]              w1=c[i]*vrot[i]-s[i]*vrot[i+1]
793          w2=s[i]*vrot[i]+c[i]*vrot[i+1]              w2=s[i]*vrot[i]+c[i]*vrot[i+1]
794              vrot[i:i+2]=w1,w2              vrot[i]=w1
795                vrot[i+1]=w2
796      return vrot      return vrot
797    
798  def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):  def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
799     h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)     h=numpy.zeros((iter_restart,iter_restart),numpy.float64)
800     c=numarray.zeros(iter_restart,numarray.Float64)     c=numpy.zeros(iter_restart,numpy.float64)
801     s=numarray.zeros(iter_restart,numarray.Float64)     s=numpy.zeros(iter_restart,numpy.float64)
802     g=numarray.zeros(iter_restart,numarray.Float64)     g=numpy.zeros(iter_restart,numpy.float64)
803     v=[]     v=[]
804    
805     rho=defect.norm(F0)     rho=defect.norm(F0)
806     if rho<=0.: return x0*0     if rho<=0.: return x0*0
807      
808     v.append(-F0/rho)     v.append(-F0/rho)
809     g[0]=rho     g[0]=rho
810     iter=0     iter=0
811     while rho > atol and iter<iter_restart-1:     while rho > atol and iter<iter_restart-1:
812            if iter  >= iter_max:
813      if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max              raise MaxIterReached("maximum number of %s steps reached."%iter_max)
814    
815          p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)          p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
816      v.append(p)          v.append(p)
817    
818            v_norm1=defect.norm(v[iter+1])
819    
820      v_norm1=defect.norm(v[iter+1])          # Modified Gram-Schmidt
821            for j in range(iter+1):
822                h[j,iter]=defect.bilinearform(v[j],v[iter+1])
823                v[iter+1]-=h[j,iter]*v[j]
824    
825          # Modified Gram-Schmidt          h[iter+1,iter]=defect.norm(v[iter+1])
826      for j in range(iter+1):          v_norm2=h[iter+1,iter]
          h[j][iter]=defect.bilinearform(v[j],v[iter+1])    
          v[iter+1]-=h[j][iter]*v[j]  
         
     h[iter+1][iter]=defect.norm(v[iter+1])  
     v_norm2=h[iter+1][iter]  
827    
828          # Reorthogonalize if needed          # Reorthogonalize if needed
829      if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)          if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
830          for j in range(iter+1):                for j in range(iter+1):
831             hr=defect.bilinearform(v[j],v[iter+1])                  hr=defect.bilinearform(v[j],v[iter+1])
832                 h[j][iter]=h[j][iter]+hr                  h[j,iter]=h[j,iter]+hr
833                 v[iter+1] -= hr*v[j]                  v[iter+1] -= hr*v[j]
834    
835          v_norm2=defect.norm(v[iter+1])              v_norm2=defect.norm(v[iter+1])
836          h[iter+1][iter]=v_norm2              h[iter+1,iter]=v_norm2
837          #   watch out for happy breakdown          #   watch out for happy breakdown
838          if not v_norm2 == 0:          if not v_norm2 == 0:
839                  v[iter+1]=v[iter+1]/h[iter+1][iter]              v[iter+1]=v[iter+1]/h[iter+1,iter]
840    
841          #   Form and store the information for the new Givens rotation          #   Form and store the information for the new Givens rotation
842      if iter > 0 :          if iter > 0 :
843          hhat=numarray.zeros(iter+1,numarray.Float64)              hhat=numpy.zeros(iter+1,numpy.float64)
844          for i in range(iter+1) : hhat[i]=h[i][iter]              for i in range(iter+1) : hhat[i]=h[i,iter]
845          hhat=__givapp(c[0:iter],s[0:iter],hhat);              hhat=__givapp(c[0:iter],s[0:iter],hhat);
846              for i in range(iter+1) : h[i][iter]=hhat[i]              for i in range(iter+1) : h[i,iter]=hhat[i]
847    
848      mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])          mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
849    
850      if mu!=0 :          if mu!=0 :
851          c[iter]=h[iter][iter]/mu              c[iter]=h[iter,iter]/mu
852          s[iter]=-h[iter+1][iter]/mu              s[iter]=-h[iter+1,iter]/mu
853          h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]              h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
854          h[iter+1][iter]=0.0              h[iter+1,iter]=0.0
855          g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])              gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
856                g[iter]=gg[0]
857                g[iter+1]=gg[1]
858    
859          # Update the residual norm          # Update the residual norm
860          rho=abs(g[iter+1])          rho=abs(g[iter+1])
861      iter+=1          iter+=1
862    
863     # At this point either iter > iter_max or rho < tol.     # At this point either iter > iter_max or rho < tol.
864     # It's time to compute x and leave.             # It's time to compute x and leave.
865     if iter > 0 :     if iter > 0 :
866       y=numarray.zeros(iter,numarray.Float64)           y=numpy.zeros(iter,numpy.float64)
867       y[iter-1] = g[iter-1] / h[iter-1][iter-1]       y[iter-1] = g[iter-1] / h[iter-1,iter-1]
868       if iter > 1 :         if iter > 1 :
869          i=iter-2            i=iter-2
870          while i>=0 :          while i>=0 :
871            y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]            y[i] = ( g[i] - numpy.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
872            i=i-1            i=i-1
873       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
874       for i in range(iter-1):       for i in range(iter-1):
875      xhat += v[i]*y[i]         xhat += v[i]*y[i]
876     else :     else :
877        xhat=v[0] * 0        xhat=v[0] * 0
878    
879     if iter<iter_restart-1:     if iter<iter_restart-1:
880        stopped=iter        stopped=iter
881     else:     else:
882        stopped=-1        stopped=-1
883    
884     return xhat,stopped     return xhat,stopped
885    
886  ##############################################  def GMRES(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100, iter_restart=20, verbose=False,P_R=None):
887  def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):     """
888  ################################################     Solver for
889    
890       *Ax=b*
891    
892       with a general operator A (more details required!).
893       It uses the generalized minimum residual method (GMRES).
894    
895       The iteration is terminated if
896    
897       *|r| <= atol+rtol*|r0|*
898    
899       where *r0* is the initial residual and *|.|* is the energy norm. In fact
900    
901       *|r| = sqrt( bilinearform(r,r))*
902    
903       :param r: initial residual *r=b-Ax*. ``r`` is altered.
904       :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
905       :param x: an initial guess for the solution
906       :type x: same like ``r``
907       :param Aprod: returns the value Ax
908       :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
909                    argument ``x``. The returned object needs to be of the same
910                    type like argument ``r``.
911       :param bilinearform: inner product ``<x,r>``
912       :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
913                           type like argument ``x`` and ``r``. The returned value is
914                           a ``float``.
915       :param atol: absolute tolerance
916       :type atol: non-negative ``float``
917       :param rtol: relative tolerance
918       :type rtol: non-negative ``float``
919       :param iter_max: maximum number of iteration steps
920       :type iter_max: ``int``
921       :param iter_restart: in order to save memory the orthogonalization process
922                            is terminated after ``iter_restart`` steps and the
923                            iteration is restarted.
924       :type iter_restart: ``int``
925       :return: the solution approximation and the corresponding residual
926       :rtype: ``tuple``
927       :warning: ``r`` and ``x`` are altered.
928       """
929     m=iter_restart     m=iter_restart
930       restarted=False
931     iter=0     iter=0
932     xc=x     if rtol>0:
933          r_dot_r = bilinearform(r, r)
934          if r_dot_r<0: raise NegativeNorm("negative norm.")
935          atol2=atol+rtol*math.sqrt(r_dot_r)
936          if verbose: print(("GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)))
937       else:
938          atol2=atol
939          if verbose: print(("GMRES: absolute tolerance = %e"%atol2))
940       if atol2<=0:
941          raise ValueError("Non-positive tolarance.")
942    
943     while True:     while True:
944        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max        if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached"%iter_max)
945        xc,stopped=__GMRESm(b*1, Aprod, Msolve, bilinearform, stoppingcriterium, x=xc*1, iter_max=iter_max-iter, iter_restart=m)        if restarted:
946             r2 = r-Aprod(x-x2)
947          else:
948             r2=1*r
949          x2=x*1.
950          x,stopped=_GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose,P_R=P_R)
951          iter+=iter_restart
952        if stopped: break        if stopped: break
953        iter+=iter_restart            if verbose: print("GMRES: restart.")
954     return xc        restarted=True
955       if verbose: print("GMRES: tolerance has been reached.")
956       return x
957    
958  def __GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):  def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False, P_R=None):
959     iter=0     iter=0
    r=Msolve(b)  
    r_dot_r = bilinearform(r, r)  
    if r_dot_r<0: raise NegativeNorm,"negative norm."  
    norm_b=math.sqrt(r_dot_r)  
960    
961     if x==None:     h=numpy.zeros((iter_restart+1,iter_restart),numpy.float64)
962        x=0*b     c=numpy.zeros(iter_restart,numpy.float64)
963     else:     s=numpy.zeros(iter_restart,numpy.float64)
964        r=Msolve(b-Aprod(x))     g=numpy.zeros(iter_restart+1,numpy.float64)
       r_dot_r = bilinearform(r, r)  
       if r_dot_r<0: raise NegativeNorm,"negative norm."  
     
    h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)  
    c=numarray.zeros(iter_restart,numarray.Float64)  
    s=numarray.zeros(iter_restart,numarray.Float64)  
    g=numarray.zeros(iter_restart,numarray.Float64)  
965     v=[]     v=[]
966    
967       r_dot_r = bilinearform(r, r)
968       if r_dot_r<0: raise NegativeNorm("negative norm.")
969     rho=math.sqrt(r_dot_r)     rho=math.sqrt(r_dot_r)
970      
971     v.append(r/rho)     v.append(r/rho)
972     g[0]=rho     g[0]=rho
973    
974     while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):     if verbose: print(("GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)))
975       while not (rho<=atol or iter==iter_restart):
976    
977      if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max          if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
978    
979      p=Msolve(Aprod(v[iter]))          if P_R!=None:
980                p=Aprod(P_R(v[iter]))
981            else:
982                p=Aprod(v[iter])
983            v.append(p)
984    
985      v.append(p)          v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
986    
987      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))    # Modified Gram-Schmidt
988            for j in range(iter+1):
989              h[j,iter]=bilinearform(v[j],v[iter+1])
990              v[iter+1]-=h[j,iter]*v[j]
991    
992  # Modified Gram-Schmidt          h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
993      for j in range(iter+1):          v_norm2=h[iter+1,iter]
       h[j][iter]=bilinearform(v[j],v[iter+1])    
       v[iter+1]-=h[j][iter]*v[j]  
         
     h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))  
     v_norm2=h[iter+1][iter]  
994    
995  # Reorthogonalize if needed  # Reorthogonalize if needed
996      if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)          if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
997       for j in range(iter+1):             for j in range(iter+1):
998          hr=bilinearform(v[j],v[iter+1])              hr=bilinearform(v[j],v[iter+1])
999              h[j][iter]=h[j][iter]+hr              h[j,iter]=h[j,iter]+hr
1000              v[iter+1] -= hr*v[j]              v[iter+1] -= hr*v[j]
1001    
1002       v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))             v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
1003       h[iter+1][iter]=v_norm2           h[iter+1,iter]=v_norm2
1004    
1005  #   watch out for happy breakdown  #   watch out for happy breakdown
1006          if not v_norm2 == 0:          if not v_norm2 == 0:
1007           v[iter+1]=v[iter+1]/h[iter+1][iter]           v[iter+1]=v[iter+1]/h[iter+1,iter]
1008    
1009  #   Form and store the information for the new Givens rotation  #   Form and store the information for the new Givens rotation
1010      if iter > 0 :          if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])
1011          hhat=numarray.zeros(iter+1,numarray.Float64)          mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
         for i in range(iter+1) : hhat[i]=h[i][iter]  
         hhat=__givapp(c[0:iter],s[0:iter],hhat);  
             for i in range(iter+1) : h[i][iter]=hhat[i]  
   
     mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])  
   
     if mu!=0 :  
         c[iter]=h[iter][iter]/mu  
         s[iter]=-h[iter+1][iter]/mu  
         h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]  
         h[iter+1][iter]=0.0  
         g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])  
1012    
1013            if mu!=0 :
1014                    c[iter]=h[iter,iter]/mu
1015                    s[iter]=-h[iter+1,iter]/mu
1016                    h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
1017                    h[iter+1,iter]=0.0
1018                    gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
1019                    g[iter]=gg[0]
1020                    g[iter+1]=gg[1]
1021  # Update the residual norm  # Update the residual norm
1022                  
1023          rho=abs(g[iter+1])          rho=abs(g[iter+1])
1024      iter+=1          if verbose: print(("GMRES: iteration step %s: residual %e"%(iter,rho)))
1025            iter+=1
1026    
1027  # At this point either iter > iter_max or rho < tol.  # At this point either iter > iter_max or rho < tol.
1028  # It's time to compute x and leave.          # It's time to compute x and leave.
1029    
1030     if iter > 0 :     if verbose: print(("GMRES: iteration stopped after %s step."%iter))
1031       y=numarray.zeros(iter,numarray.Float64)         if iter > 0 :
1032       y[iter-1] = g[iter-1] / h[iter-1][iter-1]       y=numpy.zeros(iter,numpy.float64)
1033       if iter > 1 :         y[iter-1] = g[iter-1] / h[iter-1,iter-1]
1034          i=iter-2         if iter > 1 :
1035            i=iter-2
1036          while i>=0 :          while i>=0 :
1037            y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]            y[i] = ( g[i] - numpy.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
1038            i=i-1            i=i-1
1039       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
1040       for i in range(iter-1):       for i in range(iter-1):
1041      xhat += v[i]*y[i]         xhat += v[i]*y[i]
1042     else : xhat=v[0]     else:
1043         xhat=v[0] * 0
1044     x += xhat     if P_R!=None:
1045     if iter<iter_restart-1:        x += P_R(xhat)
1046        stopped=True     else:
1047     else:        x += xhat
1048       if iter<iter_restart-1:
1049          stopped=True
1050       else:
1051        stopped=False        stopped=False
1052    
1053     return x,stopped     return x,stopped
1054    
1055  #################################################  def MINRES(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1056  def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):      """
1057  #################################################      Solver for
1058      #  
1059      #  minres solves the system of linear equations Ax = b      *Ax=b*
1060      #  where A is a symmetric matrix (possibly indefinite or singular)  
1061      #  and b is a given vector.      with a symmetric and positive definite operator A (more details required!).
1062      #        It uses the minimum residual method (MINRES) with preconditioner M
1063      #  "A" may be a dense or sparse matrix (preferably sparse!)      providing an approximation of A.
1064      #  or the name of a function such that  
1065      #               y = A(x)      The iteration is terminated if
1066      #  returns the product y = Ax for any given vector x.  
1067      #      *|r| <= atol+rtol*|r0|*
1068      #  "M" defines a positive-definite preconditioner M = C C'.  
1069      #  "M" may be a dense or sparse matrix (preferably sparse!)      where *r0* is the initial residual and *|.|* is the energy norm. In fact
1070      #  or the name of a function such that  
1071      #  solves the system My = x for any given vector x.      *|r| = sqrt( bilinearform(Msolve(r),r))*
1072      #  
1073      #      For details on the preconditioned conjugate gradient method see the book:
1074        
1075        I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
1076        T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
1077        C. Romine, and H. van der Vorst}.
1078    
1079        :param r: initial residual *r=b-Ax*. ``r`` is altered.
1080        :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1081        :param x: an initial guess for the solution
1082        :type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1083        :param Aprod: returns the value Ax
1084        :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
1085                     argument ``x``. The returned object needs to be of the same
1086                     type like argument ``r``.
1087        :param Msolve: solves Mx=r
1088        :type Msolve: function ``Msolve(r)`` where ``r`` is of the same type like
1089                      argument ``r``. The returned object needs to be of the same
1090                      type like argument ``x``.
1091        :param bilinearform: inner product ``<x,r>``
1092        :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
1093                            type like argument ``x`` and ``r`` is. The returned value
1094                            is a ``float``.
1095        :param atol: absolute tolerance
1096        :type atol: non-negative ``float``
1097        :param rtol: relative tolerance
1098        :type rtol: non-negative ``float``
1099        :param iter_max: maximum number of iteration steps
1100        :type iter_max: ``int``
1101        :return: the solution approximation and the corresponding residual
1102        :rtype: ``tuple``
1103        :warning: ``r`` and ``x`` are altered.
1104        """
1105      #------------------------------------------------------------------      #------------------------------------------------------------------
1106      # Set up y and v for the first Lanczos vector v1.      # Set up y and v for the first Lanczos vector v1.
1107      # y  =  beta1 P' v1,  where  P = C**(-1).      # y  =  beta1 P' v1,  where  P = C**(-1).
1108      # v is really P' v1.      # v is really P' v1.
1109      #------------------------------------------------------------------      #------------------------------------------------------------------
1110      if x==None:      r1    = r
1111        x=0*b      y = Msolve(r)
1112      else:      beta1 = bilinearform(y,r)
       b += (-1)*Aprod(x)  
1113    
1114      r1    = b      if beta1< 0: raise NegativeNorm("negative norm.")
     y = Msolve(b)  
     beta1 = bilinearform(y,b)  
   
     if beta1< 0: raise NegativeNorm,"negative norm."  
1115    
1116      #  If b = 0 exactly, stop with x = 0.      #  If r = 0 exactly, stop with x
1117      if beta1==0: return x*0.      if beta1==0: return x
1118    
1119      if beta1> 0:      if beta1> 0: beta1  = math.sqrt(beta1)
       beta1  = math.sqrt(beta1)        
1120    
1121      #------------------------------------------------------------------      #------------------------------------------------------------------
1122      # Initialize quantities.      # Initialize quantities.
# Line 1008  def MINRES(b, Aprod, Msolve, bilinearfor Line 1136  def MINRES(b, Aprod, Msolve, bilinearfor
1136      ynorm2 = 0      ynorm2 = 0
1137      cs     = -1      cs     = -1
1138      sn     = 0      sn     = 0
1139      w      = b*0.      w      = r*0.
1140      w2     = b*0.      w2     = r*0.
1141      r2     = r1      r2     = r1
1142      eps    = 0.0001      eps    = 0.0001
1143    
1144      #---------------------------------------------------------------------      #---------------------------------------------------------------------
1145      # Main iteration loop.      # Main iteration loop.
1146      # --------------------------------------------------------------------      # --------------------------------------------------------------------
1147      while not stoppingcriterium(rnorm,Anorm*ynorm,'MINRES'):    #  checks ||r|| < (||A|| ||x||) * TOL      while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL
1148    
1149      if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max          if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
1150          iter    = iter  +  1          iter    = iter  +  1
1151    
1152          #-----------------------------------------------------------------          #-----------------------------------------------------------------
# Line 1035  def MINRES(b, Aprod, Msolve, bilinearfor Line 1163  def MINRES(b, Aprod, Msolve, bilinearfor
1163          #-----------------------------------------------------------------          #-----------------------------------------------------------------
1164          s = 1/beta                 # Normalize previous vector (in y).          s = 1/beta                 # Normalize previous vector (in y).
1165          v = s*y                    # v = vk if P = I          v = s*y                    # v = vk if P = I
1166        
1167          y      = Aprod(v)          y      = Aprod(v)
1168        
1169          if iter >= 2:          if iter >= 2:
1170            y = y - (beta/oldb)*r1            y = y - (beta/oldb)*r1
1171    
1172          alfa   = bilinearform(v,y)              # alphak          alfa   = bilinearform(v,y)              # alphak
1173          y      += (- alfa/beta)*r2          y      += (- alfa/beta)*r2
1174          r1     = r2          r1     = r2
1175          r2     = y          r2     = y
1176          y = Msolve(r2)          y = Msolve(r2)
1177          oldb   = beta                           # oldb = betak          oldb   = beta                           # oldb = betak
1178          beta   = bilinearform(y,r2)             # beta = betak+1^2          beta   = bilinearform(y,r2)             # beta = betak+1^2
1179          if beta < 0: raise NegativeNorm,"negative norm."          if beta < 0: raise NegativeNorm("negative norm.")
1180    
1181          beta   = math.sqrt( beta )          beta   = math.sqrt( beta )
1182          tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta          tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
1183            
1184          if iter==1:                 # Initialize a few things.          if iter==1:                 # Initialize a few things.
1185            gmax   = abs( alfa )      # alpha1            gmax   = abs( alfa )      # alpha1
1186            gmin   = gmax             # alpha1            gmin   = gmax             # alpha1
# Line 1060  def MINRES(b, Aprod, Msolve, bilinearfor Line 1188  def MINRES(b, Aprod, Msolve, bilinearfor
1188          # Apply previous rotation Qk-1 to get          # Apply previous rotation Qk-1 to get
1189          #   [deltak epslnk+1] = [cs  sn][dbark    0   ]          #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
1190          #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].          #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
1191        
1192          oldeps = epsln          oldeps = epsln
1193          delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak          delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
1194          gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k          gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
# Line 1070  def MINRES(b, Aprod, Msolve, bilinearfor Line 1198  def MINRES(b, Aprod, Msolve, bilinearfor
1198          # Compute the next plane rotation Qk          # Compute the next plane rotation Qk
1199    
1200          gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak          gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
1201          gamma  = max(gamma,eps)          gamma  = max(gamma,eps)
1202          cs     = gbar / gamma             # ck          cs     = gbar / gamma             # ck
1203          sn     = beta / gamma             # sk          sn     = beta / gamma             # sk
1204          phi    = cs * phibar              # phik          phi    = cs * phibar              # phik
# Line 1078  def MINRES(b, Aprod, Msolve, bilinearfor Line 1206  def MINRES(b, Aprod, Msolve, bilinearfor
1206    
1207          # Update  x.          # Update  x.
1208    
1209          denom = 1/gamma          denom = 1/gamma
1210          w1    = w2          w1    = w2
1211          w2    = w          w2    = w
1212          w     = (v - oldeps*w1 - delta*w2) * denom          w     = (v - oldeps*w1 - delta*w2) * denom
1213          x     +=  phi*w          x     +=  phi*w
1214    
# Line 1095  def MINRES(b, Aprod, Msolve, bilinearfor Line 1223  def MINRES(b, Aprod, Msolve, bilinearfor
1223    
1224          # Estimate various norms and test for convergence.          # Estimate various norms and test for convergence.
1225    
1226          Anorm  = math.sqrt( tnorm2 )          Anorm  = math.sqrt( tnorm2 )
1227          ynorm  = math.sqrt( ynorm2 )          ynorm  = math.sqrt( ynorm2 )
1228    
1229          rnorm  = phibar          rnorm  = phibar
1230    
1231      return x      return x
1232    
1233  def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  def TFQMR(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1234      """
1235      Solver for
1236    
1237  # TFQMR solver for linear systems    *Ax=b*
 #  
 #  
 # initialization  
 #  
   errtol = math.sqrt(bilinearform(b,b))  
   norm_b=errtol  
   kmax  = iter_max  
   error = []  
   
   if math.sqrt(bilinearform(x,x)) != 0.0:  
     r = b - Aprod(x)  
   else:  
     r = b  
1238    
1239    r=Msolve(r)    with a general operator A (more details required!).
1240      It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).
1241    
1242      The iteration is terminated if
1243    
1244      *|r| <= atol+rtol*|r0|*
1245    
1246      where *r0* is the initial residual and *|.|* is the energy norm. In fact
1247    
1248      *|r| = sqrt( bilinearform(r,r))*
1249    
1250      :param r: initial residual *r=b-Ax*. ``r`` is altered.
1251      :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1252      :param x: an initial guess for the solution
1253      :type x: same like ``r``
1254      :param Aprod: returns the value Ax
1255      :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
1256                   argument ``x``. The returned object needs to be of the same type
1257                   like argument ``r``.
1258      :param bilinearform: inner product ``<x,r>``
1259      :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
1260                          type like argument ``x`` and ``r``. The returned value is
1261                          a ``float``.
1262      :param atol: absolute tolerance
1263      :type atol: non-negative ``float``
1264      :param rtol: relative tolerance
1265      :type rtol: non-negative ``float``
1266      :param iter_max: maximum number of iteration steps
1267      :type iter_max: ``int``
1268      :rtype: ``tuple``
1269      :warning: ``r`` and ``x`` are altered.
1270      """
1271    u1=0    u1=0
1272    u2=0    u2=0
1273    y1=0    y1=0
1274    y2=0    y2=0
1275    
1276    w = r    w = r
1277    y1 = r    y1 = r
1278    iter = 0    iter = 0
1279    d = 0    d = 0
1280        v = Aprod(y1)
   v = Msolve(Aprod(y1))  
1281    u1 = v    u1 = v
1282      
1283    theta = 0.0;    theta = 0.0;
1284    eta = 0.0;    eta = 0.0;
1285    tau = math.sqrt(bilinearform(r,r))    rho=bilinearform(r,r)
1286    error = [ error, tau ]    if rho < 0: raise NegativeNorm("negative norm.")
1287    rho = tau * tau    tau = math.sqrt(rho)
1288    m=1    norm_r0=tau
1289  #    while tau>atol+rtol*norm_r0:
1290  #  TFQMR iteration      if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
 #  
 #  while ( iter < kmax-1 ):  
     
   while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):  
     if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max  
1291    
1292      sigma = bilinearform(r,v)      sigma = bilinearform(r,v)
1293        if sigma == 0.0: raise IterationBreakDown('TFQMR breakdown, sigma=0')
     if ( sigma == 0.0 ):  
       raise 'TFQMR breakdown, sigma=0'  
       
1294    
1295      alpha = rho / sigma      alpha = rho / sigma
1296    
# Line 1162  def TFQMR(b, Aprod, Msolve, bilinearform Line 1300  def TFQMR(b, Aprod, Msolve, bilinearform
1300  #  #
1301        if ( j == 1 ):        if ( j == 1 ):
1302          y2 = y1 - alpha * v          y2 = y1 - alpha * v
1303          u2 = Msolve(Aprod(y2))          u2 = Aprod(y2)
1304    
1305        m = 2 * (iter+1) - 2 + (j+1)        m = 2 * (iter+1) - 2 + (j+1)
1306        if j==0:        if j==0:
1307           w = w - alpha * u1           w = w - alpha * u1
1308           d = y1 + ( theta * theta * eta / alpha ) * d           d = y1 + ( theta * theta * eta / alpha ) * d
1309        if j==1:        if j==1:
# Line 1180  def TFQMR(b, Aprod, Msolve, bilinearform Line 1318  def TFQMR(b, Aprod, Msolve, bilinearform
1318  #  #
1319  #  Try to terminate the iteration at each pass through the loop  #  Try to terminate the iteration at each pass through the loop
1320  #  #
1321       # if ( tau * math.sqrt ( m + 1 ) <= errtol ):      if rho == 0.0: raise IterationBreakDown('TFQMR breakdown, rho=0')
      #   error = [ error, tau ]  
      #   total_iters = iter  
      #   break  
         
   
     if ( rho == 0.0 ):  
       raise 'TFQMR breakdown, rho=0'  
       
1322    
1323      rhon = bilinearform(r,w)      rhon = bilinearform(r,w)
1324      beta = rhon / rho;      beta = rhon / rho;
1325      rho = rhon;      rho = rhon;
1326      y1 = w + beta * y2;      y1 = w + beta * y2;
1327      u1 = Msolve(Aprod(y1))      u1 = Aprod(y1)
1328      v = u1 + beta * ( u2 + beta * v )      v = u1 + beta * ( u2 + beta * v )
1329      error = [ error, tau ]  
1330      total_iters = iter      iter += 1
       
     iter = iter + 1  
1331    
1332    return x    return x
1333    
# Line 1208  def TFQMR(b, Aprod, Msolve, bilinearform Line 1336  def TFQMR(b, Aprod, Msolve, bilinearform
1336    
1337  class ArithmeticTuple(object):  class ArithmeticTuple(object):
1338     """     """
1339     tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.     Tuple supporting inplace update x+=y and scaling x=a*y where ``x,y`` is an
1340       ArithmeticTuple and ``a`` is a float.
1341    
1342     example of usage:     Example of usage::
1343    
1344     from esys.escript import Data         from esys.escript import Data
1345     from numarray import array         from numpy import array
1346     a=Data(...)         a=Data(...)
1347     b=array([1.,4.])         b=array([1.,4.])
1348     x=ArithmeticTuple(a,b)         x=ArithmeticTuple(a,b)
1349     y=5.*x         y=5.*x
1350    
1351     """     """
1352     def __init__(self,*args):     def __init__(self,*args):
1353         """         """
1354         initialize object with elements args.         Initializes object with elements ``args``.
1355    
1356         @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)         :param args: tuple of objects that support inplace add (x+=y) and
1357                        scaling (x=a*y)
1358         """         """
1359         self.__items=list(args)         self.__items=list(args)
1360    
1361     def __len__(self):     def __len__(self):
1362         """         """
1363         number of items         Returns the number of items.
1364    
1365         @return: number of items         :return: number of items
1366         @rtype: C{int}         :rtype: ``int``
1367         """         """
1368         return len(self.__items)         return len(self.__items)
1369    
1370     def __getitem__(self,index):     def __getitem__(self,index):
1371         """         """
1372         get an item         Returns item at specified position.
1373    
1374         @param index: item to be returned         :param index: index of item to be returned
1375         @type index: C{int}         :type index: ``int``
1376         @return: item with index C{index}         :return: item with index ``index``
1377         """         """
1378         return self.__items.__getitem__(index)         return self.__items.__getitem__(index)
1379    
1380     def __mul__(self,other):     def __mul__(self,other):
1381         """         """
1382         scaling from the right         Scales by ``other`` from the right.
1383    
1384         @param other: scaling factor         :param other: scaling factor
1385         @type other: C{float}         :type other: ``float``
1386         @return: itemwise self*other         :return: itemwise self*other
1387         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1388         """         """
1389         out=[]         out=[]
1390         try:           try:
1391             l=len(other)             l=len(other)
1392             if l!=len(self):             if l!=len(self):
1393                 raise ValueError,"length of of arguments don't match."                 raise ValueError("length of arguments don't match.")
1394             for i in range(l): out.append(self[i]*other[i])             for i in range(l): out.append(self[i]*other[i])
1395         except TypeError:         except TypeError:
1396         for i in range(len(self)): out.append(self[i]*other)             for i in range(len(self)): out.append(self[i]*other)
1397         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1398    
1399     def __rmul__(self,other):     def __rmul__(self,other):
1400         """         """
1401         scaling from the left         Scales by ``other`` from the left.
1402    
1403         @param other: scaling factor         :param other: scaling factor
1404         @type other: C{float}         :type other: ``float``
1405         @return: itemwise other*self         :return: itemwise other*self
1406         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1407         """         """
1408         out=[]         out=[]
1409         try:           try:
1410             l=len(other)             l=len(other)
1411             if l!=len(self):             if l!=len(self):
1412                 raise ValueError,"length of of arguments don't match."                 raise ValueError("length of arguments don't match.")
1413             for i in range(l): out.append(other[i]*self[i])             for i in range(l): out.append(other[i]*self[i])
1414         except TypeError:         except TypeError:
1415         for i in range(len(self)): out.append(other*self[i])             for i in range(len(self)): out.append(other*self[i])
1416         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1417    
1418     def __div__(self,other):     def __div__(self,other):
1419         """         """
1420         dividing from the right         Scales by (1/``other``) from the right.
1421    
1422         @param other: scaling factor         :param other: scaling factor
1423         @type other: C{float}         :type other: ``float``
1424         @return: itemwise self/other         :return: itemwise self/other
1425         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1426         """         """
1427         return self*(1/other)         return self*(1/other)
1428    
1429     def __rdiv__(self,other):     def __rdiv__(self,other):
1430         """         """
1431         dividing from the left         Scales by (1/``other``) from the left.
1432    
1433         @param other: scaling factor         :param other: scaling factor
1434         @type other: C{float}         :type other: ``float``
1435         @return: itemwise other/self         :return: itemwise other/self
1436         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1437         """         """
1438         out=[]         out=[]
1439         try:           try:
1440             l=len(other)             l=len(other)
1441             if l!=len(self):             if l!=len(self):
1442                 raise ValueError,"length of of arguments don't match."                 raise ValueError("length of arguments don't match.")
1443             for i in range(l): out.append(other[i]/self[i])             for i in range(l): out.append(other[i]/self[i])
1444         except TypeError:         except TypeError:
1445         for i in range(len(self)): out.append(other/self[i])             for i in range(len(self)): out.append(other/self[i])
1446         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1447      
1448     def __iadd__(self,other):     def __iadd__(self,other):
1449         """         """
1450         in-place add of other to self         Inplace addition of ``other`` to self.
1451    
1452         @param other: increment         :param other: increment
1453         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1454         """         """
1455         if len(self) != len(other):         if len(self) != len(other):
1456             raise ValueError,"tuple length must match."             raise ValueError("tuple lengths must match.")
1457         for i in range(len(self)):         for i in range(len(self)):
1458             self.__items[i]+=other[i]             self.__items[i]+=other[i]
1459         return self         return self
1460    
1461     def __add__(self,other):     def __add__(self,other):
1462         """         """
1463         add other to self         Adds ``other`` to self.
1464    
1465         @param other: increment         :param other: increment
1466         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1467         """         """
1468         out=[]         out=[]
1469         try:           try:
1470             l=len(other)             l=len(other)
1471             if l!=len(self):             if l!=len(self):
1472                 raise ValueError,"length of of arguments don't match."                 raise ValueError("length of arguments don't match.")
1473             for i in range(l): out.append(self[i]+other[i])             for i in range(l): out.append(self[i]+other[i])
1474         except TypeError:         except TypeError:
1475         for i in range(len(self)): out.append(self[i]+other)             for i in range(len(self)): out.append(self[i]+other)
1476         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1477    
1478     def __sub__(self,other):     def __sub__(self,other):
1479         """         """
1480         subtract other from self         Subtracts ``other`` from self.
1481    
1482         @param other: increment         :param other: decrement
1483         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1484         """         """
1485         out=[]         out=[]
1486         try:           try:
1487             l=len(other)             l=len(other)
1488             if l!=len(self):             if l!=len(self):
1489                 raise ValueError,"length of of arguments don't match."                 raise ValueError("length of arguments don't match.")
1490             for i in range(l): out.append(self[i]-other[i])             for i in range(l): out.append(self[i]-other[i])
1491         except TypeError:         except TypeError:
1492         for i in range(len(self)): out.append(self[i]-other)             for i in range(len(self)): out.append(self[i]-other)
1493         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1494      
1495     def __isub__(self,other):     def __isub__(self,other):
1496         """         """
1497         subtract other from self         Inplace subtraction of ``other`` from self.
1498    
1499         @param other: increment         :param other: decrement
1500         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1501         """         """
1502         if len(self) != len(other):         if len(self) != len(other):
1503             raise ValueError,"tuple length must match."             raise ValueError("tuple length must match.")
1504         for i in range(len(self)):         for i in range(len(self)):
1505             self.__items[i]-=other[i]             self.__items[i]-=other[i]
1506         return self         return self
1507    
1508     def __neg__(self):     def __neg__(self):
1509         """         """
1510         negate         Negates values.
   
1511         """         """
1512         out=[]         out=[]
1513         for i in range(len(self)):         for i in range(len(self)):
# Line 1388  class ArithmeticTuple(object): Line 1517  class ArithmeticTuple(object):
1517    
1518  class HomogeneousSaddlePointProblem(object):  class HomogeneousSaddlePointProblem(object):
1519        """        """
1520        This provides a framwork for solving linear homogeneous saddle point problem of the form        This class provides a framework for solving linear homogeneous saddle
1521          point problems of the form::
              Av+B^*p=f  
              Bv    =0  
1522    
1523        for the unknowns v and p and given operators A and B and given right hand side f.            *Av+B^*p=f*
1524        B^* is the adjoint operator of B is the given inner product.            *Bv     =0*
1525    
1526          for the unknowns *v* and *p* and given operators *A* and *B* and
1527          given right hand side *f*. *B^** is the adjoint operator of *B*.
1528          *A* may depend weakly on *v* and *p*.
1529        """        """
1530        def __init__(self,**kwargs):        def __init__(self, **kwargs):
         self.setTolerance()  
         self.setToleranceReductionFactor()  
   
       def initialize(self):  
1531          """          """
1532          initialize the problem (overwrite)          initializes the saddle point problem
1533          """          """
1534          pass          self.resetControlParameters()
1535        def B(self,v):          self.setTolerance()
1536            self.setAbsoluteTolerance()
1537          def resetControlParameters(self, K_p=1., K_v=1., rtol_max=0.01, rtol_min = 1.e-7, chi_max=0.5, reduction_factor=0.3, theta = 0.1):
1538           """           """
1539           returns Bv (overwrite)           sets a control parameter
          @rtype: equal to the type of p  
1540    
1541           @note: boundary conditions on p should be zero!           :param K_p: initial value for constant to adjust pressure tolerance
1542             :type K_p: ``float``
1543             :param K_v: initial value for constant to adjust velocity tolerance
1544             :type K_v: ``float``
1545             :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1546             :type rtol_max: ``float``
1547             :param chi_max: maximum tolerable converegence rate.
1548             :type chi_max: ``float``
1549             :param reduction_factor: reduction factor for adjustment factors.
1550             :type reduction_factor: ``float``
1551           """           """
1552           pass           self.setControlParameter(K_p, K_v, rtol_max, rtol_min, chi_max, reduction_factor, theta)
1553    
1554        def inner(self,p0,p1):        def setControlParameter(self,K_p=None, K_v=None, rtol_max=None, rtol_min=None, chi_max=None, reduction_factor=None, theta=None):
1555           """           """
1556           returns inner product of two element p0 and p1  (overwrite)           sets a control parameter
           
          @type p0: equal to the type of p  
          @type p1: equal to the type of p  
          @rtype: C{float}  
1557    
          @rtype: equal to the type of p  
          """  
          pass  
1558    
1559        def solve_A(self,u,p):           :param K_p: initial value for constant to adjust pressure tolerance
1560             :type K_p: ``float``
1561             :param K_v: initial value for constant to adjust velocity tolerance
1562             :type K_v: ``float``
1563             :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1564             :type rtol_max: ``float``
1565             :param chi_max: maximum tolerable converegence rate.
1566             :type chi_max: ``float``
1567             :type reduction_factor: ``float``
1568           """           """
1569           solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)           if not K_p == None:
1570                if K_p<1:
1571                   raise ValueError("K_p need to be greater or equal to 1.")
1572             else:
1573                K_p=self.__K_p
1574    
1575             if not K_v == None:
1576                if K_v<1:
1577                   raise ValueError("K_v need to be greater or equal to 1.")
1578             else:
1579                K_v=self.__K_v
1580    
1581             if not rtol_max == None:
1582                if rtol_max<=0 or rtol_max>=1:
1583                   raise ValueError("rtol_max needs to be positive and less than 1.")
1584             else:
1585                rtol_max=self.__rtol_max
1586    
1587             if not rtol_min == None:
1588                if rtol_min<=0 or rtol_min>=1:
1589                   raise ValueError("rtol_min needs to be positive and less than 1.")
1590             else:
1591                rtol_min=self.__rtol_min
1592    
1593             if not chi_max == None:
1594                if chi_max<=0 or chi_max>=1:
1595                   raise ValueError("chi_max needs to be positive and less than 1.")
1596             else:
1597                chi_max = self.__chi_max
1598    
1599             if not reduction_factor == None:
1600                if reduction_factor<=0 or reduction_factor>1:
1601                   raise ValueError("reduction_factor need to be between zero and one.")
1602             else:
1603                reduction_factor=self.__reduction_factor
1604    
1605             if not theta == None:
1606                if theta<=0 or theta>1:
1607                   raise ValueError("theta need to be between zero and one.")
1608             else:
1609                theta=self.__theta
1610    
1611             if rtol_min>=rtol_max:
1612                 raise ValueError("rtol_max = %e needs to be greater than rtol_min = %e"%(rtol_max,rtol_min))
1613             self.__chi_max = chi_max
1614             self.__rtol_max = rtol_max
1615             self.__K_p = K_p
1616             self.__K_v = K_v
1617             self.__reduction_factor = reduction_factor
1618             self.__theta = theta
1619             self.__rtol_min=rtol_min
1620    
1621           @rtype: equal to the type of v        #=============================================================
1622           @note: boundary conditions on v should be zero!        def inner_pBv(self,p,Bv):
1623           """           """
1624           pass           Returns inner product of element p and Bv (overwrite).
1625    
1626        def solve_prec(self,p):           :param p: a pressure increment
1627             :param Bv: a residual
1628             :return: inner product of element p and Bv
1629             :rtype: ``float``
1630             :note: used if PCG is applied.
1631           """           """
1632           provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)           raise NotImplementedError("no inner product for p and Bv implemented.")
1633    
1634           @rtype: equal to the type of p        def inner_p(self,p0,p1):
1635           """           """
1636           pass           Returns inner product of p0 and p1 (overwrite).
1637    
1638        def stoppingcriterium(self,Bv,v,p):           :param p0: a pressure
1639             :param p1: a pressure
1640             :return: inner product of p0 and p1
1641             :rtype: ``float``
1642           """           """
1643           returns a True if iteration is terminated. (overwrite)           raise NotImplementedError("no inner product for p implemented.")
1644      
1645          def norm_v(self,v):
1646             """
1647             Returns the norm of v (overwrite).
1648    
1649           @rtype: C{bool}           :param v: a velovity
1650             :return: norm of v
1651             :rtype: non-negative ``float``
1652             """
1653             raise NotImplementedError("no norm of v implemented.")
1654          def getDV(self, p, v, tol):
1655           """           """
1656           pass           return a correction to the value for a given v and a given p with accuracy `tol` (overwrite)
               
       def __inner(self,p,r):  
          return self.inner(p,r[1])  
1657    
1658        def __inner_p(self,p1,p2):           :param p: pressure
1659           return self.inner(p1,p2)           :param v: pressure
1660                   :return: dv given as *dv= A^{-1} (f-A v-B^*p)*
1661        def __inner_a(self,a1,a2):           :note: Only *A* may depend on *v* and *p*
1662           return self.inner_a(a1,a2)           """
1663             raise NotImplementedError("no dv calculation implemented.")
1664    
1665        def __inner_a1(self,a1,a2):          
1666           return self.inner(a1[1],a2[1])        def Bv(self,v, tol):
1667            """
1668            Returns Bv with accuracy `tol` (overwrite)
1669    
1670        def __stoppingcriterium(self,norm_r,r,p):          :rtype: equal to the type of p
1671            return self.stoppingcriterium(r[1],r[0],p)          :note: boundary conditions on p should be zero!
1672            """
1673            raise NotImplementedError("no operator B implemented.")
1674    
1675        def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):        def norm_Bv(self,Bv):
1676            return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)          """
1677            Returns the norm of Bv (overwrite).
1678    
1679        def setTolerance(self,tolerance=1.e-8):          :rtype: equal to the type of p
1680                self.__tol=tolerance          :note: boundary conditions on p should be zero!
1681        def getTolerance(self):          """
1682                return self.__tol          raise NotImplementedError("no norm of Bv implemented.")
       def setToleranceReductionFactor(self,reduction=0.01):  
               self.__reduction=reduction  
       def getSubProblemTolerance(self):  
               return self.__reduction*self.getTolerance()  
   
       def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):  
               """  
               solves the saddle point problem using initial guesses v and p.  
   
               @param max_iter: maximum number of iteration steps.  
               """  
               self.verbose=verbose  
               self.show_details=show_details and self.verbose  
   
               # assume p is known: then v=A^-1(f-B^*p)  
               # which leads to BA^-1B^*p = BA^-1f    
   
           # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)        
           self.__z=v+self.solve_A(v,p*0)  
               Bz=self.B(self.__z)  
               #  
           #   solve BA^-1B^*p = Bz  
               #  
               #  
               #  
               self.iter=0  
           if solver=='GMRES':        
                 if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter  
                 p=GMRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)  
                 # solve Au=f-B^*p  
                 #       A(u-v)=f-B^*p-Av  
                 #       u=v+(u-v)  
         u=v+self.solve_A(v,p)  
   
           if solver=='TFQMR':        
                 if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter  
                 p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)  
                 # solve Au=f-B^*p  
                 #       A(u-v)=f-B^*p-Av  
                 #       u=v+(u-v)  
         u=v+self.solve_A(v,p)  
   
           if solver=='MINRES':        
                 if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter  
                 p=MINRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)  
                 # solve Au=f-B^*p  
                 #       A(u-v)=f-B^*p-Av  
                 #       u=v+(u-v)  
         u=v+self.solve_A(v,p)  
                 
           if solver=='GMRESC':        
                 if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter  
                 p0=self.solve_prec1(Bz)  
             #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))  
                 #p-=alfa  
                 x=GMRES(ArithmeticTuple(self.__z*1.,p0*1),self.__Anext,self.__Mempty,self.__inner_a,self.__stoppingcriterium2,iter_max=max_iter, x=ArithmeticTuple(v*1,p*1),iter_restart=20)  
                 #x=NewtonGMRES(ArithmeticTuple(self.__z*1.,p0*1),self.__Aprod_Newton2,self.__Mempty,self.__inner_a,self.__stoppingcriterium2,iter_max=max_iter, x=ArithmeticTuple(v*1,p*1),atol=0,rtol=self.getTolerance())  
   
                 # solve Au=f-B^*p  
                 #       A(u-v)=f-B^*p-Av  
                 #       u=v+(u-v)  
             p=x[1]  
         u=v+self.solve_A(v,p)        
         #p=x[1]  
         #u=x[0]  
   
               if solver=='PCG':  
                 #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv  
                 #  
                 #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)  
                 #                           A(v-z)= f -Az - B^*p (v-z=0 on fixed_u_mask)  
                 if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter  
                 p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)  
             u=r[0]    
                 # print "DIFF=",util.integrate(p)  
   
               # print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)  
   
           return u,p  
   
       def __Msolve(self,r):  
           return self.solve_prec(r[1])  
   
       def __Msolve2(self,r):  
           return self.solve_prec(r*1)  
   
       def __Mempty(self,r):  
           return r  
   
   
       def __Aprod(self,p):  
           # return BA^-1B*p  
           #solve Av =B^*p as Av =f-Az-B^*(-p)  
           v=self.solve_A(self.__z,-p)  
           return ArithmeticTuple(v, self.B(v))  
   
       def __Anext(self,x):  
           # return next v,p  
           #solve Av =-B^*p as Av =f-Az-B^*p  
   
       pc=x[1]  
           v=self.solve_A(self.__z,-pc)  
       p=self.solve_prec1(self.B(v))  
   
           return ArithmeticTuple(v,p)  
   
   
       def __Aprod2(self,p):  
           # return BA^-1B*p  
           #solve Av =B^*p as Av =f-Az-B^*(-p)  
       v=self.solve_A(self.__z,-p)  
           return self.B(v)  
   
       def __Aprod_Newton(self,p):  
           # return BA^-1B*p - Bz  
           #solve Av =-B^*p as Av =f-Az-B^*p  
       v=self.solve_A(self.__z,-p)  
           return self.B(v-self.__z)  
   
       def __Aprod_Newton2(self,x):  
           # return BA^-1B*p - Bz  
           #solve Av =-B^*p as Av =f-Az-B^*p  
           pc=x[1]  
       v=self.solve_A(self.__z,-pc)  
           p=self.solve_prec1(self.B(v-self.__z))  
           return ArithmeticTuple(v,p)  
1683    
1684          def solve_AinvBt(self,dp, tol):
1685             """
1686             Solves *A dv=B^*dp* with accuracy `tol`
1687    
1688  def MaskFromBoundaryTag(domain,*tags):           :param dp: a pressure increment
1689     """           :return: the solution of *A dv=B^*dp*
1690     creates a mask on the Solution(domain) function space which one for samples           :note: boundary conditions on dv should be zero! *A* is the operator used in ``getDV`` and must not be altered.
1691     that touch the boundary tagged by tags.           """
1692             raise NotImplementedError("no operator A implemented.")
1693    
1694     usage: m=MaskFromBoundaryTag(domain,"left", "right")        def solve_prec(self,Bv, tol):
1695             """
1696             Provides a preconditioner for *(BA^{-1}B^ * )* applied to Bv with accuracy `tol`
1697    
1698     @param domain: a given domain           :rtype: equal to the type of p
1699     @type domain: L{escript.Domain}           :note: boundary conditions on p should be zero!
1700     @param tags: boundray tags           """
1701     @type tags: C{str}           raise NotImplementedError("no preconditioner for Schur complement implemented.")
1702     @return: a mask which marks samples that are touching the boundary tagged by any of the given tags.        #=============================================================
1703     @rtype: L{escript.Data} of rank 0        def __Aprod_PCG(self,dp):
1704     """            dv=self.solve_AinvBt(dp, self.__subtol)
1705     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)            return ArithmeticTuple(dv,self.Bv(dv, self.__subtol))
1706     d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))  
1707     for t in tags: d.setTaggedValue(t,1.)        def __inner_PCG(self,p,r):
1708     pde.setValue(y=d)           return self.inner_pBv(p,r[1])
1709     return util.whereNonZero(pde.getRightHandSide())  
1710  #============================================================================================================================================        def __Msolve_PCG(self,r):
1711  class SaddlePointProblem(object):            return self.solve_prec(r[1], self.__subtol)
1712     """        #=============================================================
1713     This implements a solver for a saddlepoint problem        def __Aprod_GMRES(self,p):
1714              return self.solve_prec(self.Bv(self.solve_AinvBt(p, self.__subtol), self.__subtol), self.__subtol)
1715          def __inner_GMRES(self,p0,p1):
1716             return self.inner_p(p0,p1)
1717    
1718          #=============================================================
1719          def norm_p(self,p):
1720              """
1721              calculates the norm of ``p``
1722              
1723              :param p: a pressure
1724              :return: the norm of ``p`` using the inner product for pressure
1725              :rtype: ``float``
1726              """
1727              f=self.inner_p(p,p)
1728              if f<0: raise ValueError("negative pressure norm.")
1729              return math.sqrt(f)
1730          
1731          def solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1732             """
1733             Solves the saddle point problem using initial guesses v and p.
1734    
1735     M{f(u,p)=0}           :param v: initial guess for velocity
1736     M{g(u)=0}           :param p: initial guess for pressure
1737             :type v: `Data`
1738             :type p: `Data`
1739             :param usePCG: indicates the usage of the PCG rather than GMRES scheme.
1740             :param max_iter: maximum number of iteration steps per correction
1741                              attempt
1742             :param verbose: if True, shows information on the progress of the
1743                             saddlepoint problem solver.
1744             :param iter_restart: restart the iteration after ``iter_restart`` steps
1745                                  (only used if useUzaw=False)
1746             :type usePCG: ``bool``
1747             :type max_iter: ``int``
1748             :type verbose: ``bool``
1749             :type iter_restart: ``int``
1750             :rtype: ``tuple`` of `Data` objects
1751             :note: typically this method is overwritten by a subclass. It provides a wrapper for the ``_solve`` method.
1752             """
1753             return self._solve(v=v,p=p,max_iter=max_iter,verbose=verbose, usePCG=usePCG, iter_restart=iter_restart, max_correction_steps=max_correction_steps)
1754    
1755     for u and p. The problem is solved with an inexact Uszawa scheme for p:        def _solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1756             """
1757             see `_solve` method.
1758             """
1759             self.verbose=verbose
1760             rtol=self.getTolerance()
1761             atol=self.getAbsoluteTolerance()
1762    
1763             K_p=self.__K_p
1764             K_v=self.__K_v
1765             correction_step=0
1766             converged=False
1767             chi=None
1768             eps=None
1769    
1770             if self.verbose: print(("HomogeneousSaddlePointProblem: start iteration: rtol= %e, atol=%e"%(rtol, atol)))
1771             while not converged:
1772    
1773                 # get tolerance for velecity increment:
1774                 if chi == None:
1775                    rtol_v=self.__rtol_max
1776                 else:
1777                    rtol_v=min(chi/K_v,self.__rtol_max)
1778                 rtol_v=max(rtol_v, self.__rtol_min)
1779                 if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: rtol_v= %e"%(correction_step,rtol_v)))
1780                 # get velocity increment:
1781                 dv1=self.getDV(p,v,rtol_v)
1782                 v1=v+dv1
1783                 Bv1=self.Bv(v1, rtol_v)
1784                 norm_Bv1=self.norm_Bv(Bv1)
1785                 norm_dv1=self.norm_v(dv1)
1786                 if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: norm_Bv1 = %e, norm_dv1 = %e"%(correction_step, norm_Bv1, norm_dv1)))
1787                 if norm_dv1*self.__theta < norm_Bv1:
1788                    # get tolerance for pressure increment:
1789                    large_Bv1=True
1790                    if chi == None or eps == None:
1791                       rtol_p=self.__rtol_max
1792                    else:
1793                       rtol_p=min(chi**2*eps/K_p/norm_Bv1, self.__rtol_max)
1794                    self.__subtol=max(rtol_p**2, self.__rtol_min)
1795                    if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: rtol_p= %e"%(correction_step,rtol_p)))
1796                    # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1
1797                    if usePCG:
1798                        dp,r,a_norm=PCG(ArithmeticTuple(v1,Bv1),self.__Aprod_PCG,0*p,self.__Msolve_PCG,self.__inner_PCG,atol=0, rtol=rtol_p,iter_max=max_iter, verbose=self.verbose)
1799                        v2=r[0]
1800                        Bv2=r[1]
1801                    else:
1802                        # don't use!!!!
1803                        dp=GMRES(self.solve_prec(Bv1,self.__subtol),self.__Aprod_GMRES, 0*p, self.__inner_GMRES,atol=0, rtol=rtol_p,iter_max=max_iter, iter_restart=iter_restart, verbose=self.verbose)
1804                        dv2=self.solve_AinvBt(dp, self.__subtol)
1805                        v2=v1-dv2
1806                        Bv2=self.Bv(v2, self.__subtol)
1807                    p2=p+dp
1808                 else:
1809                    large_Bv1=False
1810                    v2=v1
1811                    p2=p
1812                 # update business:
1813                 norm_dv2=self.norm_v(v2-v)
1814                 norm_v2=self.norm_v(v2)
1815                 if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: v2 = %e, norm_dv2 = %e"%(correction_step, norm_v2, self.norm_v(v2-v))))
1816                 eps, eps_old = max(norm_Bv1, norm_dv2), eps
1817                 if eps_old == None:
1818                      chi, chi_old = None, chi
1819                 else:
1820                      chi, chi_old = min(eps/ eps_old, self.__chi_max), chi
1821                 if eps != None:
1822                     if chi !=None:
1823                        if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: convergence rate = %e, correction = %e"%(correction_step,chi, eps)))
1824                     else:
1825                        if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: correction = %e"%(correction_step, eps)))
1826                 if eps <= rtol*norm_v2+atol :
1827                     converged = True
1828                 else:
1829                     if correction_step>=max_correction_steps:
1830                          raise CorrectionFailed("Given up after %d correction steps."%correction_step)
1831                     if chi_old!=None:
1832                        K_p=max(1,self.__reduction_factor*K_p,(chi-chi_old)/chi_old**2*K_p)
1833                        K_v=max(1,self.__reduction_factor*K_v,(chi-chi_old)/chi_old**2*K_p)
1834                        if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: new adjustment factor K = %e"%(correction_step,K_p)))
1835                 correction_step+=1
1836                 v,p =v2, p2
1837             if self.verbose: print(("HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step))
1838             return v,p
1839          #========================================================================
1840          def setTolerance(self,tolerance=1.e-4):
1841             """
1842             Sets the relative tolerance for (v,p).
1843    
1844     M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}           :param tolerance: tolerance to be used
1845     M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}           :type tolerance: non-negative ``float``
1846             """
1847             if tolerance<0:
1848                 raise ValueError("tolerance must be positive.")
1849             self.__rtol=tolerance
1850    
1851     where Q_f is an approximation of the Jacobiean A_f of f with respect to u  and Q_f is an approximation of        def getTolerance(self):
1852     A_g A_f^{-1} A_g with A_g is the jacobiean of g with respect to p. As a the construction of a 'proper'           """
1853     Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays           Returns the relative tolerance.
    in fact the role of a preconditioner.  
    """  
    def __init__(self,verbose=False,*args):  
        """  
        initializes the problem  
1854    
1855         @param verbose: switches on the printing out some information           :return: relative tolerance
1856         @type verbose: C{bool}           :rtype: ``float``
1857         @note: this method may be overwritten by a particular saddle point problem           """
1858         """           return self.__rtol
        print "SaddlePointProblem should not be used anymore!"  
        if not isinstance(verbose,bool):  
             raise TypeError("verbose needs to be of type bool.")  
        self.__verbose=verbose  
        self.relaxation=1.  
        DeprecationWarning("SaddlePointProblem should not be used anymore.")  
1859    
1860     def trace(self,text):        def setAbsoluteTolerance(self,tolerance=0.):
1861         """           """
1862         prints text if verbose has been set           Sets the absolute tolerance.
1863    
1864         @param text: a text message           :param tolerance: tolerance to be used
1865         @type text: C{str}           :type tolerance: non-negative ``float``
1866         """           """
1867         if self.__verbose: print "%s: %s"%(str(self),text)           if tolerance<0:
1868                 raise ValueError("tolerance must be non-negative.")
1869             self.__atol=tolerance
1870    
1871     def solve_f(self,u,p,tol=1.e-8):        def getAbsoluteTolerance(self):
1872         """           """
1873         solves           Returns the absolute tolerance.
1874    
1875         A_f du = f(u,p)           :return: absolute tolerance
1876             :rtype: ``float``
1877             """
1878             return self.__atol
1879    
        with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.  
1880    
1881         @param u: current approximation of u  def MaskFromBoundaryTag(domain,*tags):
1882         @type u: L{escript.Data}     """
1883         @param p: current approximation of p     Creates a mask on the Solution(domain) function space where the value is
1884         @type p: L{escript.Data}     one for samples that touch the boundary tagged by tags.
        @param tol: tolerance expected for du  
        @type tol: C{float}  
        @return: increment du  
        @rtype: L{escript.Data}  
        @note: this method has to be overwritten by a particular saddle point problem  
        """  
        pass  
1885    
1886     def solve_g(self,u,tol=1.e-8):     Usage: m=MaskFromBoundaryTag(domain, "left", "right")
        """  
        solves  
1887    
1888         Q_g dp = g(u)     :param domain: domain to be used
1889       :type domain: `escript.Domain`
1890       :param tags: boundary tags
1891       :type tags: ``str``
1892       :return: a mask which marks samples that are touching the boundary tagged
1893                by any of the given tags
1894       :rtype: `escript.Data` of rank 0
1895       """
1896       pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1897       d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
1898       for t in tags: d.setTaggedValue(t,1.)
1899       pde.setValue(y=d)
1900       return util.whereNonZero(pde.getRightHandSide())
1901    
1902         with Q_g is a preconditioner for A_g A_f^{-1} A_g with  A_g is the jacobiean of g with respect to p.  def MaskFromTag(domain,*tags):
1903       """
1904       Creates a mask on the Solution(domain) function space where the value is
1905       one for samples that touch regions tagged by tags.
1906    
1907         @param u: current approximation of u     Usage: m=MaskFromTag(domain, "ham")
        @type u: L{escript.Data}  
        @param tol: tolerance expected for dp  
        @type tol: C{float}  
        @return: increment dp  
        @rtype: L{escript.Data}  
        @note: this method has to be overwritten by a particular saddle point problem  
        """  
        pass  
1908    
1909     def inner(self,p0,p1):     :param domain: domain to be used
1910         """     :type domain: `escript.Domain`
1911         inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)     :param tags: boundary tags
1912         @return: inner product of p0 and p1     :type tags: ``str``
1913         @rtype: C{float}     :return: a mask which marks samples that are touching the boundary tagged
1914         """              by any of the given tags
1915         pass     :rtype: `escript.Data` of rank 0
1916       """
1917       pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1918       d=escript.Scalar(0.,escript.Function(domain))
1919       for t in tags: d.setTaggedValue(t,1.)
1920       pde.setValue(Y=d)
1921       return util.whereNonZero(pde.getRightHandSide())
1922    
    subiter_max=3  
    def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):  
         """  
         runs the solver  
1923    
         @param u0: initial guess for C{u}  
         @type u0: L{esys.escript.Data}  
         @param p0: initial guess for C{p}  
         @type p0: L{esys.escript.Data}  
         @param tolerance: tolerance for relative error in C{u} and C{p}  
         @type tolerance: positive C{float}  
         @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}  
         @type tolerance_u: positive C{float}  
         @param iter_max: maximum number of iteration steps.  
         @type iter_max: C{int}  
         @param accepted_reduction: if the norm  g cannot be reduced by C{accepted_reduction} backtracking to adapt the  
                                    relaxation factor. If C{accepted_reduction=None} no backtracking is used.  
         @type accepted_reduction: positive C{float} or C{None}  
         @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.  
         @type relaxation: C{float} or C{None}  
         """  
         tol=1.e-2  
         if tolerance_u==None: tolerance_u=tolerance  
         if not relaxation==None: self.relaxation=relaxation  
         if accepted_reduction ==None:  
               angle_limit=0.  
         elif accepted_reduction>=1.:  
               angle_limit=0.  
         else:  
               angle_limit=util.sqrt(1-accepted_reduction**2)  
         self.iter=0  
         u=u0  
         p=p0  
         #  
         #   initialize things:  
         #  
         converged=False  
         #  
         #  start loop:  
         #  
         #  initial search direction is g  
         #  
         while not converged :  
             if self.iter>iter_max:  
                 raise ArithmeticError("no convergence after %s steps."%self.iter)  
             f_new=self.solve_f(u,p,tol)  
             norm_f_new = util.Lsup(f_new)  
             u_new=u-f_new  
             g_new=self.solve_g(u_new,tol)  
             self.iter+=1  
             norm_g_new = util.sqrt(self.inner(g_new,g_new))  
             if norm_f_new==0. and norm_g_new==0.: return u, p  
             if self.iter>1 and not accepted_reduction==None:  
                #  
                #   did we manage to reduce the norm of G? I  
                #   if not we start a backtracking procedure  
                #  
                # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g  
                if norm_g_new > accepted_reduction * norm_g:  
                   sub_iter=0  
                   s=self.relaxation  
                   d=g  
                   g_last=g  
                   self.trace("    start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))  
                   while sub_iter < self.subiter_max and  norm_g_new > accepted_reduction * norm_g:  
                      dg= g_new-g_last  
                      norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)  
                      rad=self.inner(g_new,dg)/self.relaxation  
                      # print "   ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit  
                      # print "   ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g  
                      if abs(rad) < norm_dg*norm_g_new * angle_limit:  
                          if sub_iter>0: self.trace("    no further improvements expected from backtracking.")  
                          break  
                      r=self.relaxation  
                      self.relaxation= - rad/norm_dg**2  
                      s+=self.relaxation  
                      #####  
                      # a=g_new+self.relaxation*dg/r  
                      # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation  
                      #####  
                      g_last=g_new  
                      p+=self.relaxation*d  
                      f_new=self.solve_f(u,p,tol)  
                      u_new=u-f_new  
                      g_new=self.solve_g(u_new,tol)  
                      self.iter+=1  
                      norm_f_new = util.Lsup(f_new)  
                      norm_g_new = util.sqrt(self.inner(g_new,g_new))  
                      # print "   ",sub_iter," new g norm",norm_g_new  
                      self.trace("    %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))  
                      #  
                      #   can we expect reduction of g?  
                      #  
                      # u_last=u_new  
                      sub_iter+=1  
                   self.relaxation=s  
             #  
             #  check for convergence:  
             #  
             norm_u_new = util.Lsup(u_new)  
             p_new=p+self.relaxation*g_new  
             norm_p_new = util.sqrt(self.inner(p_new,p_new))  
             self.trace("%s th step: f = %s, f/u = %s, g = %s, g/p = %s, relaxation = %s."%(self.iter, norm_f_new ,norm_f_new/norm_u_new, norm_g_new, norm_g_new/norm_p_new, self.relaxation))  
   
             if self.iter>1:  
                dg2=g_new-g  
                df2=f_new-f  
                norm_dg2=util.sqrt(self.inner(dg2,dg2))  
                norm_df2=util.Lsup(df2)  
                # print norm_g_new, norm_g, norm_dg, norm_p, tolerance  
                tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new  
                tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new  
                if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:  
                    converged=True  
             f, norm_f, u, norm_u, g, norm_g, p, norm_p = f_new, norm_f_new, u_new, norm_u_new, g_new, norm_g_new, p_new, norm_p_new  
         self.trace("convergence after %s steps."%self.iter)  
         return u,p  

Legend:
Removed from v.1956  
changed lines
  Added in v.3981

  ViewVC Help
Powered by ViewVC 1.1.26