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

  ViewVC Help
Powered by ViewVC 1.1.26