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

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

  ViewVC Help
Powered by ViewVC 1.1.26