/[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 1481 by artak, Wed Apr 9 00:45:47 2008 UTC revision 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC
# Line 1  Line 1 
1    
2    ########################################################
3  #  #
4  # $Id$  # Copyright (c) 2003-2009 by University of Queensland
5  #  # Earth Systems Science Computational Center (ESSCC)
6  #######################################################  # http://www.uq.edu.au/esscc
 #  
 #           Copyright 2003-2007 by ACceSS MNRF  
 #       Copyright 2007 by University of Queensland  
 #  
 #                http://esscc.uq.edu.au  
 #        Primary Business: Queensland, Australia  
 #  Licensed under the Open Software License version 3.0  
 #     http://www.opensource.org/licenses/osl-3.0.php  
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-2008 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      - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme
30    
31  @var __author__: name of author  @var __author__: name of author
# Line 32  Currently includes: Line 37  Currently includes:
37  """  """
38    
39  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
 __copyright__="""  Copyright (c) 2006 by ACcESS MNRF  
                     http://www.access.edu.au  
                 Primary Business: Queensland, Australia"""  
 __license__="""Licensed under the Open Software License version 3.0  
              http://www.opensource.org/licenses/osl-3.0.php"""  
 __url__="http://www.iservo.edu.au/esys"  
 __version__="$Revision$"  
 __date__="$Date$"  
40    
41    
42  import escript  import escript
43  import linearPDEs  import linearPDEs
44  import numarray  import numpy
45  import util  import util
46  import math  import math
47    
# Line 55  import math Line 52  import math
52    
53  class TimeIntegrationManager:  class TimeIntegrationManager:
54    """    """
55    a simple mechanism to manage time dependend values.    A simple mechanism to manage time dependend values.
56    
57    typical usage is::    Typical usage is::
58    
59       dt=0.1 # time increment       dt=0.1 # time increment
60       tm=TimeIntegrationManager(inital_value,p=1)       tm=TimeIntegrationManager(inital_value,p=1)
# Line 71  class TimeIntegrationManager: Line 68  class TimeIntegrationManager:
68    """    """
69    def __init__(self,*inital_values,**kwargs):    def __init__(self,*inital_values,**kwargs):
70       """       """
71       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 C{inital_values} are the initial values
72         and p is the order used for extrapolation.
73       """       """
74       if kwargs.has_key("p"):       if kwargs.has_key("p"):
75              self.__p=kwargs["p"]              self.__p=kwargs["p"]
# Line 88  class TimeIntegrationManager: Line 86  class TimeIntegrationManager:
86    
87    def getTime(self):    def getTime(self):
88        return self.__t        return self.__t
89    
90    def getValue(self):    def getValue(self):
91        out=self.__v_mem[0]        out=self.__v_mem[0]
92        if len(out)==1:        if len(out)==1:
# Line 97  class TimeIntegrationManager: Line 96  class TimeIntegrationManager:
96    
97    def checkin(self,dt,*values):    def checkin(self,dt,*values):
98        """        """
99        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.
100        """        """
101        o=min(self.__order+1,self.__p)        o=min(self.__order+1,self.__p)
102        self.__order=min(self.__order+1,self.__p)        self.__order=min(self.__order+1,self.__p)
# Line 114  class TimeIntegrationManager: Line 113  class TimeIntegrationManager:
113    
114    def extrapolate(self,dt):    def extrapolate(self,dt):
115        """        """
116        extrapolates to dt forward in time.        Extrapolates to C{dt} forward in time.
117        """        """
118        if self.__order==0:        if self.__order==0:
119           out=self.__v_mem[0]           out=self.__v_mem[0]
# Line 129  class TimeIntegrationManager: Line 128  class TimeIntegrationManager:
128           return out[0]           return out[0]
129        else:        else:
130           return out           return out
131    
132    
133  class Projector:  class Projector:
134    """    """
135    The Projector is a factory which projects a discontiuous function onto a    The Projector is a factory which projects a discontinuous function onto a
136    continuous function on the a given domain.    continuous function on a given domain.
137    """    """
138    def __init__(self, domain, reduce = True, fast=True):    def __init__(self, domain, reduce=True, fast=True):
139      """      """
140      Create a continuous function space projector for a domain.      Creates a continuous function space projector for a domain.
141    
142      @param domain: Domain of the projection.      @param domain: Domain of the projection.
143      @param reduce: Flag to reduce projection order (default is True)      @param reduce: Flag to reduce projection order
144      @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
145      """      """
146      self.__pde = linearPDEs.LinearPDE(domain)      self.__pde = linearPDEs.LinearPDE(domain)
147      if fast:      if fast:
148        self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)          self.__pde.getSolverOptions().setSolverMethod(linearPDEs.SolverOptions.LUMPING)
149      self.__pde.setSymmetryOn()      self.__pde.setSymmetryOn()
150      self.__pde.setReducedOrderTo(reduce)      self.__pde.setReducedOrderTo(reduce)
151      self.__pde.setValue(D = 1.)      self.__pde.setValue(D = 1.)
152      return      return
153      def getSolverOptions(self):
154        """
155        Returns the solver options of the PDE solver.
156        
157        @rtype: L{linearPDEs.SolverOptions}
158        """
159    
160    def __call__(self, input_data):    def __call__(self, input_data):
161      """      """
162      Projects input_data onto a continuous function      Projects C{input_data} onto a continuous function.
163    
164      @param input_data: The input_data to be projected.      @param input_data: the data to be projected
165      """      """
166      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
167      self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())      self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())
# Line 189  class Projector: Line 194  class Projector:
194    
195  class NoPDE:  class NoPDE:
196       """       """
197       solves the following problem for u:       Solves the following problem for u:
198    
199       M{kronecker[i,j]*D[j]*u[j]=Y[i]}       M{kronecker[i,j]*D[j]*u[j]=Y[i]}
200    
201       with constraint       with constraint
202    
203       M{u[j]=r[j]}  where M{q[j]>0}       M{u[j]=r[j]}  where M{q[j]>0}
204    
205       where D, Y, r and q are given functions of rank 1.       where M{D}, M{Y}, M{r} and M{q} are given functions of rank 1.
206    
207       In the case of scalars this takes the form       In the case of scalars this takes the form
208    
209       M{D*u=Y}       M{D*u=Y}
210    
211       with constraint       with constraint
212    
213       M{u=r}  where M{q>0}       M{u=r} where M{q>0}
214    
215       where D, Y, r and q are given scalar functions.       where M{D}, M{Y}, M{r} and M{q} are given scalar functions.
216    
217       The constraint is overwriting any other condition.       The constraint overwrites any other condition.
218    
219       @note: This class is similar to the L{linearPDEs.LinearPDE} class with A=B=C=X=0 but has the intention       @note: This class is similar to the L{linearPDEs.LinearPDE} class with
220              that all input parameter are given in L{Solution} or L{ReducedSolution}. The whole              A=B=C=X=0 but has the intention that all input parameters are given
221              thing is a bit strange and I blame Robert.Woodcock@csiro.au for this.              in L{Solution} or L{ReducedSolution}.
222       """       """
223         # The whole thing is a bit strange and I blame Rob Woodcock (CSIRO) for
224         # this.
225       def __init__(self,domain,D=None,Y=None,q=None,r=None):       def __init__(self,domain,D=None,Y=None,q=None,r=None):
226           """           """
227           initialize the problem           Initializes the problem.
228    
229           @param domain: domain of the PDE.           @param domain: domain of the PDE
230           @type domain: L{Domain}           @type domain: L{Domain}
231           @param D: coefficient of the solution.           @param D: coefficient of the solution
232           @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}           @type D: C{float}, C{int}, C{numpy.ndarray}, L{Data}
233           @param Y: right hand side           @param Y: right hand side
234           @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}           @type Y: C{float}, C{int}, C{numpy.ndarray}, L{Data}
235           @param q: location of constraints           @param q: location of constraints
236           @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}           @type q: C{float}, C{int}, C{numpy.ndarray}, L{Data}
237           @param r: value of solution at locations of constraints           @param r: value of solution at locations of constraints
238           @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}           @type r: C{float}, C{int}, C{numpy.ndarray}, L{Data}
239           """           """
240           self.__domain=domain           self.__domain=domain
241           self.__D=D           self.__D=D
# Line 237  class NoPDE: Line 244  class NoPDE:
244           self.__r=r           self.__r=r
245           self.__u=None           self.__u=None
246           self.__function_space=escript.Solution(self.__domain)           self.__function_space=escript.Solution(self.__domain)
247    
248       def setReducedOn(self):       def setReducedOn(self):
249           """           """
250           sets the L{FunctionSpace} of the solution to L{ReducedSolution}           Sets the L{FunctionSpace} of the solution to L{ReducedSolution}.
251           """           """
252           self.__function_space=escript.ReducedSolution(self.__domain)           self.__function_space=escript.ReducedSolution(self.__domain)
253           self.__u=None           self.__u=None
254    
255       def setReducedOff(self):       def setReducedOff(self):
256           """           """
257           sets the L{FunctionSpace} of the solution to L{Solution}           Sets the L{FunctionSpace} of the solution to L{Solution}.
258           """           """
259           self.__function_space=escript.Solution(self.__domain)           self.__function_space=escript.Solution(self.__domain)
260           self.__u=None           self.__u=None
261            
262       def setValue(self,D=None,Y=None,q=None,r=None):       def setValue(self,D=None,Y=None,q=None,r=None):
263           """           """
264           assigns values to the parameters.           Assigns values to the parameters.
265    
266           @param D: coefficient of the solution.           @param D: coefficient of the solution
267           @type D: C{float}, C{int}, L{numarray.NumArray}, L{Data}           @type D: C{float}, C{int}, C{numpy.ndarray}, L{Data}
268           @param Y: right hand side           @param Y: right hand side
269           @type Y: C{float}, C{int}, L{numarray.NumArray}, L{Data}           @type Y: C{float}, C{int}, C{numpy.ndarray}, L{Data}
270           @param q: location of constraints           @param q: location of constraints
271           @type q: C{float}, C{int}, L{numarray.NumArray}, L{Data}           @type q: C{float}, C{int}, C{numpy.ndarray}, L{Data}
272           @param r: value of solution at locations of constraints           @param r: value of solution at locations of constraints
273           @type r: C{float}, C{int}, L{numarray.NumArray}, L{Data}           @type r: C{float}, C{int}, C{numpy.ndarray}, L{Data}
274           """           """
275           if not D==None:           if not D==None:
276              self.__D=D              self.__D=D
# Line 279  class NoPDE: Line 287  class NoPDE:
287    
288       def getSolution(self):       def getSolution(self):
289           """           """
290           returns the solution           Returns the solution.
291            
292           @return: the solution of the problem           @return: the solution of the problem
293           @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or L{ReducedSolution}.           @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or
294                     L{ReducedSolution}
295           """           """
296           if self.__u==None:           if self.__u==None:
297              if self.__D==None:              if self.__D==None:
# Line 299  class NoPDE: Line 308  class NoPDE:
308                  self.__u*=(1.-q)                  self.__u*=(1.-q)
309                  if not self.__r==None: self.__u+=q*self.__r                  if not self.__r==None: self.__u+=q*self.__r
310           return self.__u           return self.__u
311                
312  class Locator:  class Locator:
313       """       """
314       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
315       spatial coordinate x.         coordinate x.
316        
317       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
318       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.  
319       """       """
320    
321       def __init__(self,where,x=numarray.zeros((3,))):       def __init__(self,where,x=numpy.zeros((3,))):
322         """         """
323         Initializes a Locator to access values in Data objects on the Doamin         Initializes a Locator to access values in Data objects on the Doamin
324         or FunctionSpace where for the sample point which         or FunctionSpace for the sample point which is closest to the given
325         closest to the given point x.         point x.
326    
327         @param where: function space         @param where: function space
328         @type where: L{escript.FunctionSpace}         @type where: L{escript.FunctionSpace}
329         @param x: coefficient of the solution.         @param x: location(s) of the Locator
330         @type x: L{numarray.NumArray} or C{list} of L{numarray.NumArray}         @type x: C{numpy.ndarray} or C{list} of C{numpy.ndarray}
331         """         """
332         if isinstance(where,escript.FunctionSpace):         if isinstance(where,escript.FunctionSpace):
333            self.__function_space=where            self.__function_space=where
334         else:         else:
335            self.__function_space=escript.ContinuousFunction(where)            self.__function_space=escript.ContinuousFunction(where)
336           iterative=False
337         if isinstance(x, list):         if isinstance(x, list):
338               if len(x)==0:
339                  raise "ValueError", "At least one point must be given."
340               try:
341                 iter(x[0])
342                 iterative=True
343               except TypeError:
344                 iterative=False
345           if iterative:
346             self.__id=[]             self.__id=[]
347             for p in x:             for p in x:
348                self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())                self.__id.append(util.length(self.__function_space.getX()-p[:self.__function_space.getDim()]).minGlobalDataPoint())
# Line 353  class Locator: Line 370  class Locator:
370    
371       def getX(self):       def getX(self):
372          """          """
373      Returns the exact coordinates of the Locator.          Returns the exact coordinates of the Locator.
374      """          """
375          return self(self.getFunctionSpace().getX())          return self(self.getFunctionSpace().getX())
376    
377       def getFunctionSpace(self):       def getFunctionSpace(self):
378          """          """
379      Returns the function space of the Locator.          Returns the function space of the Locator.
380      """          """
381          return self.__function_space          return self.__function_space
382    
383       def getId(self,item=None):       def getId(self,item=None):
384          """          """
385      Returns the identifier of the location.          Returns the identifier of the location.
386      """          """
387          if item == None:          if item == None:
388             return self.__id             return self.__id
389          else:          else:
# Line 378  class Locator: Line 395  class Locator:
395    
396       def __call__(self,data):       def __call__(self,data):
397          """          """
398      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.
399      the object is returned.          """
     """  
