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

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

  ViewVC Help
Powered by ViewVC 1.1.26