/[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 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC revision 4446 by caltinay, Tue Jun 11 04:00:15 2013 UTC
# Line 1  Line 1 
1    
2  ########################################################  ##############################################################################
3  #  #
4  # Copyright (c) 2003-2009 by University of Queensland  # Copyright (c) 2003-2013 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # http://www.uq.edu.au
 # http://www.uq.edu.au/esscc  
6  #  #
7  # Primary Business: Queensland, Australia  # Primary Business: Queensland, Australia
8  # Licensed under the Open Software License version 3.0  # Licensed under the Open Software License version 3.0
9  # http://www.opensource.org/licenses/osl-3.0.php  # http://www.opensource.org/licenses/osl-3.0.php
10  #  #
11  ########################################################  # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    # Development since 2012 by School of Earth Sciences
13    #
14    ##############################################################################
15    
16  __copyright__="""Copyright (c) 2003-2008 by University of Queensland  __copyright__="""Copyright (c) 2003-2013 by University of Queensland
17  Earth Systems Science Computational Center (ESSCC)  http://www.uq.edu.au
 http://www.uq.edu.au/esscc  
18  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
19  __license__="""Licensed under the Open Software License version 3.0  __license__="""Licensed under the Open Software License version 3.0
20  http://www.opensource.org/licenses/osl-3.0.php"""  http://www.opensource.org/licenses/osl-3.0.php"""
# Line 28  Currently includes: Line 29  Currently includes:
29      - TimeIntegrationManager - to handle extrapolation in time      - TimeIntegrationManager - to handle extrapolation in time
30      - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme      - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme
31    
32  @var __author__: name of author  :var __author__: name of author
33  @var __copyright__: copyrights  :var __copyright__: copyrights
34  @var __license__: licence agreement  :var __license__: licence agreement
35  @var __url__: url entry point on documentation  :var __url__: url entry point on documentation
36  @var __version__: version  :var __version__: version
37  @var __date__: date of the version  :var __date__: date of the version
38  """  """
39    
40  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
41    
42    
43  import escript  from . import escriptcpp as escore
44  import linearPDEs  from . import linearPDEs
45  import numpy  from . import util
 import util  
46  import math  import math
47    import numpy
48    
49  ##### Added by Artak  class TimeIntegrationManager(object):
 # from Numeric import zeros,Int,Float64  
 ###################################  
   
   
 class TimeIntegrationManager:  
50    """    """
51    A simple mechanism to manage time dependend values.    A simple mechanism to manage time dependend values.
52    
# Line 64  class TimeIntegrationManager: Line 60  class TimeIntegrationManager:
60           tm.checkin(dt,v)           tm.checkin(dt,v)
61           t+=dt           t+=dt
62    
63    @note: currently only p=1 is supported.    :note: currently only p=1 is supported.
64    """    """
65    def __init__(self,*inital_values,**kwargs):    def __init__(self,*inital_values,**kwargs):
66       """       """
67       Sets up the value manager where C{inital_values} are the initial values       Sets up the value manager where ``inital_values`` are the initial values
68       and p is the order used for extrapolation.       and p is the order used for extrapolation.
69       """       """
70       if kwargs.has_key("p"):       if "p" in kwargs:
71              self.__p=kwargs["p"]              self.__p=kwargs["p"]
72       else:       else:
73              self.__p=1              self.__p=1
74       if kwargs.has_key("time"):       if "time" in kwargs:
75              self.__t=kwargs["time"]              self.__t=kwargs["time"]
76       else:       else:
77              self.__t=0.              self.__t=0.
# Line 113  class TimeIntegrationManager: Line 109  class TimeIntegrationManager:
109    
110    def extrapolate(self,dt):    def extrapolate(self,dt):
111        """        """
112        Extrapolates to C{dt} forward in time.        Extrapolates to ``dt`` forward in time.
113        """        """
114        if self.__order==0:        if self.__order==0:
115           out=self.__v_mem[0]           out=self.__v_mem[0]
# Line 130  class TimeIntegrationManager: Line 126  class TimeIntegrationManager:
126           return out           return out
127    
128    
129  class Projector:  class Projector(object):
130    """    """
131    The Projector is a factory which projects a discontinuous function onto a    The Projector is a factory which projects a discontinuous function onto a
132    continuous function on a given domain.    continuous function on a given domain.
# Line 139  class Projector: Line 135  class Projector:
135      """      """
136      Creates a continuous function space projector for a domain.      Creates a continuous function space projector for a domain.
137    
138      @param domain: Domain of the projection.      :param domain: Domain of the projection.
139      @param reduce: Flag to reduce projection order      :param reduce: Flag to reduce projection order
140      @param fast: Flag to use a fast method based on matrix lumping      :param fast: Flag to use a fast method based on matrix lumping
141      """      """
142      self.__pde = linearPDEs.LinearPDE(domain)      self.__pde = linearPDEs.LinearPDE(domain)
143      if fast:      if fast:
# Line 151  class Projector: Line 147  class Projector:
147      self.__pde.setValue(D = 1.)      self.__pde.setValue(D = 1.)
148      return      return
149    def getSolverOptions(self):    def getSolverOptions(self):
150      """      """
151      Returns the solver options of the PDE solver.      Returns the solver options of the PDE solver.
152            
153      @rtype: L{linearPDEs.SolverOptions}      :rtype: `linearPDEs.SolverOptions`
154      """      """
155        return self.__pde.getSolverOptions()
156    
157      def getValue(self, input_data):
158        """
159        Projects ``input_data`` onto a continuous function.
160    
161        :param input_data: the data to be projected
162        """
163        return self(input_data)
164    
165    def __call__(self, input_data):    def __call__(self, input_data):
166      """      """
167      Projects C{input_data} onto a continuous function.      Projects ``input_data`` onto a continuous function.
168    
169      @param input_data: the data to be projected      :param input_data: the data to be projected
170      """      """
171      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())      out=escore.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
172      self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())      self.__pde.setValue(Y = escore.Data(), Y_reduced = escore.Data())
173      if input_data.getRank()==0:      if input_data.getRank()==0:
174          self.__pde.setValue(Y = input_data)          self.__pde.setValue(Y = input_data)
175          out=self.__pde.getSolution()          out=self.__pde.getSolution()
# Line 192  class Projector: Line 197  class Projector:
197                      out[i0,i1,i2,i3]=self.__pde.getSolution()                      out[i0,i1,i2,i3]=self.__pde.getSolution()
198      return out      return out
199    
200  class NoPDE:  class NoPDE(object):
201       """       """
202       Solves the following problem for u:       Solves the following problem for u:
203    
204       M{kronecker[i,j]*D[j]*u[j]=Y[i]}       *kronecker[i,j]*D[j]*u[j]=Y[i]*
205    
206       with constraint       with constraint
207    
208       M{u[j]=r[j]}  where M{q[j]>0}       *u[j]=r[j]*  where *q[j]>0*
209    
210       where M{D}, M{Y}, M{r} and M{q} are given functions of rank 1.       where *D*, *Y*, *r* and *q* are given functions of rank 1.
211    
212       In the case of scalars this takes the form       In the case of scalars this takes the form
213    
214       M{D*u=Y}       *D*u=Y*
215    
216       with constraint       with constraint
217    
218       M{u=r} where M{q>0}       *u=r* where *q>0*
219    
220       where M{D}, M{Y}, M{r} and M{q} are given scalar functions.       where *D*, *Y*, *r* and *q* are given scalar functions.
221    
222       The constraint overwrites any other condition.       The constraint overwrites any other condition.
223    
224       @note: This class is similar to the L{linearPDEs.LinearPDE} class with       :note: This class is similar to the `linearPDEs.LinearPDE` class with
225              A=B=C=X=0 but has the intention that all input parameters are given              A=B=C=X=0 but has the intention that all input parameters are given
226              in L{Solution} or L{ReducedSolution}.              in `Solution` or `ReducedSolution`.
227       """       """
228       # The whole thing is a bit strange and I blame Rob Woodcock (CSIRO) for       # The whole thing is a bit strange and I blame Rob Woodcock (CSIRO) for
229       # this.       # this.
# Line 226  class NoPDE: Line 231  class NoPDE:
231           """           """
232           Initializes the problem.           Initializes the problem.
233    
234           @param domain: domain of the PDE           :param domain: domain of the PDE
235           @type domain: L{Domain}           :type domain: `Domain`
236           @param D: coefficient of the solution           :param D: coefficient of the solution
237           @type D: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type D: ``float``, ``int``, ``numpy.ndarray``, `Data`
238           @param Y: right hand side           :param Y: right hand side
239           @type Y: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type Y: ``float``, ``int``, ``numpy.ndarray``, `Data`
240           @param q: location of constraints           :param q: location of constraints
241           @type q: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type q: ``float``, ``int``, ``numpy.ndarray``, `Data`
242           @param r: value of solution at locations of constraints           :param r: value of solution at locations of constraints
243           @type r: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type r: ``float``, ``int``, ``numpy.ndarray``, `Data`
244           """           """
245           self.__domain=domain           self.__domain=domain
246           self.__D=D           self.__D=D
# Line 243  class NoPDE: Line 248  class NoPDE:
248           self.__q=q           self.__q=q
249           self.__r=r           self.__r=r
250           self.__u=None           self.__u=None
251           self.__function_space=escript.Solution(self.__domain)           self.__function_space=escore.Solution(self.__domain)
252    
253       def setReducedOn(self):       def setReducedOn(self):
254           """           """
255           Sets the L{FunctionSpace} of the solution to L{ReducedSolution}.           Sets the `FunctionSpace` of the solution to `ReducedSolution`.
256           """           """
257           self.__function_space=escript.ReducedSolution(self.__domain)           self.__function_space=escore.ReducedSolution(self.__domain)
258           self.__u=None           self.__u=None
259    
260       def setReducedOff(self):       def setReducedOff(self):
261           """           """
262           Sets the L{FunctionSpace} of the solution to L{Solution}.           Sets the `FunctionSpace` of the solution to `Solution`.
263           """           """
264           self.__function_space=escript.Solution(self.__domain)           self.__function_space=escore.Solution(self.__domain)
265           self.__u=None           self.__u=None
266    
267       def setValue(self,D=None,Y=None,q=None,r=None):       def setValue(self,D=None,Y=None,q=None,r=None):
268           """           """
269           Assigns values to the parameters.           Assigns values to the parameters.
270    
271           @param D: coefficient of the solution           :param D: coefficient of the solution
272           @type D: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type D: ``float``, ``int``, ``numpy.ndarray``, `Data`
273           @param Y: right hand side           :param Y: right hand side
274           @type Y: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type Y: ``float``, ``int``, ``numpy.ndarray``, `Data`
275           @param q: location of constraints           :param q: location of constraints
276           @type q: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type q: ``float``, ``int``, ``numpy.ndarray``, `Data`
277           @param r: value of solution at locations of constraints           :param r: value of solution at locations of constraints
278           @type r: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type r: ``float``, ``int``, ``numpy.ndarray``, `Data`
279           """           """
280           if not D==None:           if not D==None:
281              self.__D=D              self.__D=D
# Line 289  class NoPDE: Line 294  class NoPDE:
294           """           """
295           Returns the solution.           Returns the solution.
296    
297           @return: the solution of the problem           :return: the solution of the problem
298           @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or           :rtype: `Data` object in the `FunctionSpace` `Solution` or
299                   L{ReducedSolution}                   `ReducedSolution`
300           """           """
301           if self.__u==None:           if self.__u==None:
302              if self.__D==None:              if self.__D==None:
303                 raise ValueError,"coefficient D is undefined"                 raise ValueError("coefficient D is undefined")
304              D=escript.Data(self.__D,self.__function_space)              D=escore.Data(self.__D,self.__function_space)
305              if D.getRank()>1:              if D.getRank()>1:
306                 raise ValueError,"coefficient D must have rank 0 or 1"                 raise ValueError("coefficient D must have rank 0 or 1")
307              if self.__Y==None:              if self.__Y==None:
308                 self.__u=escript.Data(0.,D.getShape(),self.__function_space)                 self.__u=escore.Data(0.,D.getShape(),self.__function_space)
309              else:              else:
310                 self.__u=util.quotient(self.__Y,D)                 self.__u=1./D*self.__Y
311              if not self.__q==None:              if not self.__q==None:
312                  q=util.wherePositive(escript.Data(self.__q,self.__function_space))                  q=util.wherePositive(escore.Data(self.__q,self.__function_space))
313                  self.__u*=(1.-q)                  self.__u*=(1.-q)
314                  if not self.__r==None: self.__u+=q*self.__r                  if not self.__r==None: self.__u+=q*self.__r
315           return self.__u           return self.__u
316    
317  class Locator:  class Locator(object):
318       """       """
319       Locator provides access to the values of data objects at a given spatial       Locator provides access to the values of data objects at a given spatial
320       coordinate x.       coordinate x.
# Line 324  class Locator: Line 329  class Locator:
329         or FunctionSpace for the sample point which is closest to the given         or FunctionSpace for the sample point which is closest to the given
330         point x.         point x.
331    
332         @param where: function space         :param where: function space
333         @type where: L{escript.FunctionSpace}         :type where: `escript.FunctionSpace`
334         @param x: location(s) of the Locator         :param x: location(s) of the Locator
335         @type x: C{numpy.ndarray} or C{list} of C{numpy.ndarray}         :type x: ``numpy.ndarray`` or ``list`` of ``numpy.ndarray``
336         """         """
337         if isinstance(where,escript.FunctionSpace):         if isinstance(where,escore.FunctionSpace):
338            self.__function_space=where            self.__function_space=where
339         else:         else:
340            self.__function_space=escript.ContinuousFunction(where)            self.__function_space=escore.ContinuousFunction(where)
341         iterative=False         iterative=False
342         if isinstance(x, list):         if isinstance(x, list):
343             if len(x)==0:             if len(x)==0:
344                raise "ValueError", "At least one point must be given."                raise ValueError("At least one point must be given.")
345             try:             try:
346               iter(x[0])               iter(x[0])
347               iterative=True               iterative=True
348             except TypeError:             except TypeError:
349               iterative=False               iterative=False
350           xxx=self.__function_space.getX()
351         if iterative:         if iterative:
352             self.__id=[]             self.__id=[]
353             for p in x:             for p in x:
354                self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())                self.__id.append(util.length(xxx-p[:self.__function_space.getDim()]).minGlobalDataPoint())
355         else:         else:
356             self.__id=util.length(self.__function_space.getX()-x[:self.__function_space.getDim()]).minGlobalDataPoint()             self.__id=util.length(xxx-x[:self.__function_space.getDim()]).minGlobalDataPoint()
357    
358       def __str__(self):       def __str__(self):
359         """         """
360         Returns the coordinates of the Locator as a string.         Returns the coordinates of the Locator as a string.
361         """         """
362         x=self.getX()         x=self.getX()
363         if instance(x,list):         if isinstance(x,list):
364            out="["            out="["
365            first=True            first=True
366            for xx in x:            for xx in x:
# Line 401  class Locator: Line 407  class Locator:
407    
408       def getValue(self,data):       def getValue(self,data):
409          """          """
410          Returns the value of C{data} at the Locator if C{data} is a L{Data}          Returns the value of ``data`` at the Locator if ``data`` is a `Data`
411          object otherwise the object is returned.          object otherwise the object is returned.
412          """          """
413          if isinstance(data,escript.Data):          if isinstance(data,escore.Data):
414             dat=util.interpolate(data,self.getFunctionSpace())             dat=util.interpolate(data,self.getFunctionSpace())
415             id=self.getId()             id=self.getId()
416             r=data.getRank()             r=data.getRank()
# Line 425  class Locator: Line 431  class Locator:
431                  return out                  return out
432          else:          else:
433             return data             return data
434              
435         def setValue(self, data, v):
436          """
437          Sets the value of the ``data`` at the Locator.
438          """
439          if isinstance(data, escore.Data):
440             if data.getFunctionSpace()!=self.getFunctionSpace():
441               raise TypeError("setValue: FunctionSpace of Locator and Data object must match.")
442             data.expand()  
443             id=self.getId()
444             if isinstance(id, list):
445              for i in id:
446               data._setTupleForGlobalDataPoint(i[1], i[0], v)
447             else:
448               data._setTupleForGlobalDataPoint(id[1], id[0], v)
449          else:
450               raise TypeError("setValue: Invalid argument type.")
451    
452    
453    def getInfLocator(arg):
454        """
455        Return a Locator for a point with the inf value over all arg.
456        """
457        if not isinstance(arg, escore.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, escore.Data):
470           raise TypeError("getInfLocator: Unknown argument type.")
471        a_inf=util.sup(arg)
472        loc=util.length(arg-a_inf).minGlobalDataPoint()     # This gives us the location but not coords
473        x=arg.getFunctionSpace().getX()
474        x_min=x.getTupleForGlobalDataPoint(*loc)
475        return Locator(arg.getFunctionSpace(),x_min)
476            
477    
478  class SolverSchemeException(Exception):  class SolverSchemeException(Exception):
479     """     """
# Line 467  def PCG(r, Aprod, x, Msolve, bilinearfor Line 516  def PCG(r, Aprod, x, Msolve, bilinearfor
516     """     """
517     Solver for     Solver for
518    
519     M{Ax=b}     *Ax=b*
520    
521     with a symmetric and positive definite operator A (more details required!).     with a symmetric and positive definite operator A (more details required!).
522     It uses the conjugate gradient method with preconditioner M providing an     It uses the conjugate gradient method with preconditioner M providing an
# Line 475  def PCG(r, Aprod, x, Msolve, bilinearfor Line 524  def PCG(r, Aprod, x, Msolve, bilinearfor
524    
525     The iteration is terminated if     The iteration is terminated if
526    
527     M{|r| <= atol+rtol*|r0|}     *|r| <= atol+rtol*|r0|*
528    
529     where M{r0} is the initial residual and M{|.|} is the energy norm. In fact     where *r0* is the initial residual and *|.|* is the energy norm. In fact
530    
531     M{|r| = sqrt( bilinearform(Msolve(r),r))}     *|r| = sqrt( bilinearform(Msolve(r),r))*
532    
533     For details on the preconditioned conjugate gradient method see the book:     For details on the preconditioned conjugate gradient method see the book:
534    
# Line 487  def PCG(r, Aprod, x, Msolve, bilinearfor Line 536  def PCG(r, Aprod, x, Msolve, bilinearfor
536     T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,     T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
537     C. Romine, and H. van der Vorst}.     C. Romine, and H. van der Vorst}.
538    
539     @param r: initial residual M{r=b-Ax}. C{r} is altered.     :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)     :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
541     @param x: an initial guess for the solution     :param x: an initial guess for the solution
542     @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)     :type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
543     @param Aprod: returns the value Ax     :param Aprod: returns the value Ax
544     @type Aprod: function C{Aprod(x)} where C{x} is of the same object like     :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
545                  argument C{x}. The returned object needs to be of the same type                  argument ``x``. The returned object needs to be of the same type
546                  like argument C{r}.                  like argument ``r``.
547     @param Msolve: solves Mx=r     :param Msolve: solves Mx=r
548     @type Msolve: function C{Msolve(r)} where C{r} is of the same type like     :type Msolve: function ``Msolve(r)`` where ``r`` is of the same type like
549                   argument C{r}. The returned object needs to be of the same                   argument ``r``. The returned object needs to be of the same
550                   type like argument C{x}.                   type like argument ``x``.
551     @param bilinearform: inner product C{<x,r>}     :param bilinearform: inner product ``<x,r>``
552     @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same     :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
553                         type like argument C{x} and C{r} is. The returned value                         type like argument ``x`` and ``r`` is. The returned value
554                         is a C{float}.                         is a ``float``.
555     @param atol: absolute tolerance     :param atol: absolute tolerance
556     @type atol: non-negative C{float}     :type atol: non-negative ``float``
557     @param rtol: relative tolerance     :param rtol: relative tolerance
558     @type rtol: non-negative C{float}     :type rtol: non-negative ``float``
559     @param iter_max: maximum number of iteration steps     :param iter_max: maximum number of iteration steps
560     @type iter_max: C{int}     :type iter_max: ``int``
561     @return: the solution approximation and the corresponding residual     :return: the solution approximation and the corresponding residual
562     @rtype: C{tuple}     :rtype: ``tuple``
563     @warning: C{r} and C{x} are altered.     :warning: ``r`` and ``x`` are altered.
564     """     """
565     iter=0     iter=0
566     rhat=Msolve(r)     rhat=Msolve(r)
567     d = rhat     d = rhat
568     rhat_dot_r = bilinearform(rhat, r)     rhat_dot_r = bilinearform(rhat, r)
569     if rhat_dot_r<0: raise NegativeNorm,"negative norm."     if rhat_dot_r<0: raise NegativeNorm("negative norm.")
570     norm_r0=math.sqrt(rhat_dot_r)     norm_r0=math.sqrt(rhat_dot_r)
571     atol2=atol+rtol*norm_r0     atol2=atol+rtol*norm_r0
572     if atol2<=0:     if atol2<=0:
573        raise ValueError,"Non-positive tolarance."        raise ValueError("Non-positive tolarance.")
574     atol2=max(atol2, 100. * util.EPSILON * norm_r0)     atol2=max(atol2, 100. * util.EPSILON * norm_r0)
575    
576     if verbose: print "PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_r0, atol2)     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:     while not math.sqrt(rhat_dot_r) <= atol2:
580         iter+=1         iter+=1
581         if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max         if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
582    
583         q=Aprod(d)         q=Aprod(d)
584         alpha = rhat_dot_r / bilinearform(d, q)         alpha = rhat_dot_r / bilinearform(d, q)
585         x += alpha * d         x += alpha * d
586         if isinstance(q,ArithmeticTuple):         if isinstance(q,ArithmeticTuple):
587         r += q * (-alpha)      # Doing it the other way calls the float64.__mul__ not AT.__rmul__            r += q * (-alpha)      # Doing it the other way calls the float64.__mul__ not AT.__rmul__
588         else:         else:
589             r += (-alpha) * q             r += (-alpha) * q
590         rhat=Msolve(r)         rhat=Msolve(r)
# Line 545  def PCG(r, Aprod, x, Msolve, bilinearfor Line 594  def PCG(r, Aprod, x, Msolve, bilinearfor
594         d=rhat         d=rhat
595    
596         rhat_dot_r = rhat_dot_r_new         rhat_dot_r = rhat_dot_r_new
597         if rhat_dot_r<0: raise NegativeNorm,"negative norm."         if rhat_dot_r<0: raise NegativeNorm("negative norm.")
598         if verbose: print "PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))         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     if verbose: print(("PCG: tolerance reached after %s steps."%iter))
600     return x,r,math.sqrt(rhat_dot_r)     return x,r,math.sqrt(rhat_dot_r)
601    
602  class Defect(object):  class Defect(object):
# Line 564  class Defect(object): Line 613  class Defect(object):
613          """          """
614          Returns the inner product of x0 and x1          Returns the inner product of x0 and x1
615    
616          @param x0: value for x0          :param x0: value for x0
617          @param x1: value for x1          :param x1: value for x1
618          @return: the inner product of x0 and x1          :return: the inner product of x0 and x1
619          @rtype: C{float}          :rtype: ``float``
620          """          """
621          return 0          return 0
622    
623      def norm(self,x):      def norm(self,x):
624          """          """
625          Returns the norm of argument C{x}.          Returns the norm of argument ``x``.
626    
627          @param x: a value          :param x: a value
628          @return: norm of argument x          :return: norm of argument x
629          @rtype: C{float}          :rtype: ``float``
630          @note: by default C{sqrt(self.bilinearform(x,x)} is returned.          :note: by default ``sqrt(self.bilinearform(x,x)`` is returned.
631          """          """
632          s=self.bilinearform(x,x)          s=self.bilinearform(x,x)
633          if s<0: raise NegativeNorm,"negative norm."          if s<0: raise NegativeNorm("negative norm.")
634          return math.sqrt(s)          return math.sqrt(s)
635    
636      def eval(self,x):      def eval(self,x):
637          """          """
638          Returns the value F of a given C{x}.          Returns the value F of a given ``x``.
639    
640          @param x: value for which the defect C{F} is evaluated          :param x: value for which the defect ``F`` is evaluated
641          @return: value of the defect at C{x}          :return: value of the defect at ``x``
642          """          """
643          return 0          return 0
644    
645      def __call__(self,x):      def __call__(self,x):
646          return self.eval(x)          return self.eval(x)
647    
648      def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):      def setDerivativeIncrementLength(self,inc=1000.*math.sqrt(util.EPSILON)):
649          """          """
650          Sets the relative length of the increment used to approximate the          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          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.          direction of v with x as a starting point.
653    
654          @param inc: relative increment length          :param inc: relative increment length
655          @type inc: positive C{float}          :type inc: positive ``float``
656          """          """
657          if inc<=0: raise ValueError,"positive increment required."          if inc<=0: raise ValueError("positive increment required.")
658          self.__inc=inc          self.__inc=inc
659    
660      def getDerivativeIncrementLength(self):      def getDerivativeIncrementLength(self):
661          """          """
662          Returns the relative increment length used to approximate the          Returns the relative increment length used to approximate the
663          derivative of the defect.          derivative of the defect.
664          @return: value of the defect at C{x}          :return: value of the defect at ``x``
665          @rtype: positive C{float}          :rtype: positive ``float``
666          """          """
667          return self.__inc          return self.__inc
668    
669      def derivative(self, F0, x0, v, v_is_normalised=True):      def derivative(self, F0, x0, v, v_is_normalised=True):
670          """          """
671          Returns the directional derivative at C{x0} in the direction of C{v}.          Returns the directional derivative at ``x0`` in the direction of ``v``.
672    
673          @param F0: value of this defect at x0          :param F0: value of this defect at x0
674          @param x0: value at which derivative is calculated          :param x0: value at which derivative is calculated
675          @param v: direction          :param v: direction
676          @param v_is_normalised: True to indicate that C{v} is nomalized          :param v_is_normalised: True to indicate that ``v`` is nomalized
677                                  (self.norm(v)=0)                                  (self.norm(v)=0)
678          @return: derivative of this defect at x0 in the direction of C{v}          :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          :note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is
680                 used but this method maybe overwritten to use exact evaluation.                 used but this method maybe overwritten to use exact evaluation.
681          """          """
682          normx=self.norm(x0)          normx=self.norm(x0)
# Line 645  class Defect(object): Line 694  class Defect(object):
694          return (F1-F0)/epsnew          return (F1-F0)/epsnew
695    
696  ######################################  ######################################
697  def NewtonGMRES(defect, x, iter_max=100, sub_iter_max=20, atol=0,rtol=1.e-4, sub_tol_max=0.5, gamma=0.9, verbose=False):  def NewtonGMRES(defect, x, iter_max=100, sub_iter_max=20, atol=0,rtol=1.e-4, subtol_max=0.5, gamma=0.9, verbose=False):
698     """     """
699     Solves a non-linear problem M{F(x)=0} for unknown M{x} using the stopping     Solves a non-linear problem *F(x)=0* for unknown *x* using the stopping
700     criterion:     criterion:
701    
702     M{norm(F(x) <= atol + rtol * norm(F(x0)}     *norm(F(x) <= atol + rtol * norm(F(x0)*
703    
704     where M{x0} is the initial guess.     where *x0* is the initial guess.
705    
706     @param defect: object defining the function M{F}. C{defect.norm} defines the     :param defect: object defining the function *F*. ``defect.norm`` defines the
707                    M{norm} used in the stopping criterion.                    *norm* used in the stopping criterion.
708     @type defect: L{Defect}     :type defect: `Defect`
709     @param x: initial guess for the solution, C{x} is altered.     :param x: initial guess for the solution, ``x`` is altered.
710     @type x: any object type allowing basic operations such as     :type x: any object type allowing basic operations such as
711              C{numpy.ndarray}, L{Data}              ``numpy.ndarray``, `Data`
712     @param iter_max: maximum number of iteration steps     :param iter_max: maximum number of iteration steps
713     @type iter_max: positive C{int}     :type iter_max: positive ``int``
714     @param sub_iter_max: maximum number of inner iteration steps     :param sub_iter_max: maximum number of inner iteration steps
715     @type sub_iter_max: positive C{int}     :type sub_iter_max: positive ``int``
716     @param atol: absolute tolerance for the solution     :param atol: absolute tolerance for the solution
717     @type atol: positive C{float}     :type atol: positive ``float``
718     @param rtol: relative tolerance for the solution     :param rtol: relative tolerance for the solution
719     @type rtol: positive C{float}     :type rtol: positive ``float``
720     @param gamma: tolerance safety factor for inner iteration     :param gamma: tolerance safety factor for inner iteration
721     @type gamma: positive C{float}, less than 1     :type gamma: positive ``float``, less than 1
722     @param sub_tol_max: upper bound for inner tolerance     :param subtol_max: upper bound for inner tolerance
723     @type sub_tol_max: positive C{float}, less than 1     :type subtol_max: positive ``float``, less than 1
724     @return: an approximation of the solution with the desired accuracy     :return: an approximation of the solution with the desired accuracy
725     @rtype: same type as the initial guess     :rtype: same type as the initial guess
726     """     """
727     lmaxit=iter_max     lmaxit=iter_max
728     if atol<0: raise ValueError,"atol needs to be non-negative."     if atol<0: raise ValueError("atol needs to be non-negative.")
729     if rtol<0: raise ValueError,"rtol needs to be non-negative."     if rtol<0: raise ValueError("rtol needs to be non-negative.")
730     if rtol+atol<=0: raise ValueError,"rtol or atol needs to be non-negative."     if rtol+atol<=0: raise ValueError("rtol or atol needs to be non-negative.")
731     if gamma<=0 or gamma>=1: raise ValueError,"tolerance safety factor for inner iteration (gamma =%s) needs to be positive and less than 1."%gamma     if gamma<=0 or gamma>=1: raise ValueError("tolerance safety factor for inner iteration (gamma =%s) needs to be positive and less than 1."%gamma)
732     if sub_tol_max<=0 or sub_tol_max>=1: raise ValueError,"upper bound for inner tolerance for inner iteration (sub_tol_max =%s) needs to be positive and less than 1."%sub_tol_max     if subtol_max<=0 or subtol_max>=1: raise ValueError("upper bound for inner tolerance for inner iteration (subtol_max =%s) needs to be positive and less than 1."%subtol_max)
733    
734     F=defect(x)     F=defect(x)
735     fnrm=defect.norm(F)     fnrm=defect.norm(F)
736     stop_tol=atol + rtol*fnrm     stop_tol=atol + rtol*fnrm
737     sub_tol=sub_tol_max     subtol=subtol_max
738     if verbose: print "NewtonGMRES: initial residual = %e."%fnrm     if verbose: print(("NewtonGMRES: initial residual = %e."%fnrm))
739     if verbose: print "             tolerance = %e."%sub_tol     if verbose: print(("             tolerance = %e."%subtol))
740     iter=1     iter=1
741     #     #
742     # main iteration loop     # main iteration loop
743     #     #
744     while not fnrm<=stop_tol:     while not fnrm<=stop_tol:
745              if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max              if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
746                #
747                #   adjust subtol_
748              #              #
         #   adjust sub_tol_  
         #  
749              if iter > 1:              if iter > 1:
750             rat=fnrm/fnrmo                 rat=fnrm/fnrmo
751                 sub_tol_old=sub_tol                 subtol_old=subtol
752             sub_tol=gamma*rat**2                 subtol=gamma*rat**2
753             if gamma*sub_tol_old**2 > .1: sub_tol=max(sub_tol,gamma*sub_tol_old**2)                 if gamma*subtol_old**2 > .1: subtol=max(subtol,gamma*subtol_old**2)
754             sub_tol=max(min(sub_tol,sub_tol_max), .5*stop_tol/fnrm)                 subtol=max(min(subtol,subtol_max), .5*stop_tol/fnrm)
755          #              #
756          # calculate newton increment xc              # calculate newton increment xc
757              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown
758              #     if iter_restart -1 is returned as sub_iter              #     if iter_restart -1 is returned as sub_iter
759              #     if  atol is reached sub_iter returns the numer of steps performed to get there              #     if  atol is reached sub_iter returns the numer of steps performed to get there
760              #              #
761              #              #
762              if verbose: print "             subiteration (GMRES) is called with relative tolerance %e."%sub_tol              if verbose: print(("             subiteration (GMRES) is called with relative tolerance %e."%subtol))
763              try:              try:
764                 xc, sub_iter=__FDGMRES(F, defect, x, sub_tol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)                 xc, sub_iter=__FDGMRES(F, defect, x, subtol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)
765              except MaxIterReached:              except MaxIterReached:
766                 raise MaxIterReached,"maximum number of %s steps reached."%iter_max                 raise MaxIterReached("maximum number of %s steps reached."%iter_max)
767              if sub_iter<0:              if sub_iter<0:
768                 iter+=sub_iter_max                 iter+=sub_iter_max
769              else:              else:
770                 iter+=sub_iter                 iter+=sub_iter
771              # ====              # ====
772          x+=xc              x+=xc
773              F=defect(x)              F=defect(x)
774          iter+=1              iter+=1
775              fnrmo, fnrm=fnrm, defect.norm(F)              fnrmo, fnrm=fnrm, defect.norm(F)
776              if verbose: print "             step %s: residual %e."%(iter,fnrm)              if verbose: print(("             step %s: residual %e."%(iter,fnrm)))
777     if verbose: print "NewtonGMRES: completed after %s steps."%iter     if verbose: print(("NewtonGMRES: completed after %s steps."%iter))
778     return x     return x
779    
780  def __givapp(c,s,vin):  def __givapp(c,s,vin):
781      """      """
782      Applies a sequence of Givens rotations (c,s) recursively to the vector      Applies a sequence of Givens rotations (c,s) recursively to the vector
783      C{vin}      ``vin``
784    
785      @warning: C{vin} is altered.      :warning: ``vin`` is altered.
786      """      """
787      vrot=vin      vrot=vin
788      if isinstance(c,float):      if isinstance(c,float):
# Line 741  def __givapp(c,s,vin): Line 790  def __givapp(c,s,vin):
790      else:      else:
791          for i in range(len(c)):          for i in range(len(c)):
792              w1=c[i]*vrot[i]-s[i]*vrot[i+1]              w1=c[i]*vrot[i]-s[i]*vrot[i+1]
793          w2=s[i]*vrot[i]+c[i]*vrot[i+1]              w2=s[i]*vrot[i]+c[i]*vrot[i+1]
794              vrot[i]=w1              vrot[i]=w1
795              vrot[i+1]=w2              vrot[i+1]=w2
796      return vrot      return vrot
# Line 761  def __FDGMRES(F0, defect, x0, atol, iter Line 810  def __FDGMRES(F0, defect, x0, atol, iter
810     iter=0     iter=0
811     while rho > atol and iter<iter_restart-1:     while rho > atol and iter<iter_restart-1:
812          if iter  >= iter_max:          if iter  >= iter_max:
813              raise MaxIterReached,"maximum number of %s steps reached."%iter_max              raise MaxIterReached("maximum number of %s steps reached."%iter_max)
814    
815          p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)          p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
816          v.append(p)          v.append(p)
# Line 823  def __FDGMRES(F0, defect, x0, atol, iter Line 872  def __FDGMRES(F0, defect, x0, atol, iter
872            i=i-1            i=i-1
873       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
874       for i in range(iter-1):       for i in range(iter-1):
875      xhat += v[i]*y[i]         xhat += v[i]*y[i]
876     else :     else :
877        xhat=v[0] * 0        xhat=v[0] * 0
878    
# Line 838  def GMRES(r, Aprod, x, bilinearform, ato Line 887  def GMRES(r, Aprod, x, bilinearform, ato
887     """     """
888     Solver for     Solver for
889    
890     M{Ax=b}     *Ax=b*
891    
892     with a general operator A (more details required!).     with a general operator A (more details required!).
893     It uses the generalized minimum residual method (GMRES).     It uses the generalized minimum residual method (GMRES).
894    
895     The iteration is terminated if     The iteration is terminated if
896    
897     M{|r| <= atol+rtol*|r0|}     *|r| <= atol+rtol*|r0|*
898    
899     where M{r0} is the initial residual and M{|.|} is the energy norm. In fact     where *r0* is the initial residual and *|.|* is the energy norm. In fact
900    
901     M{|r| = sqrt( bilinearform(r,r))}     *|r| = sqrt( bilinearform(r,r))*
902    
903     @param r: initial residual M{r=b-Ax}. C{r} is altered.     :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)     :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
905     @param x: an initial guess for the solution     :param x: an initial guess for the solution
906     @type x: same like C{r}     :type x: same like ``r``
907     @param Aprod: returns the value Ax     :param Aprod: returns the value Ax
908     @type Aprod: function C{Aprod(x)} where C{x} is of the same object like     :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
909                  argument C{x}. The returned object needs to be of the same                  argument ``x``. The returned object needs to be of the same
910                  type like argument C{r}.                  type like argument ``r``.
911     @param bilinearform: inner product C{<x,r>}     :param bilinearform: inner product ``<x,r>``
912     @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same     :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
913                         type like argument C{x} and C{r}. The returned value is                         type like argument ``x`` and ``r``. The returned value is
914                         a C{float}.                         a ``float``.
915     @param atol: absolute tolerance     :param atol: absolute tolerance
916     @type atol: non-negative C{float}     :type atol: non-negative ``float``
917     @param rtol: relative tolerance     :param rtol: relative tolerance
918     @type rtol: non-negative C{float}     :type rtol: non-negative ``float``
919     @param iter_max: maximum number of iteration steps     :param iter_max: maximum number of iteration steps
920     @type iter_max: C{int}     :type iter_max: ``int``
921     @param iter_restart: in order to save memory the orthogonalization process     :param iter_restart: in order to save memory the orthogonalization process
922                          is terminated after C{iter_restart} steps and the                          is terminated after ``iter_restart`` steps and the
923                          iteration is restarted.                          iteration is restarted.
924     @type iter_restart: C{int}     :type iter_restart: ``int``
925     @return: the solution approximation and the corresponding residual     :return: the solution approximation and the corresponding residual
926     @rtype: C{tuple}     :rtype: ``tuple``
927     @warning: C{r} and C{x} are altered.     :warning: ``r`` and ``x`` are altered.
928     """     """
929     m=iter_restart     m=iter_restart
930     restarted=False     restarted=False
931     iter=0     iter=0
932     if rtol>0:     if rtol>0:
933        r_dot_r = bilinearform(r, r)        r_dot_r = bilinearform(r, r)
934        if r_dot_r<0: raise NegativeNorm,"negative norm."        if r_dot_r<0: raise NegativeNorm("negative norm.")
935        atol2=atol+rtol*math.sqrt(r_dot_r)        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)        if verbose: print(("GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)))
937     else:     else:
938        atol2=atol        atol2=atol
939        if verbose: print "GMRES: absolute tolerance = %e"%atol2        if verbose: print(("GMRES: absolute tolerance = %e"%atol2))
940     if atol2<=0:     if atol2<=0:
941        raise ValueError,"Non-positive tolarance."        raise ValueError("Non-positive tolarance.")
942    
943     while True:     while True:
944        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max        if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached"%iter_max)
945        if restarted:        if restarted:
946           r2 = r-Aprod(x-x2)           r2 = r-Aprod(x-x2)
947        else:        else:
# Line 901  def GMRES(r, Aprod, x, bilinearform, ato Line 950  def GMRES(r, Aprod, x, bilinearform, ato
950        x,stopped=_GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose,P_R=P_R)        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        iter+=iter_restart
952        if stopped: break        if stopped: break
953        if verbose: print "GMRES: restart."        if verbose: print("GMRES: restart.")
954        restarted=True        restarted=True
955     if verbose: print "GMRES: tolerance has been reached."     if verbose: print("GMRES: tolerance has been reached.")
956     return x     return x
957    
958  def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False, P_R=None):  def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False, P_R=None):
# Line 916  def _GMRESm(r, Aprod, x, bilinearform, a Line 965  def _GMRESm(r, Aprod, x, bilinearform, a
965     v=[]     v=[]
966    
967     r_dot_r = bilinearform(r, r)     r_dot_r = bilinearform(r, r)
968     if r_dot_r<0: raise NegativeNorm,"negative norm."     if r_dot_r<0: raise NegativeNorm("negative norm.")
969     rho=math.sqrt(r_dot_r)     rho=math.sqrt(r_dot_r)
970    
971     v.append(r/rho)     v.append(r/rho)
972     g[0]=rho     g[0]=rho
973    
974     if verbose: print "GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)     if verbose: print(("GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)))
975     while not (rho<=atol or iter==iter_restart):     while not (rho<=atol or iter==iter_restart):
976    
977      if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max          if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
978    
979          if P_R!=None:          if P_R!=None:
980              p=Aprod(P_R(v[iter]))              p=Aprod(P_R(v[iter]))
981          else:          else:
982          p=Aprod(v[iter])              p=Aprod(v[iter])
983      v.append(p)          v.append(p)
984    
985      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))          v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
986    
987  # Modified Gram-Schmidt  # Modified Gram-Schmidt
988      for j in range(iter+1):          for j in range(iter+1):
989        h[j,iter]=bilinearform(v[j],v[iter+1])            h[j,iter]=bilinearform(v[j],v[iter+1])
990        v[iter+1]-=h[j,iter]*v[j]            v[iter+1]-=h[j,iter]*v[j]
991    
992      h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))          h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
993      v_norm2=h[iter+1,iter]          v_norm2=h[iter+1,iter]
994    
995  # Reorthogonalize if needed  # Reorthogonalize if needed
996      if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)          if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
997       for j in range(iter+1):           for j in range(iter+1):
998          hr=bilinearform(v[j],v[iter+1])              hr=bilinearform(v[j],v[iter+1])
999              h[j,iter]=h[j,iter]+hr              h[j,iter]=h[j,iter]+hr
1000              v[iter+1] -= hr*v[j]              v[iter+1] -= hr*v[j]
1001    
1002       v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))           v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
1003       h[iter+1,iter]=v_norm2           h[iter+1,iter]=v_norm2
1004    
1005  #   watch out for happy breakdown  #   watch out for happy breakdown
1006          if not v_norm2 == 0:          if not v_norm2 == 0:
1007           v[iter+1]=v[iter+1]/h[iter+1,iter]           v[iter+1]=v[iter+1]/h[iter+1,iter]
1008    
1009  #   Form and store the information for the new Givens rotation  #   Form and store the information for the new Givens rotation
1010      if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])          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])          mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
1012    
1013      if mu!=0 :          if mu!=0 :
1014          c[iter]=h[iter,iter]/mu                  c[iter]=h[iter,iter]/mu
1015          s[iter]=-h[iter+1,iter]/mu                  s[iter]=-h[iter+1,iter]/mu
1016          h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]                  h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
1017          h[iter+1,iter]=0.0                  h[iter+1,iter]=0.0
1018                  gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])                  gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
1019                  g[iter]=gg[0]                  g[iter]=gg[0]
1020                  g[iter+1]=gg[1]                  g[iter+1]=gg[1]
1021  # Update the residual norm  # Update the residual norm
1022    
1023          rho=abs(g[iter+1])          rho=abs(g[iter+1])
1024          if verbose: print "GMRES: iteration step %s: residual %e"%(iter,rho)          if verbose: print(("GMRES: iteration step %s: residual %e"%(iter,rho)))
1025      iter+=1          iter+=1
1026    
1027  # At this point either iter > iter_max or rho < tol.  # At this point either iter > iter_max or rho < tol.
1028  # It's time to compute x and leave.  # It's time to compute x and leave.
1029    
1030     if verbose: print "GMRES: iteration stopped after %s step."%iter     if verbose: print(("GMRES: iteration stopped after %s step."%iter))
1031     if iter > 0 :     if iter > 0 :
1032       y=numpy.zeros(iter,numpy.float64)       y=numpy.zeros(iter,numpy.float64)
1033       y[iter-1] = g[iter-1] / h[iter-1,iter-1]       y[iter-1] = g[iter-1] / h[iter-1,iter-1]
# Line 989  def _GMRESm(r, Aprod, x, bilinearform, a Line 1038  def _GMRESm(r, Aprod, x, bilinearform, a
1038            i=i-1            i=i-1
1039       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
1040       for i in range(iter-1):       for i in range(iter-1):
1041      xhat += v[i]*y[i]         xhat += v[i]*y[i]
1042     else:     else:
1043       xhat=v[0] * 0       xhat=v[0] * 0
1044     if P_R!=None:     if P_R!=None:
# Line 1007  def MINRES(r, Aprod, x, Msolve, bilinear Line 1056  def MINRES(r, Aprod, x, Msolve, bilinear
1056      """      """
1057      Solver for      Solver for
1058    
1059      M{Ax=b}      *Ax=b*
1060    
1061      with a symmetric and positive definite operator A (more details required!).      with a symmetric and positive definite operator A (more details required!).
1062      It uses the minimum residual method (MINRES) with preconditioner M      It uses the minimum residual method (MINRES) with preconditioner M
# Line 1015  def MINRES(r, Aprod, x, Msolve, bilinear Line 1064  def MINRES(r, Aprod, x, Msolve, bilinear
1064    
1065      The iteration is terminated if      The iteration is terminated if
1066    
1067      M{|r| <= atol+rtol*|r0|}      *|r| <= atol+rtol*|r0|*
1068    
1069      where M{r0} is the initial residual and M{|.|} is the energy norm. In fact      where *r0* is the initial residual and *|.|* is the energy norm. In fact
1070    
1071      M{|r| = sqrt( bilinearform(Msolve(r),r))}      *|r| = sqrt( bilinearform(Msolve(r),r))*
1072    
1073      For details on the preconditioned conjugate gradient method see the book:      For details on the preconditioned conjugate gradient method see the book:
1074    
# Line 1027  def MINRES(r, Aprod, x, Msolve, bilinear Line 1076  def MINRES(r, Aprod, x, Msolve, bilinear
1076      T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,      T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
1077      C. Romine, and H. van der Vorst}.      C. Romine, and H. van der Vorst}.
1078    
1079      @param r: initial residual M{r=b-Ax}. C{r} is altered.      :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)      :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1081      @param x: an initial guess for the solution      :param x: an initial guess for the solution
1082      @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)      :type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1083      @param Aprod: returns the value Ax      :param Aprod: returns the value Ax
1084      @type Aprod: function C{Aprod(x)} where C{x} is of the same object like      :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
1085                   argument C{x}. The returned object needs to be of the same                   argument ``x``. The returned object needs to be of the same
1086                   type like argument C{r}.                   type like argument ``r``.
1087      @param Msolve: solves Mx=r      :param Msolve: solves Mx=r
1088      @type Msolve: function C{Msolve(r)} where C{r} is of the same type like      :type Msolve: function ``Msolve(r)`` where ``r`` is of the same type like
1089                    argument C{r}. The returned object needs to be of the same                    argument ``r``. The returned object needs to be of the same
1090                    type like argument C{x}.                    type like argument ``x``.
1091      @param bilinearform: inner product C{<x,r>}      :param bilinearform: inner product ``<x,r>``
1092      @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same      :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
1093                          type like argument C{x} and C{r} is. The returned value                          type like argument ``x`` and ``r`` is. The returned value
1094                          is a C{float}.                          is a ``float``.
1095      @param atol: absolute tolerance      :param atol: absolute tolerance
1096      @type atol: non-negative C{float}      :type atol: non-negative ``float``
1097      @param rtol: relative tolerance      :param rtol: relative tolerance
1098      @type rtol: non-negative C{float}      :type rtol: non-negative ``float``
1099      @param iter_max: maximum number of iteration steps      :param iter_max: maximum number of iteration steps
1100      @type iter_max: C{int}      :type iter_max: ``int``
1101      @return: the solution approximation and the corresponding residual      :return: the solution approximation and the corresponding residual
1102      @rtype: C{tuple}      :rtype: ``tuple``
1103      @warning: C{r} and C{x} are altered.      :warning: ``r`` and ``x`` are altered.
1104      """      """
1105      #------------------------------------------------------------------      #------------------------------------------------------------------
1106      # Set up y and v for the first Lanczos vector v1.      # Set up y and v for the first Lanczos vector v1.
# Line 1062  def MINRES(r, Aprod, x, Msolve, bilinear Line 1111  def MINRES(r, Aprod, x, Msolve, bilinear
1111      y = Msolve(r)      y = Msolve(r)
1112      beta1 = bilinearform(y,r)      beta1 = bilinearform(y,r)
1113    
1114      if beta1< 0: raise NegativeNorm,"negative norm."      if beta1< 0: raise NegativeNorm("negative norm.")
1115    
1116      #  If r = 0 exactly, stop with x      #  If r = 0 exactly, stop with x
1117      if beta1==0: return x      if beta1==0: return x
# Line 1097  def MINRES(r, Aprod, x, Msolve, bilinear Line 1146  def MINRES(r, Aprod, x, Msolve, bilinear
1146      # --------------------------------------------------------------------      # --------------------------------------------------------------------
1147      while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL      while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL
1148    
1149      if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max          if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
1150          iter    = iter  +  1          iter    = iter  +  1
1151    
1152          #-----------------------------------------------------------------          #-----------------------------------------------------------------
# Line 1127  def MINRES(r, Aprod, x, Msolve, bilinear Line 1176  def MINRES(r, Aprod, x, Msolve, bilinear
1176          y = Msolve(r2)          y = Msolve(r2)
1177          oldb   = beta                           # oldb = betak          oldb   = beta                           # oldb = betak
1178          beta   = bilinearform(y,r2)             # beta = betak+1^2          beta   = bilinearform(y,r2)             # beta = betak+1^2
1179          if beta < 0: raise NegativeNorm,"negative norm."          if beta < 0: raise NegativeNorm("negative norm.")
1180    
1181          beta   = math.sqrt( beta )          beta   = math.sqrt( beta )
1182          tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta          tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
# Line 1185  def TFQMR(r, Aprod, x, bilinearform, ato Line 1234  def TFQMR(r, Aprod, x, bilinearform, ato
1234    """    """
1235    Solver for    Solver for
1236    
1237    M{Ax=b}    *Ax=b*
1238    
1239    with a general operator A (more details required!).    with a general operator A (more details required!).
1240    It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).    It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).
1241    
1242    The iteration is terminated if    The iteration is terminated if
1243    
1244    M{|r| <= atol+rtol*|r0|}    *|r| <= atol+rtol*|r0|*
1245    
1246    where M{r0} is the initial residual and M{|.|} is the energy norm. In fact    where *r0* is the initial residual and *|.|* is the energy norm. In fact
1247    
1248    M{|r| = sqrt( bilinearform(r,r))}    *|r| = sqrt( bilinearform(r,r))*
1249    
1250    @param r: initial residual M{r=b-Ax}. C{r} is altered.    :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)    :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1252    @param x: an initial guess for the solution    :param x: an initial guess for the solution
1253    @type x: same like C{r}    :type x: same like ``r``
1254    @param Aprod: returns the value Ax    :param Aprod: returns the value Ax
1255    @type Aprod: function C{Aprod(x)} where C{x} is of the same object like    :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
1256                 argument C{x}. The returned object needs to be of the same type                 argument ``x``. The returned object needs to be of the same type
1257                 like argument C{r}.                 like argument ``r``.
1258    @param bilinearform: inner product C{<x,r>}    :param bilinearform: inner product ``<x,r>``
1259    @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same    :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
1260                        type like argument C{x} and C{r}. The returned value is                        type like argument ``x`` and ``r``. The returned value is
1261                        a C{float}.                        a ``float``.
1262    @param atol: absolute tolerance    :param atol: absolute tolerance
1263    @type atol: non-negative C{float}    :type atol: non-negative ``float``
1264    @param rtol: relative tolerance    :param rtol: relative tolerance
1265    @type rtol: non-negative C{float}    :type rtol: non-negative ``float``
1266    @param iter_max: maximum number of iteration steps    :param iter_max: maximum number of iteration steps
1267    @type iter_max: C{int}    :type iter_max: ``int``
1268    @rtype: C{tuple}    :rtype: ``tuple``
1269    @warning: C{r} and C{x} are altered.    :warning: ``r`` and ``x`` are altered.
1270    """    """
1271    u1=0    u1=0
1272    u2=0    u2=0
# Line 1234  def TFQMR(r, Aprod, x, bilinearform, ato Line 1283  def TFQMR(r, Aprod, x, bilinearform, ato
1283    theta = 0.0;    theta = 0.0;
1284    eta = 0.0;    eta = 0.0;
1285    rho=bilinearform(r,r)    rho=bilinearform(r,r)
1286    if rho < 0: raise NegativeNorm,"negative norm."    if rho < 0: raise NegativeNorm("negative norm.")
1287    tau = math.sqrt(rho)    tau = math.sqrt(rho)
1288    norm_r0=tau    norm_r0=tau
1289    while tau>atol+rtol*norm_r0:    while tau>atol+rtol*norm_r0:
1290      if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max      if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
1291    
1292      sigma = bilinearform(r,v)      sigma = bilinearform(r,v)
1293      if sigma == 0.0: raise IterationBreakDown,'TFQMR breakdown, sigma=0'      if sigma == 0.0: raise IterationBreakDown('TFQMR breakdown, sigma=0')
1294    
1295      alpha = rho / sigma      alpha = rho / sigma
1296    
# Line 1269  def TFQMR(r, Aprod, x, bilinearform, ato Line 1318  def TFQMR(r, Aprod, x, bilinearform, ato
1318  #  #
1319  #  Try to terminate the iteration at each pass through the loop  #  Try to terminate the iteration at each pass through the loop
1320  #  #
1321      if rho == 0.0: raise IterationBreakDown,'TFQMR breakdown, rho=0'      if rho == 0.0: raise IterationBreakDown('TFQMR breakdown, rho=0')
1322    
1323      rhon = bilinearform(r,w)      rhon = bilinearform(r,w)
1324      beta = rhon / rho;      beta = rhon / rho;
# Line 1287  def TFQMR(r, Aprod, x, bilinearform, ato Line 1336  def TFQMR(r, Aprod, x, bilinearform, ato
1336    
1337  class ArithmeticTuple(object):  class ArithmeticTuple(object):
1338     """     """
1339     Tuple supporting inplace update x+=y and scaling x=a*y where C{x,y} is an     Tuple supporting inplace update x+=y and scaling x=a*y where ``x,y`` is an
1340     ArithmeticTuple and C{a} is a float.     ArithmeticTuple and ``a`` is a float.
1341    
1342     Example of usage::     Example of usage::
1343    
1344         from esys.escript import Data         from esys.escript import Data
1345         from numpy import array         from numpy import array
1346         a=Data(...)         a=eData(...)
1347         b=array([1.,4.])         b=array([1.,4.])
1348         x=ArithmeticTuple(a,b)         x=ArithmeticTuple(a,b)
1349         y=5.*x         y=5.*x
# Line 1302  class ArithmeticTuple(object): Line 1351  class ArithmeticTuple(object):
1351     """     """
1352     def __init__(self,*args):     def __init__(self,*args):
1353         """         """
1354         Initializes object with elements C{args}.         Initializes object with elements ``args``.
1355    
1356         @param args: tuple of objects that support inplace add (x+=y) and         :param args: tuple of objects that support inplace add (x+=y) and
1357                      scaling (x=a*y)                      scaling (x=a*y)
1358         """         """
1359         self.__items=list(args)         self.__items=list(args)
# Line 1313  class ArithmeticTuple(object): Line 1362  class ArithmeticTuple(object):
1362         """         """
1363         Returns the number of items.         Returns the number of items.
1364    
1365         @return: number of items         :return: number of items
1366         @rtype: C{int}         :rtype: ``int``
1367         """         """
1368         return len(self.__items)         return len(self.__items)
1369    
# Line 1322  class ArithmeticTuple(object): Line 1371  class ArithmeticTuple(object):
1371         """         """
1372         Returns item at specified position.         Returns item at specified position.
1373    
1374         @param index: index of item to be returned         :param index: index of item to be returned
1375         @type index: C{int}         :type index: ``int``
1376         @return: item with index C{index}         :return: item with index ``index``
1377         """         """
1378         return self.__items.__getitem__(index)         return self.__items.__getitem__(index)
1379    
1380     def __mul__(self,other):     def __mul__(self,other):
1381         """         """
1382         Scales by C{other} from the right.         Scales by ``other`` from the right.
1383    
1384         @param other: scaling factor         :param other: scaling factor
1385         @type other: C{float}         :type other: ``float``
1386         @return: itemwise self*other         :return: itemwise self*other
1387         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1388         """         """
1389         out=[]         out=[]
1390         try:         try:
1391             l=len(other)             l=len(other)
1392             if l!=len(self):             if l!=len(self):
1393                 raise ValueError,"length of arguments don't match."                 raise ValueError("length of arguments don't match.")
1394             for i in range(l): out.append(self[i]*other[i])             for i in range(l):
1395                    if self.__isEmpty(self[i]) or self.__isEmpty(other[i]):
1396                        out.append(escore.Data())
1397                    else:
1398                        out.append(self[i]*other[i])
1399         except TypeError:         except TypeError:
1400             for i in range(len(self)): out.append(self[i]*other)             for i in range(len(self)):  
1401                    if self.__isEmpty(self[i]) or self.__isEmpty(other):
1402                        out.append(escore.Data())
1403                    else:
1404                        out.append(self[i]*other)
1405         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1406    
1407     def __rmul__(self,other):     def __rmul__(self,other):
1408         """        """
1409         Scales by C{other} from the left.        Scales by ``other`` from the left.
1410    
1411         @param other: scaling factor        :param other: scaling factor
1412         @type other: C{float}        :type other: ``float``
1413         @return: itemwise other*self        :return: itemwise other*self
1414         @rtype: L{ArithmeticTuple}        :rtype: `ArithmeticTuple`
1415         """        """
1416         out=[]        out=[]
1417         try:        try:
1418             l=len(other)            l=len(other)
1419             if l!=len(self):            if l!=len(self):
1420                 raise ValueError,"length of arguments don't match."                raise ValueError("length of arguments don't match.")
1421             for i in range(l): out.append(other[i]*self[i])            for i in range(l):
1422         except TypeError:                  if self.__isEmpty(self[i]) or self.__isEmpty(other[i]):
1423             for i in range(len(self)): out.append(other*self[i])                      out.append(escore.Data())
1424         return ArithmeticTuple(*tuple(out))                  else:
1425                        out.append(other[i]*self[i])
1426          except TypeError:
1427              for i in range(len(self)):  
1428                    if self.__isEmpty(self[i]) or self.__isEmpty(other):
1429                        out.append(escore.Data())
1430                    else:
1431                        out.append(other*self[i])
1432          return ArithmeticTuple(*tuple(out))
1433    
1434     def __div__(self,other):     def __div__(self,other):
1435         """         """
1436         Scales by (1/C{other}) from the right.         Scales by (1/``other``) from the right.
1437    
1438         @param other: scaling factor         :param other: scaling factor
1439         @type other: C{float}         :type other: ``float``
1440         @return: itemwise self/other         :return: itemwise self/other
1441         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1442         """         """
1443         return self*(1/other)         return self*(1/other)
1444    
1445     def __rdiv__(self,other):     def __rdiv__(self,other):
1446         """        """
1447         Scales by (1/C{other}) from the left.        Scales by (1/``other``) from the left.
1448    
1449         @param other: scaling factor        :param other: scaling factor
1450         @type other: C{float}        :type other: ``float``
1451         @return: itemwise other/self        :return: itemwise other/self
1452         @rtype: L{ArithmeticTuple}        :rtype: `ArithmeticTuple`
1453         """        """
1454         out=[]        out=[]
1455         try:        try:
1456             l=len(other)            l=len(other)
1457             if l!=len(self):            if l!=len(self):
1458                 raise ValueError,"length of arguments don't match."                raise ValueError("length of arguments don't match.")
1459             for i in range(l): out.append(other[i]/self[i])            
1460         except TypeError:            for i in range(l):
1461             for i in range(len(self)): out.append(other/self[i])                  if self.__isEmpty(self[i]):
1462         return ArithmeticTuple(*tuple(out))                      raise ZeroDivisionError("in component %s"%i)
1463                    else:
1464                        if self.__isEmpty(other[i]):
1465                            out.append(escore.Data())
1466                        else:
1467                            out.append(other[i]/self[i])
1468          except TypeError:
1469              for i in range(len(self)):
1470                    if self.__isEmpty(self[i]):
1471                        raise ZeroDivisionError("in component %s"%i)
1472                    else:
1473                        if self.__isEmpty(other):
1474                            out.append(escore.Data())
1475                        else:
1476                            out.append(other/self[i])
1477          return ArithmeticTuple(*tuple(out))
1478    
1479     def __iadd__(self,other):     def __iadd__(self,other):
1480         """        """
1481         Inplace addition of C{other} to self.        Inplace addition of ``other`` to self.
1482    
1483         @param other: increment        :param other: increment
1484         @type other: C{ArithmeticTuple}        :type other: ``ArithmeticTuple``
1485         """        """
1486         if len(self) != len(other):        if len(self) != len(other):
1487             raise ValueError,"tuple lengths must match."            raise ValueError("tuple lengths must match.")
1488         for i in range(len(self)):        for i in range(len(self)):
1489             self.__items[i]+=other[i]            if self.__isEmpty(self.__items[i]):
1490         return self                self.__items[i]=other[i]
1491              else:
1492                  self.__items[i]+=other[i]
1493                  
1494          return self
1495    
1496     def __add__(self,other):     def __add__(self,other):
1497         """        """
1498         Adds C{other} to self.        Adds ``other`` to self.
1499    
1500         @param other: increment        :param other: increment
1501         @type other: C{ArithmeticTuple}        :type other: ``ArithmeticTuple``
1502         """        """
1503         out=[]        out=[]
1504         try:        try:
1505             l=len(other)            l=len(other)
1506             if l!=len(self):            if l!=len(self):
1507                 raise ValueError,"length of arguments don't match."                raise ValueError("length of arguments don't match.")
1508             for i in range(l): out.append(self[i]+other[i])            for i in range(l):
1509         except TypeError:                  if self.__isEmpty(self[i]):
1510             for i in range(len(self)): out.append(self[i]+other)                      out.append(other[i])
1511         return ArithmeticTuple(*tuple(out))                  elif self.__isEmpty(other[i]):
1512                        out.append(self[i])
1513                    else:
1514                        out.append(self[i]+other[i])
1515          except TypeError:
1516                for i in range(len(self)):    
1517                    if self.__isEmpty(self[i]):
1518                        out.append(other)
1519                    elif self.__isEmpty(other):
1520                        out.append(self[i])
1521                    else:
1522                        out.append(self[i]+other)
1523          return ArithmeticTuple(*tuple(out))
1524    
1525     def __sub__(self,other):     def __sub__(self,other):
1526         """        """
1527         Subtracts C{other} from self.        Subtracts ``other`` from self.
1528    
1529         @param other: decrement        :param other: decrement
1530         @type other: C{ArithmeticTuple}        :type other: ``ArithmeticTuple``
1531         """        """
1532         out=[]        out=[]
1533         try:        try:
1534             l=len(other)            l=len(other)
1535             if l!=len(self):            if l!=len(self):
1536                 raise ValueError,"length of arguments don't match."                raise ValueError("length of arguments don't match.")
1537             for i in range(l): out.append(self[i]-other[i])            for i in range(l):
1538         except TypeError:                  if self.__isEmpty(other[i]):
1539             for i in range(len(self)): out.append(self[i]-other)                      out.append(self[i])
1540         return ArithmeticTuple(*tuple(out))                  elif self.__isEmpty(self[i]):
1541                        out.append(-other[i])
1542                    else:
1543                        out.append(self[i]-other[i])
1544          except TypeError:
1545                for i in range(len(self)):    
1546                    if  self.__isEmpty(other):
1547                        out.append(self[i])
1548                    elif self.__isEmpty(self[i]):
1549                        out.append(-other)
1550                    else:
1551                        out.append(self[i]-other)
1552                        
1553          return ArithmeticTuple(*tuple(out))
1554    
1555     def __isub__(self,other):     def __isub__(self,other):
1556         """        """
1557         Inplace subtraction of C{other} from self.        Inplace subtraction of ``other`` from self.
1558    
1559         @param other: decrement        :param other: decrement
1560         @type other: C{ArithmeticTuple}        :type other: ``ArithmeticTuple``
1561         """        """
1562         if len(self) != len(other):        if len(self) != len(other):
1563             raise ValueError,"tuple length must match."            raise ValueError("tuple length must match.")
1564         for i in range(len(self)):        for i in range(len(self)):
1565             self.__items[i]-=other[i]            if not self.__isEmpty(other[i]):
1566         return self                if self.__isEmpty(self.__items[i]):
1567                      self.__items[i]=-other[i]
1568                  else:
1569                      self.__items[i]=other[i]
1570          return self
1571    
1572     def __neg__(self):     def __neg__(self):
1573         """        """
1574         Negates values.        Negates values.
1575         """        """
1576         out=[]        out=[]
1577         for i in range(len(self)):        for i in range(len(self)):
1578             out.append(-self[i])            if self.__isEmpty(self[i]):
1579         return ArithmeticTuple(*tuple(out))                out.append(escore.Data())
1580              else:
1581                  out.append(-self[i])
1582              
1583          return ArithmeticTuple(*tuple(out))
1584       def __isEmpty(self, d):
1585        if isinstance(d, escore.Data):
1586            return d.isEmpty()
1587        else:
1588            return False
1589                    
1590       def __str__(self):
1591        s="("
1592        for i in self:
1593          s=s+str(i)+", "
1594        s=s+")"
1595        return s
1596    
1597  class HomogeneousSaddlePointProblem(object):  class HomogeneousSaddlePointProblem(object):
1598        """        """
1599        This class provides a framework for solving linear homogeneous saddle        This class provides a framework for solving linear homogeneous saddle
1600        point problems of the form::        point problems of the form::
1601    
1602            M{Av+B^*p=f}            *Av+B^*p=f*
1603            M{Bv     =0}            *Bv     =0*
1604    
1605        for the unknowns M{v} and M{p} and given operators M{A} and M{B} and        for the unknowns *v* and *p* and given operators *A* and *B* and
1606        given right hand side M{f}. M{B^*} is the adjoint operator of M{B}.        given right hand side *f*. *B^** is the adjoint operator of *B*.
1607          *A* may depend weakly on *v* and *p*.
1608        """        """
1609        def __init__(self, adaptSubTolerance=True, **kwargs):        def __init__(self, **kwargs):
     """  
     initializes the saddle point problem  
       
     @param adaptSubTolerance: If True the tolerance for subproblem is set automatically.  
     @type adaptSubTolerance: C{bool}  
     """  
         self.setTolerance()  
         self.setAbsoluteTolerance()  
     self.__adaptSubTolerance=adaptSubTolerance  
       #=============================================================  
       def initialize(self):  