400          return self.getValue(data)          return self.getValue(data)
401    
402       def getValue(self,data):       def getValue(self,data):
403          """          """
404      Returns the value of data at the Locator if data is a Data object          Returns the value of C{data} at the Locator if C{data} is a L{Data}
405      otherwise the object is returned.          object otherwise the object is returned.
406      """          """
407          if isinstance(data,escript.Data):          if isinstance(data,escript.Data):
408             if data.getFunctionSpace()==self.getFunctionSpace():             dat=util.interpolate(data,self.getFunctionSpace())
              dat=data  
            else:  
              dat=data.interpolate(self.getFunctionSpace())  
409             id=self.getId()             id=self.getId()
410             r=data.getRank()             r=data.getRank()
411             if isinstance(id,list):             if isinstance(id,list):
412                 out=[]                 out=[]
413                 for i in id:                 for i in id:
414                    o=data.getValueOfGlobalDataPoint(*i)                    o=numpy.array(dat.getTupleForGlobalDataPoint(*i))
415                    if data.getRank()==0:                    if data.getRank()==0:
416                       out.append(o[0])                       out.append(o[0])
417                    else:                    else:
418                       out.append(o)                       out.append(o)
419                 return out                 return out
420             else:             else:
421               out=data.getValueOfGlobalDataPoint(*id)               out=numpy.array(dat.getTupleForGlobalDataPoint(*id))
422               if data.getRank()==0:               if data.getRank()==0:
423                  return out[0]                  return out[0]
424               else:               else:
# Line 415  class Locator: Line 428  class Locator:
428    
429  class SolverSchemeException(Exception):  class SolverSchemeException(Exception):
430     """     """
431     exceptions thrown by solvers     This is a generic exception thrown by solvers.
432     """     """
433     pass     pass
434    
435  class IndefinitePreconditioner(SolverSchemeException):  class IndefinitePreconditioner(SolverSchemeException):
436     """     """
437     the preconditioner is not positive definite.     Exception thrown if the preconditioner is not positive definite.
438     """     """
439     pass     pass
440    
441  class MaxIterReached(SolverSchemeException):  class MaxIterReached(SolverSchemeException):
442     """     """
443     maxium number of iteration steps is reached.     Exception thrown if the maximum number of iteration steps is reached.
444     """     """
445     pass     pass
446  class IterationBreakDown(SolverSchemeException):  
447    class CorrectionFailed(SolverSchemeException):
448     """     """
449     iteration scheme econouters an incurable breakdown.     Exception thrown if no convergence has been achieved in the solution
450       correction scheme.
451     """     """
452     pass     pass
453  class NegativeNorm(SolverSchemeException):  
454    class IterationBreakDown(SolverSchemeException):
455     """     """
456     a norm calculation returns a negative norm.     Exception thrown if the iteration scheme encountered an incurable breakdown.
457     """     """
458     pass     pass
459    
460  class IterationHistory(object):  class NegativeNorm(SolverSchemeException):
461     """     """
462     The IterationHistory class is used to define a stopping criterium. It keeps track of the     Exception thrown if a norm calculation returns a negative norm.
    residual norms. The stoppingcriterium indicates termination if the residual norm has been reduced by  
    a given tolerance.  
463     """     """
464     def __init__(self,tolerance=math.sqrt(util.EPSILON),verbose=False):     pass
       """  
       Initialization  
   
       @param tolerance: tolerance  
       @type tolerance: positive C{float}  
       @param verbose: switches on the printing out some information  
       @type verbose: C{bool}  
       """  
       if not tolerance>0.:  
           raise ValueError,"tolerance needs to be positive."  
       self.tolerance=tolerance  
       self.verbose=verbose  
       self.history=[]  
    def stoppingcriterium(self,norm_r,r,x):  
        """  
        returns True if the C{norm_r} is C{tolerance}*C{norm_r[0]} where C{norm_r[0]}  is the residual norm at the first call.  
   
         
        @param norm_r: current residual norm  
        @type norm_r: non-negative C{float}  
        @param r: current residual (not used)  
        @param x: current solution approximation (not used)  
        @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.  
        @rtype: C{bool}  
   
        """  
        self.history.append(norm_r)  
        if self.verbose: print "iter: %s:  inner(rhat,r) = %e"%(len(self.history)-1, self.history[-1])  
        return self.history[-1]<=self.tolerance * self.history[0]  
465    
466     def stoppingcriterium2(self,norm_r,norm_b):  def PCG(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100, initial_guess=True, verbose=False):
467         """     """
468         returns True if the C{norm_r} is C{tolerance}*C{norm_b}     Solver for
469    
470             M{Ax=b}
        @param norm_r: current residual norm  
        @type norm_r: non-negative C{float}  
        @param norm_b: norm of right hand side  
        @type norm_b: non-negative C{float}  
        @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.  
        @rtype: C{bool}  
471    
472         """     with a symmetric and positive definite operator A (more details required!).
473         self.history.append(norm_r)     It uses the conjugate gradient method with preconditioner M providing an
474         if self.verbose: print "iter: %s:  norm(r) = %e"%(len(self.history)-1, self.history[-1])     approximation of A.
        return self.history[-1]<=self.tolerance * norm_b  
475    
476  def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):     The iteration is terminated if
    """  
    Solver for  
477    
478     M{Ax=b}     M{|r| <= atol+rtol*|r0|}
479    
480     with a symmetric and positive definite operator A (more details required!).     where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
    It uses the conjugate gradient method with preconditioner M providing an approximation of A.  
481    
482     The iteration is terminated if the C{stoppingcriterium} function return C{True}.     M{|r| = sqrt( bilinearform(Msolve(r),r))}
483    
484     For details on the preconditioned conjugate gradient method see the book:     For details on the preconditioned conjugate gradient method see the book:
485    
486     Templates for the Solution of Linear Systems by R. Barrett, M. Berry,     I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
487     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,
488     C. Romine, and H. van der Vorst.     C. Romine, and H. van der Vorst}.
489    
490     @param b: the right hand side of the liner system. C{b} is altered.     @param r: initial residual M{r=b-Ax}. C{r} is altered.
491     @type b: 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)
492     @param Aprod: returns the value Ax     @param x: an initial guess for the solution
    @type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}. The returned object needs to be of the same type like argument C{b}.  
    @param Msolve: solves Mx=r  
    @type Msolve: function C{Msolve(r)} where C{r} is of the same type like argument C{b}. The returned object needs to be of the same  
 type like argument C{x}.  
    @param bilinearform: inner product C{<x,r>}  
    @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is . The returned value is a C{float}.  
    @param stoppingcriterium: function which returns True if a stopping criterium is meet. C{stoppingcriterium} has the arguments C{norm_r}, C{r} and C{x} giving the current norm of the residual (=C{sqrt(bilinearform(Msolve(r),r)}), the current residual and the current solution approximation. C{stoppingcriterium} is called in each iteration step.  
    @type stoppingcriterium: function that returns C{True} or C{False}  
    @param x: an initial guess for the solution. If no C{x} is given 0*b is used.  
493     @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)
494     @param iter_max: maximum number of iteration steps.     @param Aprod: returns the value Ax
495       @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
496                    argument C{x}. The returned object needs to be of the same type
497                    like argument C{r}.
498       @param Msolve: solves Mx=r
499       @type Msolve: function C{Msolve(r)} where C{r} is of the same type like
500                     argument C{r}. The returned object needs to be of the same
501                     type like argument C{x}.
502       @param bilinearform: inner product C{<x,r>}
503       @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
504                           type like argument C{x} and C{r} is. The returned value
505                           is a C{float}.
506       @param atol: absolute tolerance
507       @type atol: non-negative C{float}
508       @param rtol: relative tolerance
509       @type rtol: non-negative C{float}
510       @param iter_max: maximum number of iteration steps
511     @type iter_max: C{int}     @type iter_max: C{int}
512     @return: the solution approximation and the corresponding residual     @return: the solution approximation and the corresponding residual
513     @rtype: C{tuple}     @rtype: C{tuple}
514     @warning: C{b} and C{x} are altered.     @warning: C{r} and C{x} are altered.
515     """     """
516     iter=0     iter=0
    if x==None:  
       x=0*b  
    else:  
       b += (-1)*Aprod(x)  
    r=b  
