/[escript]/trunk/escript/py_src/pdetools.py
ViewVC logotype

Diff of /trunk/escript/py_src/pdetools.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

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

  ViewVC Help
Powered by ViewVC 1.1.26