/[escript]/trunk/escript/py_src/pdetools.py
ViewVC logotype

Diff of /trunk/escript/py_src/pdetools.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

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

  ViewVC Help
Powered by ViewVC 1.1.26