517     rhat=Msolve(r)     rhat=Msolve(r)
518     d = rhat     d = rhat
519     rhat_dot_r = bilinearform(rhat, r)     rhat_dot_r = bilinearform(rhat, r)
520     if rhat_dot_r<0: raise NegativeNorm,"negative norm."     if rhat_dot_r<0: raise NegativeNorm,"negative norm."
521       norm_r0=math.sqrt(rhat_dot_r)
522       atol2=atol+rtol*norm_r0
523       if atol2<=0:
524          raise ValueError,"Non-positive tolarance."
525       atol2=max(atol2, 100. * util.EPSILON * norm_r0)
526    
527       if verbose: print "PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_r0, atol2)
528    
529    
530     while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):     while not math.sqrt(rhat_dot_r) <= atol2:
531         iter+=1         iter+=1
532         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
533    
534         q=Aprod(d)         q=Aprod(d)
535         alpha = rhat_dot_r / bilinearform(d, q)         alpha = rhat_dot_r / bilinearform(d, q)
536         x += alpha * d         x += alpha * d
537         r += (-alpha) * q         if isinstance(q,ArithmeticTuple):
538           r += q * (-alpha)      # Doing it the other way calls the float64.__mul__ not AT.__rmul__
539           else:
540               r += (-alpha) * q
541         rhat=Msolve(r)         rhat=Msolve(r)
542         rhat_dot_r_new = bilinearform(rhat, r)         rhat_dot_r_new = bilinearform(rhat, r)
543         beta = rhat_dot_r_new / rhat_dot_r         beta = rhat_dot_r_new / rhat_dot_r
# Line 558  type like argument C{x}. Line 546  type like argument C{x}.
546    
547         rhat_dot_r = rhat_dot_r_new         rhat_dot_r = rhat_dot_r_new
548         if rhat_dot_r<0: raise NegativeNorm,"negative norm."         if rhat_dot_r<0: raise NegativeNorm,"negative norm."
549           if verbose: print "PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))
550       if verbose: print "PCG: tolerance reached after %s steps."%iter
551       return x,r,math.sqrt(rhat_dot_r)
552    
553    class Defect(object):
554        """
555        Defines a non-linear defect F(x) of a variable x.
556        """
557        def __init__(self):
558            """
559            Initializes defect.
560            """
561            self.setDerivativeIncrementLength()
562    
563        def bilinearform(self, x0, x1):
564            """
565            Returns the inner product of x0 and x1
566    
567            @param x0: value for x0
568            @param x1: value for x1
569            @return: the inner product of x0 and x1
570            @rtype: C{float}
571            """
572            return 0
573    
574        def norm(self,x):
575            """
576            Returns the norm of argument C{x}.
577    
578            @param x: a value
579            @return: norm of argument x
580            @rtype: C{float}
581            @note: by default C{sqrt(self.bilinearform(x,x)} is returned.
582            """
583            s=self.bilinearform(x,x)
584            if s<0: raise NegativeNorm,"negative norm."
585            return math.sqrt(s)
586    
587     return x,r      def eval(self,x):
588            """
589            Returns the value F of a given C{x}.
590    
591            @param x: value for which the defect C{F} is evaluated
592            @return: value of the defect at C{x}
593            """
594            return 0
595    
596        def __call__(self,x):
597            return self.eval(x)
598    
599        def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):
600            """
601            Sets the relative length of the increment used to approximate the
602            derivative of the defect. The increment is inc*norm(x)/norm(v)*v in the
603            direction of v with x as a starting point.
604    
605            @param inc: relative increment length
606            @type inc: positive C{float}
607            """
608            if inc<=0: raise ValueError,"positive increment required."
609            self.__inc=inc
610    
611  ############################      def getDerivativeIncrementLength(self):
612  # Added by Artak          """
613  #################################3          Returns the relative increment length used to approximate the
614            derivative of the defect.
615            @return: value of the defect at C{x}
616            @rtype: positive C{float}
617            """
618            return self.__inc
619    
620  #Apply a sequence of k Givens rotations, used within gmres codes      def derivative(self, F0, x0, v, v_is_normalised=True):
621  # vrot=givapp(c, s, vin, k)          """
622  def givapp(c,s,vin):          Returns the directional derivative at C{x0} in the direction of C{v}.
623      vrot=vin # warning: vin is altered!!!!  
624            @param F0: value of this defect at x0
625            @param x0: value at which derivative is calculated
626            @param v: direction
627            @param v_is_normalised: True to indicate that C{v} is nomalized
628                                    (self.norm(v)=0)
629            @return: derivative of this defect at x0 in the direction of C{v}
630            @note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is
631                   used but this method maybe overwritten to use exact evaluation.
632            """
633            normx=self.norm(x0)
634            if normx>0:
635                 epsnew = self.getDerivativeIncrementLength() * normx
636            else:
637                 epsnew = self.getDerivativeIncrementLength()
638            if not v_is_normalised:
639                normv=self.norm(v)
640                if normv<=0:
641                   return F0*0
642                else:
643                   epsnew /= normv
644            F1=self.eval(x0 + epsnew * v)
645            return (F1-F0)/epsnew
646    
647    ######################################
648    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):
649       """
650       Solves a non-linear problem M{F(x)=0} for unknown M{x} using the stopping
651       criterion:
652    
653       M{norm(F(x) <= atol + rtol * norm(F(x0)}
654    
655       where M{x0} is the initial guess.
656    
657       @param defect: object defining the function M{F}. C{defect.norm} defines the
658                      M{norm} used in the stopping criterion.
659       @type defect: L{Defect}
660       @param x: initial guess for the solution, C{x} is altered.
661       @type x: any object type allowing basic operations such as
662                C{numpy.ndarray}, L{Data}
663       @param iter_max: maximum number of iteration steps
664       @type iter_max: positive C{int}
665       @param sub_iter_max: maximum number of inner iteration steps
666       @type sub_iter_max: positive C{int}
667       @param atol: absolute tolerance for the solution
668       @type atol: positive C{float}
669       @param rtol: relative tolerance for the solution
670       @type rtol: positive C{float}
671       @param gamma: tolerance safety factor for inner iteration
672       @type gamma: positive C{float}, less than 1
673       @param sub_tol_max: upper bound for inner tolerance
674       @type sub_tol_max: positive C{float}, less than 1
675       @return: an approximation of the solution with the desired accuracy
676       @rtype: same type as the initial guess
677       """
678       lmaxit=iter_max
679       if atol<0: raise ValueError,"atol needs to be non-negative."
680       if rtol<0: raise ValueError,"rtol needs to be non-negative."
681       if rtol+atol<=0: raise ValueError,"rtol or atol needs to be non-negative."
682       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
683       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
684    
685       F=defect(x)
686       fnrm=defect.norm(F)
687       stop_tol=atol + rtol*fnrm
688       sub_tol=sub_tol_max
689       if verbose: print "NewtonGMRES: initial residual = %e."%fnrm
690       if verbose: print "             tolerance = %e."%sub_tol
691       iter=1
692       #
693       # main iteration loop
694       #
695       while not fnrm<=stop_tol:
696                if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
697                #
698            #   adjust sub_tol_
699            #
700                if iter > 1:
701               rat=fnrm/fnrmo
702                   sub_tol_old=sub_tol
703               sub_tol=gamma*rat**2
704               if gamma*sub_tol_old**2 > .1: sub_tol=max(sub_tol,gamma*sub_tol_old**2)
705               sub_tol=max(min(sub_tol,sub_tol_max), .5*stop_tol/fnrm)
706            #
707            # calculate newton increment xc
708                #     if iter_max in __FDGMRES is reached MaxIterReached is thrown
709                #     if iter_restart -1 is returned as sub_iter
710                #     if  atol is reached sub_iter returns the numer of steps performed to get there
711                #
712                #
713                if verbose: print "             subiteration (GMRES) is called with relative tolerance %e."%sub_tol
714                try:
715                   xc, sub_iter=__FDGMRES(F, defect, x, sub_tol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)
716                except MaxIterReached:
717                   raise MaxIterReached,"maximum number of %s steps reached."%iter_max
718                if sub_iter<0:
719                   iter+=sub_iter_max
720                else:
721                   iter+=sub_iter
722                # ====
723            x+=xc
724                F=defect(x)
725            iter+=1
726                fnrmo, fnrm=fnrm, defect.norm(F)
727                if verbose: print "             step %s: residual %e."%(iter,fnrm)
728       if verbose: print "NewtonGMRES: completed after %s steps."%iter
729       return x
730    
731    def __givapp(c,s,vin):
732        """
733        Applies a sequence of Givens rotations (c,s) recursively to the vector
734        C{vin}
735    
736        @warning: C{vin} is altered.
737        """
738        vrot=vin
739      if isinstance(c,float):      if isinstance(c,float):
740          vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]          vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
741      else:      else:
742          for i in range(len(c)):          for i in range(len(c)):
743              w1=c[i]*vrot[i]-s[i]*vrot[i+1]              w1=c[i]*vrot[i]-s[i]*vrot[i+1]
744          w2=s[i]*vrot[i]+c[i]*vrot[i+1]          w2=s[i]*vrot[i]+c[i]*vrot[i+1]
745              vrot[i:i+2]=w1,w2              vrot[i]=w1
746                vrot[i+1]=w2
747      return vrot      return vrot
748    
749  def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=10):  def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
750       h=numpy.zeros((iter_restart,iter_restart),numpy.float64)
751       c=numpy.zeros(iter_restart,numpy.float64)
752       s=numpy.zeros(iter_restart,numpy.float64)
753       g=numpy.zeros(iter_restart,numpy.float64)
754       v=[]
755    
756       rho=defect.norm(F0)
757       if rho<=0.: return x0*0
758    
759       v.append(-F0/rho)
760       g[0]=rho
761       iter=0
762       while rho > atol and iter<iter_restart-1:
763            if iter  >= iter_max:
764                raise MaxIterReached,"maximum number of %s steps reached."%iter_max
765    
766            p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
767            v.append(p)
768    
769            v_norm1=defect.norm(v[iter+1])
770    
771            # Modified Gram-Schmidt
772            for j in range(iter+1):
773                h[j,iter]=defect.bilinearform(v[j],v[iter+1])
774                v[iter+1]-=h[j,iter]*v[j]
775    
776            h[iter+1,iter]=defect.norm(v[iter+1])
777            v_norm2=h[iter+1,iter]
778    
779            # Reorthogonalize if needed
780            if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
781                for j in range(iter+1):
782                    hr=defect.bilinearform(v[j],v[iter+1])
783                    h[j,iter]=h[j,iter]+hr
784                    v[iter+1] -= hr*v[j]
785    
786                v_norm2=defect.norm(v[iter+1])
787                h[iter+1,iter]=v_norm2
788            #   watch out for happy breakdown
789            if not v_norm2 == 0:
790                v[iter+1]=v[iter+1]/h[iter+1,iter]
791    
792            #   Form and store the information for the new Givens rotation
793            if iter > 0 :
794                hhat=numpy.zeros(iter+1,numpy.float64)
795                for i in range(iter+1) : hhat[i]=h[i,iter]
796                hhat=__givapp(c[0:iter],s[0:iter],hhat);
797                for i in range(iter+1) : h[i,iter]=hhat[i]
798    
799            mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
800    
801            if mu!=0 :
802                c[iter]=h[iter,iter]/mu
803                s[iter]=-h[iter+1,iter]/mu
804                h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
805                h[iter+1,iter]=0.0
806                gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
807                g[iter]=gg[0]
808                g[iter+1]=gg[1]
809    
810            # Update the residual norm
811            rho=abs(g[iter+1])
812            iter+=1
813    
814       # At this point either iter > iter_max or rho < tol.
815       # It's time to compute x and leave.
816       if iter > 0 :
817         y=numpy.zeros(iter,numpy.float64)
818         y[iter-1] = g[iter-1] / h[iter-1,iter-1]
819         if iter > 1 :
820            i=iter-2
821            while i>=0 :
822              y[i] = ( g[i] - numpy.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
823              i=i-1
824         xhat=v[iter-1]*y[iter-1]
825         for i in range(iter-1):
826        xhat += v[i]*y[i]
827       else :
828          xhat=v[0] * 0
829    
830       if iter<iter_restart-1:
831          stopped=iter
832       else:
833          stopped=-1
834    
835       return xhat,stopped
836    
837    def GMRES(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100, iter_restart=20, verbose=False,P_R=None):
838       """
839       Solver for
840    
841       M{Ax=b}
842    
843       with a general operator A (more details required!).
844       It uses the generalized minimum residual method (GMRES).
845    
846       The iteration is terminated if
847    
848       M{|r| <= atol+rtol*|r0|}
849    
850       where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
851    
852       M{|r| = sqrt( bilinearform(r,r))}
853    
854       @param r: initial residual M{r=b-Ax}. C{r} is altered.
855       @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
856       @param x: an initial guess for the solution
857       @type x: same like C{r}
858       @param Aprod: returns the value Ax
859       @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
860                    argument C{x}. The returned object needs to be of the same
861                    type like argument C{r}.
862       @param bilinearform: inner product C{<x,r>}
863       @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
864                           type like argument C{x} and C{r}. The returned value is
865                           a C{float}.
866       @param atol: absolute tolerance
867       @type atol: non-negative C{float}
868       @param rtol: relative tolerance
869       @type rtol: non-negative C{float}
870       @param iter_max: maximum number of iteration steps
871       @type iter_max: C{int}
872       @param iter_restart: in order to save memory the orthogonalization process
873                            is terminated after C{iter_restart} steps and the
874                            iteration is restarted.
875       @type iter_restart: C{int}
876       @return: the solution approximation and the corresponding residual
877       @rtype: C{tuple}
878       @warning: C{r} and C{x} are altered.
879       """
880     m=iter_restart     m=iter_restart
881       restarted=False
882     iter=0     iter=0
883       if rtol>0:
884          r_dot_r = bilinearform(r, r)
885          if r_dot_r<0: raise NegativeNorm,"negative norm."
886          atol2=atol+rtol*math.sqrt(r_dot_r)
887          if verbose: print "GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)
888       else:
889          atol2=atol
890          if verbose: print "GMRES: absolute tolerance = %e"%atol2
891       if atol2<=0:
892          raise ValueError,"Non-positive tolarance."
893    
894     while True:     while True:
895        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
896        x,stopped=GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=x, iter_max=iter_max-iter, iter_restart=m)        if restarted:
897             r2 = r-Aprod(x-x2)
898          else:
899             r2=1*r
900          x2=x*1.
901          x,stopped=_GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose,P_R=P_R)
902        iter+=iter_restart        iter+=iter_restart
903        if stopped: break        if stopped: break
904          if verbose: print "GMRES: restart."
905          restarted=True
906       if verbose: print "GMRES: tolerance has been reached."
907     return x     return x
908    
909  def GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=10):  def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False, P_R=None):
910     iter=0     iter=0
    r=Msolve(b)  
    r_dot_r = bilinearform(r, r)  
    if r_dot_r<0: raise NegativeNorm,"negative norm."  
    norm_b=math.sqrt(r_dot_r)  
911    
912     if x==None:     h=numpy.zeros((iter_restart+1,iter_restart),numpy.float64)
913        x=0*b     c=numpy.zeros(iter_restart,numpy.float64)
914     else:     s=numpy.zeros(iter_restart,numpy.float64)
915        r=Msolve(b-Aprod(x))     g=numpy.zeros(iter_restart+1,numpy.float64)
       r_dot_r = bilinearform(r, r)  
       if r_dot_r<0: raise NegativeNorm,"negative norm."  
     
    h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)  
    c=numarray.zeros(iter_restart,numarray.Float64)  
    s=numarray.zeros(iter_restart,numarray.Float64)  
    g=numarray.zeros(iter_restart,numarray.Float64)  
916     v=[]     v=[]
917    
918       r_dot_r = bilinearform(r, r)
919       if r_dot_r<0: raise NegativeNorm,"negative norm."
920     rho=math.sqrt(r_dot_r)     rho=math.sqrt(r_dot_r)
921    
922     v.append(r/rho)     v.append(r/rho)
923     g[0]=rho     g[0]=rho
924    
925     while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):     if verbose: print "GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)
926       while not (rho<=atol or iter==iter_restart):
927    
928      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
929    
930                if P_R!=None:
931      p=Msolve(Aprod(v[iter]))              p=Aprod(P_R(v[iter]))
932            else:
933            p=Aprod(v[iter])
934      v.append(p)      v.append(p)
935    
936      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
937    
938  # Modified Gram-Schmidt  # Modified Gram-Schmidt
939      for j in range(iter+1):      for j in range(iter+1):
940        h[j][iter]=bilinearform(v[j],v[iter+1])          h[j,iter]=bilinearform(v[j],v[iter+1])
941        v[iter+1]+=(-1.)*h[j][iter]*v[j]        v[iter+1]-=h[j,iter]*v[j]
         
     h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))  
     v_norm2=h[iter+1][iter]  
942    
943        h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
944        v_norm2=h[iter+1,iter]
945    
946  # Reorthogonalize if needed  # Reorthogonalize if needed
947      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)
948       for j in range(iter+1):       for j in range(iter+1):
949          hr=bilinearform(v[j],v[iter+1])          hr=bilinearform(v[j],v[iter+1])
950              h[j][iter]=h[j][iter]+hr #vhat              h[j,iter]=h[j,iter]+hr
951              v[iter+1] +=(-1.)*hr*v[j]              v[iter+1] -= hr*v[j]
952    
953       v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
954       h[iter+1][iter]=v_norm2       h[iter+1,iter]=v_norm2
955    
956  #   watch out for happy breakdown  #   watch out for happy breakdown
957          if v_norm2 != 0:          if not v_norm2 == 0:
958           v[iter+1]=v[iter+1]/h[iter+1][iter]           v[iter+1]=v[iter+1]/h[iter+1,iter]
959    
960  #   Form and store the information for the new Givens rotation  #   Form and store the information for the new Givens rotation
961      if iter > 0 :      if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])
962          hhat=[]      mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
         for i in range(iter+1) : hhat.append(h[i][iter])  
         hhat=givapp(c[0:iter],s[0:iter],hhat);  
             for i in range(iter+1) : h[i][iter]=hhat[i]  
963    
     mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])  
964      if mu!=0 :      if mu!=0 :
965          c[iter]=h[iter][iter]/mu          c[iter]=h[iter,iter]/mu
966          s[iter]=-h[iter+1][iter]/mu          s[iter]=-h[iter+1,iter]/mu
967          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]
968          h[iter+1][iter]=0.0          h[iter+1,iter]=0.0
969          g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])                  gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
970                    g[iter]=gg[0]
971                    g[iter+1]=gg[1]
972  # Update the residual norm  # Update the residual norm
973    
974          rho=abs(g[iter+1])          rho=abs(g[iter+1])
975            if verbose: print "GMRES: iteration step %s: residual %e"%(iter,rho)
976      iter+=1      iter+=1
977    
978  # At this point either iter > iter_max or rho < tol.  # At this point either iter > iter_max or rho < tol.
979  # It's time to compute x and leave.          # It's time to compute x and leave.
980    
981     if iter > 0 :     if verbose: print "GMRES: iteration stopped after %s step."%iter
982       y=numarray.zeros(iter,numarray.Float64)         if iter > 0 :
983       y[iter-1] = g[iter-1] / h[iter-1][iter-1]       y=numpy.zeros(iter,numpy.float64)
984       if iter > 1 :         y[iter-1] = g[iter-1] / h[iter-1,iter-1]
985          i=iter-2         if iter > 1 :
986            i=iter-2
987          while i>=0 :          while i>=0 :
988            y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]            y[i] = ( g[i] - numpy.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
989            i=i-1            i=i-1
990       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
991       for i in range(iter-1):       for i in range(iter-1):
992      xhat += v[i]*y[i]      xhat += v[i]*y[i]
993     else : xhat=v[0]     else:
994             xhat=v[0] * 0
995     x += xhat     if P_R!=None:
996     if iter!=iter_restart-1:        x += P_R(xhat)
997        stopped=True     else:
998     else:        x += xhat
999       if iter<iter_restart-1:
1000          stopped=True
1001       else:
1002        stopped=False        stopped=False
1003    
1004     return x,stopped     return x,stopped
1005    
1006  def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  def MINRES(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1007        """
1008      #      Solver for
     #  minres solves the system of linear equations Ax = b  
     #  where A is a symmetric matrix (possibly indefinite or singular)  
     #  and b is a given vector.  
     #    
     #  "A" may be a dense or sparse matrix (preferably sparse!)  
     #  or the name of a function such that  
     #               y = A(x)  
     #  returns the product y = Ax for any given vector x.  
     #  
     #  "M" defines a positive-definite preconditioner M = C C'.  
     #  "M" may be a dense or sparse matrix (preferably sparse!)  
     #  or the name of a function such that  
     #  solves the system My = x for any given vector x.  
     #  
     #  
1009    
1010      #  Initialize                                    M{Ax=b}
1011    
1012      iter   = 0      with a symmetric and positive definite operator A (more details required!).
1013      Anorm = 0      It uses the minimum residual method (MINRES) with preconditioner M
1014      ynorm = 0      providing an approximation of A.
1015      x=x*0  
1016        The iteration is terminated if
1017    
1018        M{|r| <= atol+rtol*|r0|}
1019    
1020        where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
1021    
1022        M{|r| = sqrt( bilinearform(Msolve(r),r))}
1023    
1024        For details on the preconditioned conjugate gradient method see the book:
1025    
1026        I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
1027        T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
1028        C. Romine, and H. van der Vorst}.
1029    
1030        @param r: initial residual M{r=b-Ax}. C{r} is altered.
1031        @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1032        @param x: an initial guess for the solution
1033        @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1034        @param Aprod: returns the value Ax
1035        @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
1036                     argument C{x}. The returned object needs to be of the same
1037                     type like argument C{r}.
1038        @param Msolve: solves Mx=r
1039        @type Msolve: function C{Msolve(r)} where C{r} is of the same type like
1040                      argument C{r}. The returned object needs to be of the same
1041                      type like argument C{x}.
1042        @param bilinearform: inner product C{<x,r>}
1043        @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
1044                            type like argument C{x} and C{r} is. The returned value
1045                            is a C{float}.
1046        @param atol: absolute tolerance
1047        @type atol: non-negative C{float}
1048        @param rtol: relative tolerance
1049        @type rtol: non-negative C{float}
1050        @param iter_max: maximum number of iteration steps
1051        @type iter_max: C{int}
1052        @return: the solution approximation and the corresponding residual
1053        @rtype: C{tuple}
1054        @warning: C{r} and C{x} are altered.
1055        """
1056      #------------------------------------------------------------------      #------------------------------------------------------------------
1057      # Set up y and v for the first Lanczos vector v1.      # Set up y and v for the first Lanczos vector v1.
1058      # y  =  beta1 P' v1,  where  P = C**(-1).      # y  =  beta1 P' v1,  where  P = C**(-1).
1059      # v is really P' v1.      # v is really P' v1.
1060      #------------------------------------------------------------------      #------------------------------------------------------------------
1061      r1    = b      r1    = r
1062      y = Msolve(b)      y = Msolve(r)
1063      beta1 = bilinearform(b,y)      beta1 = bilinearform(y,r)
1064    
1065      if beta1< 0: raise NegativeNorm,"negative norm."      if beta1< 0: raise NegativeNorm,"negative norm."
1066    
1067      #  If b = 0 exactly, stop with x = 0.      #  If r = 0 exactly, stop with x
1068      if beta1==0: return x*0.      if beta1==0: return x
1069    
1070      if beta1> 0:      if beta1> 0: beta1  = math.sqrt(beta1)
       beta1  = math.sqrt(beta1)       # Normalize y to get v1 later.  
1071    
1072      #------------------------------------------------------------------      #------------------------------------------------------------------
1073      # Initialize other quantities.      # Initialize quantities.
1074      # ------------------------------------------------------------------      # ------------------------------------------------------------------
1075        iter   = 0
1076        Anorm = 0
1077        ynorm = 0
1078      oldb   = 0      oldb   = 0
1079      beta   = beta1      beta   = beta1
1080      dbar   = 0      dbar   = 0
# Line 747  def MINRES(b, Aprod, Msolve, bilinearfor Line 1087  def MINRES(b, Aprod, Msolve, bilinearfor
1087      ynorm2 = 0      ynorm2 = 0
1088      cs     = -1      cs     = -1
1089      sn     = 0      sn     = 0
1090      w      = b*0.      w      = r*0.
1091      w2     = b*0.      w2     = r*0.
1092      r2     = r1      r2     = r1
1093      eps    = 0.0001      eps    = 0.0001
1094    
1095      #---------------------------------------------------------------------      #---------------------------------------------------------------------
1096      # Main iteration loop.      # Main iteration loop.
1097      # --------------------------------------------------------------------      # --------------------------------------------------------------------
1098      while not stoppingcriterium(rnorm,Anorm*ynorm):    #  ||r|| / (||A|| ||x||)      while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL
1099    
1100      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
1101          iter    = iter  +  1          iter    = iter  +  1
# Line 774  def MINRES(b, Aprod, Msolve, bilinearfor Line 1114  def MINRES(b, Aprod, Msolve, bilinearfor
1114          #-----------------------------------------------------------------          #-----------------------------------------------------------------
1115          s = 1/beta                 # Normalize previous vector (in y).          s = 1/beta                 # Normalize previous vector (in y).
1116          v = s*y                    # v = vk if P = I          v = s*y                    # v = vk if P = I
1117        
1118          y      = Aprod(v)          y      = Aprod(v)
1119        
1120          if iter >= 2:          if iter >= 2:
1121            y = y - (beta/oldb)*r1            y = y - (beta/oldb)*r1
1122    
1123          alfa   = bilinearform(v,y)              # alphak          alfa   = bilinearform(v,y)              # alphak
1124          y      = (- alfa/beta)*r2 + y          y      += (- alfa/beta)*r2
1125          r1     = r2          r1     = r2
1126          r2     = y          r2     = y
1127          y = Msolve(r2)          y = Msolve(r2)
1128          oldb   = beta                           # oldb = betak          oldb   = beta                           # oldb = betak
1129          beta   = bilinearform(r2,y)             # beta = betak+1^2          beta   = bilinearform(y,r2)             # beta = betak+1^2
1130          if beta < 0: raise NegativeNorm,"negative norm."          if beta < 0: raise NegativeNorm,"negative norm."
1131    
1132          beta   = math.sqrt( beta )          beta   = math.sqrt( beta )
1133          tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta          tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
1134            
1135          if iter==1:                 # Initialize a few things.          if iter==1:                 # Initialize a few things.
1136            gmax   = abs( alfa )      # alpha1            gmax   = abs( alfa )      # alpha1
1137            gmin   = gmax             # alpha1            gmin   = gmax             # alpha1
# Line 799  def MINRES(b, Aprod, Msolve, bilinearfor Line 1139  def MINRES(b, Aprod, Msolve, bilinearfor
1139          # Apply previous rotation Qk-1 to get          # Apply previous rotation Qk-1 to get
1140          #   [deltak epslnk+1] = [cs  sn][dbark    0   ]          #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
1141          #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].          #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
1142        
1143          oldeps = epsln          oldeps = epsln
1144          delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak          delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
1145          gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k          gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
# Line 809  def MINRES(b, Aprod, Msolve, bilinearfor Line 1149  def MINRES(b, Aprod, Msolve, bilinearfor
1149          # Compute the next plane rotation Qk          # Compute the next plane rotation Qk
1150    
1151          gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak          gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
1152          gamma  = max(gamma,eps)          gamma  = max(gamma,eps)
1153          cs     = gbar / gamma             # ck          cs     = gbar / gamma             # ck
1154          sn     = beta / gamma             # sk          sn     = beta / gamma             # sk
1155          phi    = cs * phibar              # phik          phi    = cs * phibar              # phik
# Line 817  def MINRES(b, Aprod, Msolve, bilinearfor Line 1157  def MINRES(b, Aprod, Msolve, bilinearfor
1157    
1158          # Update  x.          # Update  x.
1159    
1160          denom = 1/gamma          denom = 1/gamma
1161          w1    = w2          w1    = w2
1162          w2    = w          w2    = w
1163          w     = (v - oldeps*w1 - delta*w2) * denom          w     = (v - oldeps*w1 - delta*w2) * denom
1164          x     = x  +  phi*w          x     +=  phi*w
1165    
1166          # Go round again.          # Go round again.
1167    
# Line 834  def MINRES(b, Aprod, Msolve, bilinearfor Line 1174  def MINRES(b, Aprod, Msolve, bilinearfor
1174    
1175          # Estimate various norms and test for convergence.          # Estimate various norms and test for convergence.
1176    
1177          Anorm  = math.sqrt( tnorm2 )          Anorm  = math.sqrt( tnorm2 )
1178          ynorm  = math.sqrt( ynorm2 )          ynorm  = math.sqrt( ynorm2 )
1179    
1180          rnorm  = phibar          rnorm  = phibar
1181    
     # Return final answer.  
1182      return x      return x
1183        
1184    def TFQMR(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1185      """
1186      Solver for
1187    
1188      M{Ax=b}
1189    
1190      with a general operator A (more details required!).
1191      It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).
1192    
1193      The iteration is terminated if
1194    
1195      M{|r| <= atol+rtol*|r0|}
1196    
1197      where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
1198    
1199      M{|r| = sqrt( bilinearform(r,r))}
1200    
1201      @param r: initial residual M{r=b-Ax}. C{r} is altered.
1202      @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1203      @param x: an initial guess for the solution
1204      @type x: same like C{r}
1205      @param Aprod: returns the value Ax
1206      @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
1207                   argument C{x}. The returned object needs to be of the same type
1208                   like argument C{r}.
1209      @param bilinearform: inner product C{<x,r>}
1210      @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
1211                          type like argument C{x} and C{r}. The returned value is
1212                          a C{float}.
1213      @param atol: absolute tolerance
1214      @type atol: non-negative C{float}
1215      @param rtol: relative tolerance
1216      @type rtol: non-negative C{float}
1217      @param iter_max: maximum number of iteration steps
1218      @type iter_max: C{int}
1219      @rtype: C{tuple}
1220      @warning: C{r} and C{x} are altered.
1221      """
1222      u1=0
1223      u2=0
1224      y1=0
1225      y2=0
1226    
1227      w = r
1228      y1 = r
1229      iter = 0
1230      d = 0
1231      v = Aprod(y1)
1232      u1 = v
1233    
1234      theta = 0.0;
1235      eta = 0.0;
1236      rho=bilinearform(r,r)
1237      if rho < 0: raise NegativeNorm,"negative norm."
1238      tau = math.sqrt(rho)
1239      norm_r0=tau
1240      while tau>atol+rtol*norm_r0:
1241        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1242    
1243        sigma = bilinearform(r,v)
1244        if sigma == 0.0: raise IterationBreakDown,'TFQMR breakdown, sigma=0'
1245    
1246        alpha = rho / sigma
1247    
1248        for j in range(2):
1249    #
1250    #   Compute y2 and u2 only if you have to
1251    #
1252          if ( j == 1 ):
1253            y2 = y1 - alpha * v
1254            u2 = Aprod(y2)
1255    
1256          m = 2 * (iter+1) - 2 + (j+1)
1257          if j==0:
1258             w = w - alpha * u1
1259             d = y1 + ( theta * theta * eta / alpha ) * d
1260          if j==1:
1261             w = w - alpha * u2
1262             d = y2 + ( theta * theta * eta / alpha ) * d
1263    
1264          theta = math.sqrt(bilinearform(w,w))/ tau
1265          c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1266          tau = tau * theta * c
1267          eta = c * c * alpha
1268          x = x + eta * d
1269    #
1270    #  Try to terminate the iteration at each pass through the loop
1271    #
1272        if rho == 0.0: raise IterationBreakDown,'TFQMR breakdown, rho=0'
1273    
1274        rhon = bilinearform(r,w)
1275        beta = rhon / rho;
1276        rho = rhon;
1277        y1 = w + beta * y2;
1278        u1 = Aprod(y1)
1279        v = u1 + beta * ( u2 + beta * v )
1280    
1281        iter += 1
1282    
1283      return x
1284    
1285    
1286  #############################################  #############################################
1287    
1288  class ArithmeticTuple(object):  class ArithmeticTuple(object):
1289     """     """
1290     tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.     Tuple supporting inplace update x+=y and scaling x=a*y where C{x,y} is an
1291       ArithmeticTuple and C{a} is a float.
1292    
1293     example of usage:     Example of usage::
1294    
1295     from esys.escript import Data         from esys.escript import Data
1296     from numarray import array         from numpy import array
1297     a=Data(...)         a=Data(...)
1298     b=array([1.,4.])         b=array([1.,4.])
1299     x=ArithmeticTuple(a,b)         x=ArithmeticTuple(a,b)
1300     y=5.*x         y=5.*x
1301    
1302     """     """
1303     def __init__(self,*args):     def __init__(self,*args):
1304         """         """
1305         initialize object with elements args.         Initializes object with elements C{args}.
1306    
1307         @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)         @param args: tuple of objects that support inplace add (x+=y) and
1308                        scaling (x=a*y)
1309         """         """
1310         self.__items=list(args)         self.__items=list(args)
1311    
1312     def __len__(self):     def __len__(self):
1313         """         """
1314         number of items         Returns the number of items.
1315    
1316         @return: number of items         @return: number of items
1317         @rtype: C{int}         @rtype: C{int}
# Line 877  class ArithmeticTuple(object): Line 1320  class ArithmeticTuple(object):
1320    
1321     def __getitem__(self,index):     def __getitem__(self,index):
1322         """         """
1323         get an item         Returns item at specified position.
1324    
1325         @param index: item to be returned         @param index: index of item to be returned
1326         @type index: C{int}         @type index: C{int}
1327         @return: item with index C{index}         @return: item with index C{index}
1328         """         """
# Line 887  class ArithmeticTuple(object): Line 1330  class ArithmeticTuple(object):
1330    
1331     def __mul__(self,other):     def __mul__(self,other):
1332         """         """
1333         scaling from the right         Scales by C{other} from the right.
1334    
1335         @param other: scaling factor         @param other: scaling factor
1336         @type other: C{float}         @type other: C{float}
# Line 895  class ArithmeticTuple(object): Line 1338  class ArithmeticTuple(object):
1338         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1339         """         """
1340         out=[]         out=[]
1341         for i in range(len(self)):         try:
1342             out.append(self[i]*other)             l=len(other)
1343               if l!=len(self):
1344                   raise ValueError,"length of arguments don't match."
1345               for i in range(l): out.append(self[i]*other[i])
1346           except TypeError:
1347               for i in range(len(self)): out.append(self[i]*other)
1348         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1349    
1350     def __rmul__(self,other):     def __rmul__(self,other):
1351         """         """
1352         scaling from the left         Scales by C{other} from the left.
1353    
1354         @param other: scaling factor         @param other: scaling factor
1355         @type other: C{float}         @type other: C{float}
# Line 909  class ArithmeticTuple(object): Line 1357  class ArithmeticTuple(object):
1357         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1358         """         """
1359         out=[]         out=[]
1360         for i in range(len(self)):         try:
1361             out.append(other*self[i])             l=len(other)
1362               if l!=len(self):
1363                   raise ValueError,"length of arguments don't match."
1364               for i in range(l): out.append(other[i]*self[i])
1365           except TypeError:
1366               for i in range(len(self)): out.append(other*self[i])
1367         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1368    
 #########################  
 # Added by Artak  
 #########################  
1369     def __div__(self,other):     def __div__(self,other):
1370         """         """
1371         dividing from the right         Scales by (1/C{other}) from the right.
1372    
1373         @param other: scaling factor         @param other: scaling factor
1374         @type other: C{float}         @type other: C{float}
1375         @return: itemwise self/other         @return: itemwise self/other
1376         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1377         """         """
1378         out=[]         return self*(1/other)
        for i in range(len(self)):  
            out.append(self[i]/other)  
        return ArithmeticTuple(*tuple(out))  
1379    
1380     def __rdiv__(self,other):     def __rdiv__(self,other):
1381         """         """
1382         dividing from the left         Scales by (1/C{other}) from the left.
1383    
1384         @param other: scaling factor         @param other: scaling factor
1385         @type other: C{float}         @type other: C{float}
# Line 940  class ArithmeticTuple(object): Line 1387  class ArithmeticTuple(object):
1387         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1388         """         """
1389         out=[]         out=[]
1390         for i in range(len(self)):         try:
1391             out.append(other/self[i])             l=len(other)
1392               if l!=len(self):
1393                   raise ValueError,"length of arguments don't match."
1394               for i in range(l): out.append(other[i]/self[i])
1395           except TypeError:
1396               for i in range(len(self)): out.append(other/self[i])
1397         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
     
 ##########################################33  
1398    
1399     def __iadd__(self,other):     def __iadd__(self,other):
1400         """         """
1401         in-place add of other to self         Inplace addition of C{other} to self.
1402    
1403         @param other: increment         @param other: increment
1404         @type other: C{ArithmeticTuple}         @type other: C{ArithmeticTuple}
1405         """         """
1406         if len(self) != len(other):         if len(self) != len(other):
1407             raise ValueError,"tuple length must match."             raise ValueError,"tuple lengths must match."
1408         for i in range(len(self)):         for i in range(len(self)):
1409             self.__items[i]+=other[i]             self.__items[i]+=other[i]
1410         return self         return self
1411    
1412       def __add__(self,other):
1413           """
1414           Adds C{other} to self.
1415    
1416           @param other: increment
1417           @type other: C{ArithmeticTuple}
1418           """
1419           out=[]
1420           try:
1421               l=len(other)
1422               if l!=len(self):
1423                   raise ValueError,"length of arguments don't match."
1424               for i in range(l): out.append(self[i]+other[i])
1425           except TypeError:
1426               for i in range(len(self)): out.append(self[i]+other)
1427           return ArithmeticTuple(*tuple(out))
1428    
1429       def __sub__(self,other):
1430           """
1431           Subtracts C{other} from self.
1432    
1433           @param other: decrement
1434           @type other: C{ArithmeticTuple}
1435           """
1436           out=[]
1437           try:
1438               l=len(other)
1439               if l!=len(self):
1440                   raise ValueError,"length of arguments don't match."
1441               for i in range(l): out.append(self[i]-other[i])
1442           except TypeError:
1443               for i in range(len(self)): out.append(self[i]-other)
1444           return ArithmeticTuple(*tuple(out))
1445    
1446       def __isub__(self,other):
1447           """
1448           Inplace subtraction of C{other} from self.
1449    
1450           @param other: decrement
1451           @type other: C{ArithmeticTuple}
1452           """
1453           if len(self) != len(other):
1454               raise ValueError,"tuple length must match."
1455           for i in range(len(self)):
1456               self.__items[i]-=other[i]
1457           return self
1458    
1459       def __neg__(self):
1460           """
1461           Negates values.
1462           """
1463           out=[]
1464           for i in range(len(self)):
1465               out.append(-self[i])
1466           return ArithmeticTuple(*tuple(out))
1467    
1468    
1469  class HomogeneousSaddlePointProblem(object):  class HomogeneousSaddlePointProblem(object):
1470        """        """
1471        This provides a framwork for solving homogeneous saddle point problem of the form        This class provides a framework for solving linear homogeneous saddle
1472          point problems of the form::
              Av+B^*p=f  
              Bv    =0  
1473    
1474        for the unknowns v and p and given operators A and B and given right hand side f.            M{Av+B^*p=f}
1475        B^* is the adjoint operator of B is the given inner product.            M{Bv     =0}
1476    
1477          for the unknowns M{v} and M{p} and given operators M{A} and M{B} and
1478          given right hand side M{f}. M{B^*} is the adjoint operator of M{B}.
1479        """        """
1480        def __init__(self,**kwargs):        def __init__(self, adaptSubTolerance=True, **kwargs):
1481        """
1482        initializes the saddle point problem
1483        
1484        @param adaptSubTolerance: If True the tolerance for subproblem is set automatically.
1485        @type adaptSubTolerance: C{bool}
1486        """
1487          self.setTolerance()          self.setTolerance()
1488          self.setToleranceReductionFactor()          self.setAbsoluteTolerance()
1489        self.__adaptSubTolerance=adaptSubTolerance
1490          #=============================================================
1491        def initialize(self):        def initialize(self):
1492          """          """
1493          initialize the problem (overwrite)          Initializes the problem (overwrite).
1494          """          """
1495          pass          pass
       def B(self,v):  
          """  
          returns Bv (overwrite)  
          @rtype: equal to the type of p  