1610          """          """
1611          Initializes the problem (overwrite).          initializes the saddle point problem
1612          """          """
1613          pass          self.resetControlParameters()
1614            self.setTolerance()
1615            self.setAbsoluteTolerance()
1616          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):
1617             """
1618             sets a control parameter
1619    
1620             :param K_p: initial value for constant to adjust pressure tolerance
1621             :type K_p: ``float``
1622             :param K_v: initial value for constant to adjust velocity tolerance
1623             :type K_v: ``float``
1624             :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1625             :type rtol_max: ``float``
1626             :param chi_max: maximum tolerable converegence rate.
1627             :type chi_max: ``float``
1628             :param reduction_factor: reduction factor for adjustment factors.
1629             :type reduction_factor: ``float``
1630             """
1631             self.setControlParameter(K_p, K_v, rtol_max, rtol_min, chi_max, reduction_factor, theta)
1632    
1633          def setControlParameter(self,K_p=None, K_v=None, rtol_max=None, rtol_min=None, chi_max=None, reduction_factor=None, theta=None):
1634             """
1635             sets a control parameter
1636    
1637    
1638             :param K_p: initial value for constant to adjust pressure tolerance
1639             :type K_p: ``float``
1640             :param K_v: initial value for constant to adjust velocity tolerance
1641             :type K_v: ``float``
1642             :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1643             :type rtol_max: ``float``
1644             :param chi_max: maximum tolerable converegence rate.
1645             :type chi_max: ``float``
1646             :type reduction_factor: ``float``
1647             """
1648             if not K_p == None:
1649                if K_p<1:
1650                   raise ValueError("K_p need to be greater or equal to 1.")
1651             else:
1652                K_p=self.__K_p
1653    
1654             if not K_v == None:
1655                if K_v<1:
1656                   raise ValueError("K_v need to be greater or equal to 1.")
1657             else:
1658                K_v=self.__K_v
1659    
1660             if not rtol_max == None:
1661                if rtol_max<=0 or rtol_max>=1:
1662                   raise ValueError("rtol_max needs to be positive and less than 1.")
1663             else:
1664                rtol_max=self.__rtol_max
1665    
1666             if not rtol_min == None:
1667                if rtol_min<=0 or rtol_min>=1:
1668                   raise ValueError("rtol_min needs to be positive and less than 1.")
1669             else:
1670                rtol_min=self.__rtol_min
1671    
1672             if not chi_max == None:
1673                if chi_max<=0 or chi_max>=1:
1674                   raise ValueError("chi_max needs to be positive and less than 1.")
1675             else:
1676                chi_max = self.__chi_max
1677    
1678             if not reduction_factor == None:
1679                if reduction_factor<=0 or reduction_factor>1:
1680                   raise ValueError("reduction_factor need to be between zero and one.")
1681             else:
1682                reduction_factor=self.__reduction_factor
1683    
1684             if not theta == None:
1685                if theta<=0 or theta>1:
1686                   raise ValueError("theta need to be between zero and one.")
1687             else:
1688                theta=self.__theta
1689    
1690             if rtol_min>=rtol_max:
1691                 raise ValueError("rtol_max = %e needs to be greater than rtol_min = %e"%(rtol_max,rtol_min))
1692             self.__chi_max = chi_max
1693             self.__rtol_max = rtol_max
1694             self.__K_p = K_p
1695             self.__K_v = K_v
1696             self.__reduction_factor = reduction_factor
1697             self.__theta = theta
1698             self.__rtol_min=rtol_min
1699    
1700          #=============================================================
1701        def inner_pBv(self,p,Bv):        def inner_pBv(self,p,Bv):
1702           """           """
1703           Returns inner product of element p and Bv (overwrite).           Returns inner product of element p and Bv (overwrite).
1704    
1705           @param p: a pressure increment           :param p: a pressure increment
1706           @param v: a residual           :param Bv: a residual
1707           @return: inner product of element p and Bv           :return: inner product of element p and Bv
1708           @rtype: C{float}           :rtype: ``float``
1709           @note: used if PCG is applied.           :note: used if PCG is applied.
1710           """           """
1711           raise NotImplementedError,"no inner product for p and Bv implemented."           raise NotImplementedError("no inner product for p and Bv implemented.")
1712    
1713        def inner_p(self,p0,p1):        def inner_p(self,p0,p1):
1714           """           """
1715           Returns inner product of p0 and p1 (overwrite).           Returns inner product of p0 and p1 (overwrite).
1716    
1717           @param p0: a pressure           :param p0: a pressure
1718           @param p1: a pressure           :param p1: a pressure
1719           @return: inner product of p0 and p1           :return: inner product of p0 and p1
1720           @rtype: C{float}           :rtype: ``float``
1721           """           """
1722           raise NotImplementedError,"no inner product for p implemented."           raise NotImplementedError("no inner product for p implemented.")
1723        
1724        def norm_v(self,v):        def norm_v(self,v):
1725           """           """
1726           Returns the norm of v (overwrite).           Returns the norm of v (overwrite).
1727    
1728           @param v: a velovity           :param v: a velovity
1729           @return: norm of v           :return: norm of v
1730           @rtype: non-negative C{float}           :rtype: non-negative ``float``
1731           """           """
1732           raise NotImplementedError,"no norm of v implemented."           raise NotImplementedError("no norm of v implemented.")
1733        def getV(self, p, v0):        def getDV(self, p, v, tol):
1734           """           """
1735           return the value for v for a given p (overwrite)           return a correction to the value for a given v and a given p with accuracy `tol` (overwrite)
1736    
1737           @param p: a pressure           :param p: pressure
1738           @param v0: a initial guess for the value v to return.           :param v: pressure
1739           @return: v given as M{v= A^{-1} (f-B^*p)}           :return: dv given as *dv= A^{-1} (f-A v-B^*p)*
1740             :note: Only *A* may depend on *v* and *p*
1741           """           """
1742           raise NotImplementedError,"no v calculation implemented."           raise NotImplementedError("no dv calculation implemented.")
1743    
1744                    
1745        def Bv(self,v):        def Bv(self,v, tol):
1746          """          """
1747          Returns Bv (overwrite).          Returns Bv with accuracy `tol` (overwrite)
1748    
1749          @rtype: equal to the type of p          :rtype: equal to the type of p
1750          @note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1751          """          """
1752          raise NotImplementedError, "no operator B implemented."          raise NotImplementedError("no operator B implemented.")
1753    
1754        def norm_Bv(self,Bv):        def norm_Bv(self,Bv):
1755          """          """
1756          Returns the norm of Bv (overwrite).          Returns the norm of Bv (overwrite).
1757    
1758          @rtype: equal to the type of p          :rtype: equal to the type of p
1759          @note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1760          """          """
1761          raise NotImplementedError, "no norm of Bv implemented."          raise NotImplementedError("no norm of Bv implemented.")
1762    
1763        def solve_AinvBt(self,p):        def solve_AinvBt(self,dp, tol):
1764           """           """
1765           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}           Solves *A dv=B^*dp* with accuracy `tol`
          (overwrite).  
