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

Legend:
Removed from v.351  
changed lines
  Added in v.2719

  ViewVC Help
Powered by ViewVC 1.1.26