1496    
1497           @note: boundary conditions on p should be zero!        def inner_pBv(self,p,Bv):
1498           """           """
1499           pass           Returns inner product of element p and Bv (overwrite).
1500    
1501        def inner(self,p0,p1):           @param p: a pressure increment
1502           """           @param v: a residual
1503           returns inner product of two element p0 and p1  (overwrite)           @return: inner product of element p and Bv
           
          @type p0: equal to the type of p  
          @type p1: equal to the type of p  
1504           @rtype: C{float}           @rtype: C{float}
1505             @note: used if PCG is applied.
          @rtype: equal to the type of p  
1506           """           """
1507           pass           raise NotImplementedError,"no inner product for p and Bv implemented."
1508    
1509        def solve_A(self,u,p):        def inner_p(self,p0,p1):
1510           """           """
1511           solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)           Returns inner product of p0 and p1 (overwrite).
1512    
1513           @rtype: equal to the type of v           @param p0: a pressure
1514           @note: boundary conditions on v should be zero!           @param p1: a pressure
1515             @return: inner product of p0 and p1
1516             @rtype: C{float}
1517           """           """
1518           pass           raise NotImplementedError,"no inner product for p implemented."
1519      
1520        def solve_prec(self,p):        def norm_v(self,v):
1521           """           """
1522           provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)           Returns the norm of v (overwrite).
1523    
1524           @rtype: equal to the type of p           @param v: a velovity
1525             @return: norm of v
1526             @rtype: non-negative C{float}
1527           """           """
1528           pass           raise NotImplementedError,"no norm of v implemented."
1529          def getV(self, p, v0):
       def stoppingcriterium(self,Bv,v,p):  
1530           """           """
1531           returns a True if iteration is terminated. (overwrite)           return the value for v for a given p (overwrite)
1532    
1533           @rtype: C{bool}           @param p: a pressure
1534             @param v0: a initial guess for the value v to return.
1535             @return: v given as M{v= A^{-1} (f-B^*p)}
1536           """           """
1537           pass           raise NotImplementedError,"no v calculation implemented."
               
       def __inner(self,p,r):  
          return self.inner(p,r[1])  
1538    
1539        def __inner_p(self,p1,p2):          
1540           return self.inner(p1,p2)        def Bv(self,v):
1541            """
1542            Returns Bv (overwrite).
1543    
1544        def __stoppingcriterium(self,norm_r,r,p):          @rtype: equal to the type of p
1545            return self.stoppingcriterium(r[1],r[0],p)          @note: boundary conditions on p should be zero!
1546            """
1547            raise NotImplementedError, "no operator B implemented."
1548    
1549        def __stoppingcriterium_GMRES(self,norm_r,norm_b):        def norm_Bv(self,Bv):
1550            return self.stoppingcriterium_GMRES(norm_r,norm_b)          """
1551            Returns the norm of Bv (overwrite).
1552    
1553        def __stoppingcriterium_MINRES(self,norm_r,norm_Ax):          @rtype: equal to the type of p
1554            return self.stoppingcriterium_MINRES(norm_r,norm_Ax)          @note: boundary conditions on p should be zero!
1555            """
1556            raise NotImplementedError, "no norm of Bv implemented."
1557    
1558          def solve_AinvBt(self,p):
1559             """
1560             Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
1561             (overwrite).
1562    
1563        def setTolerance(self,tolerance=1.e-8):           @param p: a pressure increment
1564                self.__tol=tolerance           @return: the solution of M{Av=B^*p}
1565        def getTolerance(self):           @note: boundary conditions on v should be zero!
1566                return self.__tol           """
1567        def setToleranceReductionFactor(self,reduction=0.01):           raise NotImplementedError,"no operator A implemented."
               self.__reduction=reduction  
       def getSubProblemTolerance(self):  
               return self.__reduction*self.getTolerance()  
1568    
1569        def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG'):        def solve_prec(self,Bv):
1570                """           """
1571                solves the saddle point problem using initial guesses v and p.           Provides a preconditioner for M{BA^{-1}B^*} applied to Bv with accuracy
1572             L{self.getSubProblemTolerance()} (overwrite).
               @param max_iter: maximum number of iteration steps.  
               """  
               self.verbose=verbose  
               self.show_details=show_details and self.verbose  
   
               # assume p is known: then v=A^-1(f-B^*p)  
               # which leads to BA^-1B^*p = BA^-1f    
   
           # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)        
   
             
           self.__z=v+self.solve_A(v,p*0)  
   
               Bz=self.B(self.__z)  
               #  
           #   solve BA^-1B^*p = Bz  
               #  
               #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv  
               #  
               #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)  
               #                           A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)  
               #  
               self.iter=0  
           if solver=='GMRES':        
                 if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter  
                 p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1.)  
                 # solve Au=f-B^*p  
                 #       A(u-v)=f-B^*p-Av  
                 #       u=v+(u-v)  
         u=v+self.solve_A(v,p)  
   
           if solver=='MINRES':        
                 if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter  
                 p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_MINRES,iter_max=max_iter, x=p*1.)  
                 # solve Au=f-B^*p  
                 #       A(u-v)=f-B^*p-Av  
                 #       u=v+(u-v)  
         u=v+self.solve_A(v,p)  
       
               if solver=='PCG':  
                 if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter  
                 p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)  
             u=r[0]    