1766    
1767           @param p: a pressure increment           :param dp: a pressure increment
1768           @return: the solution of M{Av=B^*p}           :return: the solution of *A dv=B^*dp*
1769           @note: boundary conditions on v should be zero!           :note: boundary conditions on dv should be zero! *A* is the operator used in ``getDV`` and must not be altered.
1770           """           """
1771           raise NotImplementedError,"no operator A implemented."           raise NotImplementedError("no operator A implemented.")
1772    
1773        def solve_prec(self,Bv):        def solve_prec(self,Bv, tol):
1774           """           """
1775           Provides a preconditioner for M{BA^{-1}B^*} applied to Bv with accuracy           Provides a preconditioner for *(BA^{-1}B^ * )* applied to Bv with accuracy `tol`
          L{self.getSubProblemTolerance()} (overwrite).  
1776    
1777           @rtype: equal to the type of p           :rtype: equal to the type of p
1778           @note: boundary conditions on p should be zero!           :note: boundary conditions on p should be zero!
          """  
          raise NotImplementedError,"no preconditioner for Schur complement implemented."  
       def setSubProblemTolerance(self):  
          """  
      Updates the tolerance for subproblems  
      @note: method is typically the method is overwritten.  
1779           """           """
1780           pass           raise NotImplementedError("no preconditioner for Schur complement implemented.")
1781        #=============================================================        #=============================================================
1782        def __Aprod_PCG(self,p):        def __Aprod_PCG(self,dp):
1783            dv=self.solve_AinvBt(p)            dv=self.solve_AinvBt(dp, self.__subtol)
1784            return ArithmeticTuple(dv,self.Bv(dv))            return ArithmeticTuple(dv,self.Bv(dv, self.__subtol))
1785    
1786        def __inner_PCG(self,p,r):        def __inner_PCG(self,p,r):
1787           return self.inner_pBv(p,r[1])           return self.inner_pBv(p,r[1])
1788    
1789        def __Msolve_PCG(self,r):        def __Msolve_PCG(self,r):
1790            return self.solve_prec(r[1])            return self.solve_prec(r[1], self.__subtol)
1791        #=============================================================        #=============================================================
1792        def __Aprod_GMRES(self,p):        def __Aprod_GMRES(self,p):
1793            return self.solve_prec(self.Bv(self.solve_AinvBt(p)))            return self.solve_prec(self.Bv(self.solve_AinvBt(p, self.__subtol), self.__subtol), self.__subtol)
1794        def __inner_GMRES(self,p0,p1):        def __inner_GMRES(self,p0,p1):
1795           return self.inner_p(p0,p1)           return self.inner_p(p0,p1)
1796    
1797        #=============================================================        #=============================================================
1798        def norm_p(self,p):        def norm_p(self,p):
1799            """            """
1800            calculates the norm of C{p}            calculates the norm of ``p``
1801                        
1802            @param p: a pressure            :param p: a pressure
1803            @return: the norm of C{p} using the inner product for pressure            :return: the norm of ``p`` using the inner product for pressure
1804            @rtype: C{float}            :rtype: ``float``
1805            """            """
1806            f=self.inner_p(p,p)            f=self.inner_p(p,p)
1807            if f<0: raise ValueError,"negative pressure norm."            if f<0: raise ValueError("negative pressure norm.")
1808            return math.sqrt(f)            return math.sqrt(f)
1809        def adaptSubTolerance(self):            
       """  
       Returns True if tolerance adaption for subproblem is choosen.  
       """  
           self.__adaptSubTolerance  
         
