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

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

  ViewVC Help
Powered by ViewVC 1.1.26