1573    
1574                print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)           @rtype: equal to the type of p
1575             @note: boundary conditions on p should be zero!
1576             """
1577             raise NotImplementedError,"no preconditioner for Schur complement implemented."
1578          def setSubProblemTolerance(self):
1579             """
1580         Updates the tolerance for subproblems
1581         @note: method is typically the method is overwritten.
1582             """
1583             pass
1584          #=============================================================
1585          def __Aprod_PCG(self,p):
1586              dv=self.solve_AinvBt(p)
1587              return ArithmeticTuple(dv,self.Bv(dv))
1588    
1589            return u,p        def __inner_PCG(self,p,r):
1590             return self.inner_pBv(p,r[1])
1591    
1592        def __Msolve(self,r):        def __Msolve_PCG(self,r):
1593            return self.solve_prec(r[1])            return self.solve_prec(r[1])
1594          #=============================================================
       def __Msolve_GMRES(self,r):  
           return self.solve_prec(r)  
   
   
       def __Aprod(self,p):  
           # return BA^-1B*p  
           #solve Av =-B^*p as Av =f-Az-B^*p  
           v=self.solve_A(self.__z,-p)  
           return ArithmeticTuple(v, self.B(v))  
   
1595        def __Aprod_GMRES(self,p):        def __Aprod_GMRES(self,p):
1596            # return BA^-1B*p            return self.solve_prec(self.Bv(self.solve_AinvBt(p)))
1597            #solve Av =-B^*p as Av =f-Az-B^*p        def __inner_GMRES(self,p0,p1):
1598        v=self.solve_A(self.__z,-p)           return self.inner_p(p0,p1)
1599            return self.B(v)  
1600          #=============================================================
1601  class SaddlePointProblem(object):        def norm_p(self,p):
1602     """            """
1603     This implements a solver for a saddlepoint problem            calculates the norm of C{p}
1604              
1605     M{f(u,p)=0}            @param p: a pressure
1606     M{g(u)=0}            @return: the norm of C{p} using the inner product for pressure
1607              @rtype: C{float}
1608     for u and p. The problem is solved with an inexact Uszawa scheme for p:            """
1609              f=self.inner_p(p,p)
1610     M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}            if f<0: raise ValueError,"negative pressure norm."
1611     M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}            return math.sqrt(f)
1612          def adaptSubTolerance(self):
1613          """
1614          Returns True if tolerance adaption for subproblem is choosen.
1615          """
1616              self.__adaptSubTolerance
1617          
1618          def solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1619             """
1620             Solves the saddle point problem using initial guesses v and p.
1621    
1622     where Q_f is an approximation of the Jacobiean A_f of f with respect to u  and Q_f is an approximation of           @param v: initial guess for velocity
1623     A_g A_f^{-1} A_g with A_g is the jacobiean of g with respect to p. As a the construction of a 'proper'           @param p: initial guess for pressure
1624     Q_g can be difficult, non-linear conjugate gradient method is applied to solve for p, so Q_g plays           @type v: L{Data}
1625     in fact the role of a preconditioner.           @type p: L{Data}
1626     """           @param usePCG: indicates the usage of the PCG rather than GMRES scheme.
1627     def __init__(self,verbose=False,*args):           @param max_iter: maximum number of iteration steps per correction
1628         """                            attempt
1629         initializes the problem           @param verbose: if True, shows information on the progress of the
1630                             saddlepoint problem solver.
1631             @param iter_restart: restart the iteration after C{iter_restart} steps
1632                                  (only used if useUzaw=False)
1633             @type usePCG: C{bool}
1634             @type max_iter: C{int}
1635             @type verbose: C{bool}
1636             @type iter_restart: C{int}
1637             @rtype: C{tuple} of L{Data} objects
1638             """
1639             self.verbose=verbose
1640             rtol=self.getTolerance()
1641             atol=self.getAbsoluteTolerance()
1642         if self.adaptSubTolerance(): self.setSubProblemTolerance()
1643             correction_step=0
1644             converged=False
1645             while not converged:
1646                  # calculate velocity for current pressure:
1647                  v=self.getV(p,v)
1648                  Bv=self.Bv(v)
1649                  norm_v=self.norm_v(v)
1650                  norm_Bv=self.norm_Bv(Bv)
1651                  ATOL=norm_v*rtol+atol
1652                  if self.verbose: print "HomogeneousSaddlePointProblem: norm v= %e, norm_Bv= %e, tolerance = %e."%(norm_v, norm_Bv,ATOL)
1653                  if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."
1654                  if norm_Bv <= ATOL:
1655                     converged=True
1656                  else:
1657                     correction_step+=1
1658                     if correction_step>max_correction_steps:
1659                          raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1660                     dp=self.solve_prec(Bv)
1661                     if usePCG:
1662                       norm2=self.inner_pBv(dp,Bv)
1663                       if norm2<0: raise ValueError,"negative PCG norm."
1664                       norm2=math.sqrt(norm2)
1665                     else:
1666                       norm2=self.norm_p(dp)
1667                     ATOL_ITER=ATOL/norm_Bv*norm2*0.5
1668                     if self.verbose: print "HomogeneousSaddlePointProblem: tolerance for solver: %e"%ATOL_ITER
1669                     if usePCG:
1670                           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)
1671                     else:
1672                           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)
1673             if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached."
1674         return v,p
1675    
1676         @param verbose: switches on the printing out some information        #========================================================================
1677         @type verbose: C{bool}        def setTolerance(self,tolerance=1.e-4):
1678         @note: this method may be overwritten by a particular saddle point problem           """
1679         """           Sets the relative tolerance for (v,p).
        if not isinstance(verbose,bool):  
             raise TypeError("verbose needs to be of type bool.")  
        self.__verbose=verbose  
        self.relaxation=1.  