1810        def solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):        def solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1811           """           """
1812           Solves the saddle point problem using initial guesses v and p.           Solves the saddle point problem using initial guesses v and p.
1813    
1814           @param v: initial guess for velocity           :param v: initial guess for velocity
1815           @param p: initial guess for pressure           :param p: initial guess for pressure
1816           @type v: L{Data}           :type v: `Data`
1817           @type p: L{Data}           :type p: `Data`
1818           @param usePCG: indicates the usage of the PCG rather than GMRES scheme.           :param usePCG: indicates the usage of the PCG rather than GMRES scheme.
1819           @param max_iter: maximum number of iteration steps per correction           :param max_iter: maximum number of iteration steps per correction
1820                            attempt                            attempt
1821           @param verbose: if True, shows information on the progress of the           :param verbose: if True, shows information on the progress of the
1822                           saddlepoint problem solver.                           saddlepoint problem solver.
1823           @param iter_restart: restart the iteration after C{iter_restart} steps           :param iter_restart: restart the iteration after ``iter_restart`` steps
1824                                (only used if useUzaw=False)                                (only used if useUzaw=False)
1825           @type usePCG: C{bool}           :type usePCG: ``bool``
1826           @type max_iter: C{int}           :type max_iter: ``int``
1827           @type verbose: C{bool}           :type verbose: ``bool``
1828           @type iter_restart: C{int}           :type iter_restart: ``int``
1829           @rtype: C{tuple} of L{Data} objects           :rtype: ``tuple`` of `Data` objects
1830             :note: typically this method is overwritten by a subclass. It provides a wrapper for the ``_solve`` method.
1831             """
1832             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)
1833    
1834          def _solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1835             """
1836             see `_solve` method.
1837           """           """
1838           self.verbose=verbose           self.verbose=verbose
1839           rtol=self.getTolerance()           rtol=self.getTolerance()
1840           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
1841       if self.adaptSubTolerance(): self.setSubProblemTolerance()  
1842             K_p=self.__K_p
1843             K_v=self.__K_v
1844           correction_step=0           correction_step=0
1845           converged=False           converged=False
1846             chi=None
1847             eps=None
1848    
1849             if self.verbose: print(("HomogeneousSaddlePointProblem: start iteration: rtol= %e, atol=%e"%(rtol, atol)))
1850           while not converged:           while not converged:
               # calculate velocity for current pressure:  
               v=self.getV(p,v)  
               Bv=self.Bv(v)  
               norm_v=self.norm_v(v)  
               norm_Bv=self.norm_Bv(Bv)  
               ATOL=norm_v*rtol+atol  
               if self.verbose: print "HomogeneousSaddlePointProblem: norm v= %e, norm_Bv= %e, tolerance = %e."%(norm_v, norm_Bv,ATOL)  
               if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."  
               if norm_Bv <= ATOL:  
                  converged=True  
               else:  
                  correction_step+=1  
                  if correction_step>max_correction_steps:  
                       raise CorrectionFailed,"Given up after %d correction steps."%correction_step  
                  dp=self.solve_prec(Bv)  
                  if usePCG:  
                    norm2=self.inner_pBv(dp,Bv)  
                    if norm2<0: raise ValueError,"negative PCG norm."  
                    norm2=math.sqrt(norm2)  
                  else:  
                    norm2=self.norm_p(dp)  
                  ATOL_ITER=ATOL/norm_Bv*norm2*0.5  
                  if self.verbose: print "HomogeneousSaddlePointProblem: tolerance for solver: %e"%ATOL_ITER  
                  if usePCG:  
                        p,v0,a_norm=PCG(ArithmeticTuple(v,Bv),self.__Aprod_PCG,p,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL_ITER, rtol=0.,iter_max=max_iter, verbose=self.verbose)  
                  else:  
                        p=GMRES(dp,self.__Aprod_GMRES, p, self.__inner_GMRES,atol=ATOL_ITER, rtol=0.,iter_max=max_iter, iter_restart=iter_restart, verbose=self.verbose)  
          if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached."  
      return v,p  