1680    
1681     def trace(self,text):           @param tolerance: tolerance to be used
1682         """           @type tolerance: non-negative C{float}
1683         prints text if verbose has been set           """
1684             if tolerance<0:
1685                 raise ValueError,"tolerance must be positive."
1686             self.__rtol=tolerance
1687    
1688         @param text: a text message        def getTolerance(self):
1689         @type text: C{str}           """
1690         """           Returns the relative tolerance.
        if self.__verbose: print "%s: %s"%(str(self),text)  
1691    
1692     def solve_f(self,u,p,tol=1.e-8):           @return: relative tolerance
1693         """           @rtype: C{float}
1694         solves           """
1695             return self.__rtol
1696    
1697         A_f du = f(u,p)        def setAbsoluteTolerance(self,tolerance=0.):
1698             """
1699             Sets the absolute tolerance.
1700    
1701         with tolerance C{tol} and return du. A_f is Jacobiean of f with respect to u.           @param tolerance: tolerance to be used
1702             @type tolerance: non-negative C{float}
1703             """
1704             if tolerance<0:
1705                 raise ValueError,"tolerance must be non-negative."
1706             self.__atol=tolerance
1707    
1708         @param u: current approximation of u        def getAbsoluteTolerance(self):
1709         @type u: L{escript.Data}           """
1710         @param p: current approximation of p           Returns the absolute tolerance.
        @type p: L{escript.Data}  
        @param tol: tolerance expected for du  
        @type tol: C{float}  
        @return: increment du  
        @rtype: L{escript.Data}  
        @note: this method has to be overwritten by a particular saddle point problem  
        """  
        pass  
1711    
1712     def solve_g(self,u,tol=1.e-8):           @return: absolute tolerance
1713         """           @rtype: C{float}
1714         solves           """
1715             return self.__atol
1716    
1717         Q_g dp = g(u)        def getSubProblemTolerance(self):
1718             """
1719             Sets the relative tolerance to solve the subproblem(s).
1720    
1721         with Q_g is a preconditioner for A_g A_f^{-1} A_g with  A_g is the jacobiean of g with respect to p.           @param rtol: relative tolerence
1722             @type rtol: positive C{float}
1723             """
1724             return max(200.*util.EPSILON,self.getTolerance()**2)
1725    
1726         @param u: current approximation of u  def MaskFromBoundaryTag(domain,*tags):
1727         @type u: L{escript.Data}     """
1728         @param tol: tolerance expected for dp     Creates a mask on the Solution(domain) function space where the value is
1729         @type tol: C{float}     one for samples that touch the boundary tagged by tags.
        @return: increment dp  
        @rtype: L{escript.Data}  
        @note: this method has to be overwritten by a particular saddle point problem  
        """  
        pass  
1730    
1731     def inner(self,p0,p1):     Usage: m=MaskFromBoundaryTag(domain, "left", "right")
        """  
        inner product of p0 and p1 approximating p. Typically this returns integrate(p0*p1)  
        @return: inner product of p0 and p1  
        @rtype: C{float}  
        """  
        pass  