1851    
1852                 # get tolerance for velecity increment:
1853                 if chi == None:
1854                    rtol_v=self.__rtol_max
1855                 else:
1856                    rtol_v=min(chi/K_v,self.__rtol_max)
1857                 rtol_v=max(rtol_v, self.__rtol_min)
1858                 if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: rtol_v= %e"%(correction_step,rtol_v)))
1859                 # get velocity increment:
1860                 dv1=self.getDV(p,v,rtol_v)
1861                 v1=v+dv1
1862                 Bv1=self.Bv(v1, rtol_v)
1863                 norm_Bv1=self.norm_Bv(Bv1)
1864                 norm_dv1=self.norm_v(dv1)
1865                 if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: norm_Bv1 = %e, norm_dv1 = %e"%(correction_step, norm_Bv1, norm_dv1)))
1866                 if norm_dv1*self.__theta < norm_Bv1:
1867                    # get tolerance for pressure increment:
1868                    large_Bv1=True
1869                    if chi == None or eps == None:
1870                       rtol_p=self.__rtol_max
1871                    else:
1872                       rtol_p=min(chi**2*eps/K_p/norm_Bv1, self.__rtol_max)
1873                    self.__subtol=max(rtol_p**2, self.__rtol_min)
1874                    if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: rtol_p= %e"%(correction_step,rtol_p)))
1875                    # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1
1876                    if usePCG:
1877                        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)
1878                        v2=r[0]
1879                        Bv2=r[1]
1880                    else:
1881                        # don't use!!!!
1882                        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)
1883                        dv2=self.solve_AinvBt(dp, self.__subtol)
1884                        v2=v1-dv2
1885                        Bv2=self.Bv(v2, self.__subtol)
1886                    p2=p+dp
1887                 else:
1888                    large_Bv1=False
1889                    v2=v1
1890                    p2=p
1891                 # update business:
1892                 norm_dv2=self.norm_v(v2-v)
1893                 norm_v2=self.norm_v(v2)
1894                 if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: v2 = %e, norm_dv2 = %e"%(correction_step, norm_v2, self.norm_v(v2-v))))
1895                 eps, eps_old = max(norm_Bv1, norm_dv2), eps
1896                 if eps_old == None:
1897                      chi, chi_old = None, chi
1898                 else:
1899                      chi, chi_old = min(eps/ eps_old, self.__chi_max), chi
1900                 if eps != None:
1901                     if chi !=None:
1902                        if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: convergence rate = %e, correction = %e"%(correction_step,chi, eps)))
1903                     else:
1904                        if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: correction = %e"%(correction_step, eps)))
1905                 if eps <= rtol*norm_v2+atol :
1906                     converged = True
1907                 else:
1908                     if correction_step>=max_correction_steps:
1909                          raise CorrectionFailed("Given up after %d correction steps."%correction_step)
1910                     if chi_old!=None:
1911                        K_p=max(1,self.__reduction_factor*K_p,(chi-chi_old)/chi_old**2*K_p)
1912                        K_v=max(1,self.__reduction_factor*K_v,(chi-chi_old)/chi_old**2*K_p)
1913                        if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: new adjustment factor K = %e"%(correction_step,K_p)))
1914                 correction_step+=1
1915                 v,p =v2, p2
1916             if self.verbose: print(("HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step))
1917             return v,p
1918        #========================================================================        #========================================================================
1919        def setTolerance(self,tolerance=1.e-4):        def setTolerance(self,tolerance=1.e-4):
1920           """           """
1921           Sets the relative tolerance for (v,p).           Sets the relative tolerance for (v,p).
1922    
1923           @param tolerance: tolerance to be used           :param tolerance: tolerance to be used
1924           @type tolerance: non-negative C{float}           :type tolerance: non-negative ``float``
1925           """           """
1926           if tolerance<0:           if tolerance<0:
1927               raise ValueError,"tolerance must be positive."               raise ValueError("tolerance must be positive.")
1928           self.__rtol=tolerance           self.__rtol=tolerance
1929    
1930        def getTolerance(self):        def getTolerance(self):
1931           """           """
1932           Returns the relative tolerance.           Returns the relative tolerance.
1933    
1934           @return: relative tolerance           :return: relative tolerance
1935           @rtype: C{float}           :rtype: ``float``
1936           """           """
1937           return self.__rtol           return self.__rtol
1938    
# Line 1698  class HomogeneousSaddlePointProblem(obje Line 1940  class HomogeneousSaddlePointProblem(obje
1940           """           """
1941           Sets the absolute tolerance.           Sets the absolute tolerance.
1942    
1943           @param tolerance: tolerance to be used           :param tolerance: tolerance to be used
1944           @type tolerance: non-negative C{float}           :type tolerance: non-negative ``float``
1945           """           """
1946           if tolerance<0:           if tolerance<0:
1947               raise ValueError,"tolerance must be non-negative."               raise ValueError("tolerance must be non-negative.")
1948           self.__atol=tolerance           self.__atol=tolerance
1949    
1950        def getAbsoluteTolerance(self):        def getAbsoluteTolerance(self):
1951           """           """
1952           Returns the absolute tolerance.           Returns the absolute tolerance.
1953    
1954           @return: absolute tolerance           :return: absolute tolerance
1955           @rtype: C{float}           :rtype: ``float``
1956           """           """
1957           return self.__atol           return self.__atol
1958    
       def getSubProblemTolerance(self):  
          """  
          Sets the relative tolerance to solve the subproblem(s).  
   
          @param rtol: relative tolerence  
          @type rtol: positive C{float}  
          """  
          return max(200.*util.EPSILON,self.getTolerance()**2)  
1959    
1960  def MaskFromBoundaryTag(domain,*tags):  def MaskFromBoundaryTag(domain,*tags):
1961     """     """
# Line 1730  def MaskFromBoundaryTag(domain,*tags): Line 1964  def MaskFromBoundaryTag(domain,*tags):
1964    
1965     Usage: m=MaskFromBoundaryTag(domain, "left", "right")     Usage: m=MaskFromBoundaryTag(domain, "left", "right")
1966    
1967     @param domain: domain to be used     :param domain: domain to be used
1968     @type domain: L{escript.Domain}     :type domain: `escript.Domain`
1969     @param tags: boundary tags     :param tags: boundary tags
1970     @type tags: C{str}     :type tags: ``str``
1971     @return: a mask which marks samples that are touching the boundary tagged     :return: a mask which marks samples that are touching the boundary tagged
1972              by any of the given tags              by any of the given tags
1973     @rtype: L{escript.Data} of rank 0     :rtype: `escript.Data` of rank 0
1974     """     """
1975     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1976     d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))     d=escore.Scalar(0.,escore.FunctionOnBoundary(domain))
1977     for t in tags: d.setTaggedValue(t,1.)     for t in tags: d.setTaggedValue(t,1.)
1978     pde.setValue(y=d)     pde.setValue(y=d)
1979     return util.whereNonZero(pde.getRightHandSide())     return util.whereNonZero(pde.getRightHandSide())
# Line 1751  def MaskFromTag(domain,*tags): Line 1985  def MaskFromTag(domain,*tags):
1985    
1986     Usage: m=MaskFromTag(domain, "ham")     Usage: m=MaskFromTag(domain, "ham")
1987    
1988     @param domain: domain to be used     :param domain: domain to be used
1989     @type domain: L{escript.Domain}     :type domain: `escript.Domain`
1990     @param tags: boundary tags     :param tags: boundary tags
1991     @type tags: C{str}     :type tags: ``str``
1992     @return: a mask which marks samples that are touching the boundary tagged     :return: a mask which marks samples that are touching the boundary tagged
1993              by any of the given tags              by any of the given tags
1994     @rtype: L{escript.Data} of rank 0     :rtype: `escript.Data` of rank 0
1995     """     """
1996     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1997     d=escript.Scalar(0.,escript.Function(domain))     d=escore.Scalar(0.,escore.Function(domain))
1998     for t in tags: d.setTaggedValue(t,1.)     for t in tags: d.setTaggedValue(t,1.)
1999     pde.setValue(Y=d)     pde.setValue(Y=d)
2000     return util.whereNonZero(pde.getRightHandSide())     return util.whereNonZero(pde.getRightHandSide())

Legend:
Removed from v.2548  
changed lines
  Added in v.4446

  ViewVC Help
Powered by ViewVC 1.1.26