1732    
1733     subiter_max=3     @param domain: domain to be used
1734     def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):     @type domain: L{escript.Domain}
1735          """     @param tags: boundary tags
1736          runs the solver     @type tags: C{str}
1737       @return: a mask which marks samples that are touching the boundary tagged
1738                by any of the given tags
1739       @rtype: L{escript.Data} of rank 0
1740       """
1741       pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1742       d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
1743       for t in tags: d.setTaggedValue(t,1.)
1744       pde.setValue(y=d)
1745       return util.whereNonZero(pde.getRightHandSide())
1746    
1747          @param u0: initial guess for C{u}  def MaskFromTag(domain,*tags):
         @type u0: L{esys.escript.Data}  
         @param p0: initial guess for C{p}  
         @type p0: L{esys.escript.Data}  
         @param tolerance: tolerance for relative error in C{u} and C{p}  
         @type tolerance: positive C{float}  
         @param tolerance_u: tolerance for relative error in C{u} if different from C{tolerance}  
         @type tolerance_u: positive C{float}  
         @param iter_max: maximum number of iteration steps.  
         @type iter_max: C{int}  
         @param accepted_reduction: if the norm  g cannot be reduced by C{accepted_reduction} backtracking to adapt the  
                                    relaxation factor. If C{accepted_reduction=None} no backtracking is used.  
         @type accepted_reduction: positive C{float} or C{None}  
         @param relaxation: initial relaxation factor. If C{relaxation==None}, the last relaxation factor is used.  
         @type relaxation: C{float} or C{None}  
         """  
         tol=1.e-2  
         if tolerance_u==None: tolerance_u=tolerance  
         if not relaxation==None: self.relaxation=relaxation  
         if accepted_reduction ==None:  
               angle_limit=0.  
         elif accepted_reduction>=1.:  
               angle_limit=0.  
         else:  
               angle_limit=util.sqrt(1-accepted_reduction**2)  
         self.iter=0  
         u=u0  
         p=p0  
         #  
         #   initialize things:  
         #  
         converged=False  
         #  
         #  start loop:  
         #  
         #  initial search direction is g  
         #  
         while not converged :  
             if self.iter>iter_max:  
                 raise ArithmeticError("no convergence after %s steps."%self.iter)  
             f_new=self.solve_f(u,p,tol)  
             norm_f_new = util.Lsup(f_new)  
             u_new=u-f_new  
             g_new=self.solve_g(u_new,tol)  
             self.iter+=1  
             norm_g_new = util.sqrt(self.inner(g_new,g_new))  
             if norm_f_new==0. and norm_g_new==0.: return u, p  
             if self.iter>1 and not accepted_reduction==None:  
                #  
                #   did we manage to reduce the norm of G? I  
                #   if not we start a backtracking procedure  
                #  
                # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g  
                if norm_g_new > accepted_reduction * norm_g:  
                   sub_iter=0  
                   s=self.relaxation  
                   d=g  
                   g_last=g  
                   self.trace("    start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))  
                   while sub_iter < self.subiter_max and  norm_g_new > accepted_reduction * norm_g:  
                      dg= g_new-g_last  
                      norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)  
                      rad=self.inner(g_new,dg)/self.relaxation  
                      # print "   ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit  
                      # print "   ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g  
                      if abs(rad) < norm_dg*norm_g_new * angle_limit:  
                          if sub_iter>0: self.trace("    no further improvements expected from backtracking.")  
                          break  
                      r=self.relaxation  
                      self.relaxation= - rad/norm_dg**2  
                      s+=self.relaxation  
                      #####  
                      # a=g_new+self.relaxation*dg/r  
                      # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation  
                      #####  
                      g_last=g_new  
                      p+=self.relaxation*d  
                      f_new=self.solve_f(u,p,tol)  
                      u_new=u-f_new  
                      g_new=self.solve_g(u_new,tol)  
                      self.iter+=1  
                      norm_f_new = util.Lsup(f_new)  
                      norm_g_new = util.sqrt(self.inner(g_new,g_new))  
                      # print "   ",sub_iter," new g norm",norm_g_new  
                      self.trace("    %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))  
                      #  
                      #   can we expect reduction of g?  
                      #  
                      # u_last=u_new  
                      sub_iter+=1  
                   self.relaxation=s  
             #  
             #  check for convergence:  
             #  
             norm_u_new = util.Lsup(u_new)  
             p_new=p+self.relaxation*g_new  
             norm_p_new = util.sqrt(self.inner(p_new,p_new))  
             self.trace("%s th step: f = %s, f/u = %s, g = %s, g/p = %s, relaxation = %s."%(self.iter, norm_f_new ,norm_f_new/norm_u_new, norm_g_new, norm_g_new/norm_p_new, self.relaxation))  
   
             if self.iter>1:  
                dg2=g_new-g  
                df2=f_new-f  
                norm_dg2=util.sqrt(self.inner(dg2,dg2))  
                norm_df2=util.Lsup(df2)  
                # print norm_g_new, norm_g, norm_dg, norm_p, tolerance  
                tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new  
                tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new  
                if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:  
                    converged=True  
             f, norm_f, u, norm_u, g, norm_g, p, norm_p = f_new, norm_f_new, u_new, norm_u_new, g_new, norm_g_new, p_new, norm_p_new  
         self.trace("convergence after %s steps."%self.iter)  
         return u,p  
 #   def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):  
 #      tol=1.e-2  
 #      iter=0  
 #      converged=False  
 #      u=u0*1.  
 #      p=p0*1.  
 #      while not converged and iter<iter_max:  
 #          du=self.solve_f(u,p,tol)  
 #          u-=du  
 #          norm_du=util.Lsup(du)  
 #          norm_u=util.Lsup(u)  
 #          
 #          dp=self.relaxation*self.solve_g(u,tol)  
 #          p+=dp  
 #          norm_dp=util.sqrt(self.inner(dp,dp))  
 #          norm_p=util.sqrt(self.inner(p,p))  
 #          print iter,"-th step rel. errror u,p= (%s,%s) (%s,%s)(%s,%s)"%(norm_du,norm_dp,norm_du/norm_u,norm_dp/norm_p,norm_u,norm_p)  
 #          iter+=1  
 #  
 #          converged = (norm_du <= tolerance*norm_u) and  (norm_dp <= tolerance*norm_p)  
 #      if converged:  
 #          print "convergence after %s steps."%iter  
 #      else:  
 #          raise ArithmeticError("no convergence after %s steps."%iter)  
 #  
 #      return u,p  
             
 def MaskFromBoundaryTag(function_space,*tags):  
1748     """     """
1749     create a mask on the given function space which one for samples     Creates a mask on the Solution(domain) function space where the value is
1750     that touch the boundary tagged by tags.     one for samples that touch regions tagged by tags.
1751    
1752     usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")     Usage: m=MaskFromTag(domain, "ham")
1753    
1754     @param function_space: a given function space     @param domain: domain to be used
1755     @type function_space: L{escript.FunctionSpace}     @type domain: L{escript.Domain}
1756     @param tags: boundray tags     @param tags: boundary tags
1757     @type tags: C{str}     @type tags: C{str}
1758     @return: a mask which marks samples used by C{function_space} that are touching the     @return: a mask which marks samples that are touching the boundary tagged
1759              boundary tagged by any of the given tags.              by any of the given tags
1760     @rtype: L{escript.Data} of rank 0     @rtype: L{escript.Data} of rank 0
1761     """     """
1762     pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1763     d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))     d=escript.Scalar(0.,escript.Function(domain))
1764     for t in tags: d.setTaggedValue(t,1.)     for t in tags: d.setTaggedValue(t,1.)
1765     pde.setValue(y=d)     pde.setValue(Y=d)
1766     out=util.whereNonZero(pde.getRightHandSide())     return util.whereNonZero(pde.getRightHandSide())
    if out.getFunctionSpace() == function_space:  
       return out  
    else:  
       return util.whereNonZero(util.interpolate(out,function_space))  
   
1767    
1768    

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

  ViewVC Help
Powered by ViewVC 1.1.26