/[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 2245 by gross, Wed Feb 4 06:27:59 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__="http://www.uq.edu.au/esscc/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 numarray
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{numarray.NumArray}, L{Data}
227             @param Y: right hand side
228             @type Y: C{float}, C{int}, C{numarray.NumArray}, L{Data}
229             @param q: location of constraints
230             @type q: C{float}, C{int}, C{numarray.NumArray}, L{Data}
231             @param r: value of solution at locations of constraints
232             @type r: C{float}, C{int}, C{numarray.NumArray}, 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{numarray.NumArray}, L{Data}
262             @param Y: right hand side
263             @type Y: C{float}, C{int}, C{numarray.NumArray}, L{Data}
264             @param q: location of constraints
265             @type q: C{float}, C{int}, C{numarray.NumArray}, L{Data}
266             @param r: value of solution at locations of constraints
267             @type r: C{float}, C{int}, C{numarray.NumArray}, 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=numarray.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{numarray.NumArray} or C{list} of C{numarray.NumArray}
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=data.getValueOfGlobalDataPoint(*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=data.getValueOfGlobalDataPoint(*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           r += (-alpha) * q
526    
527           rhat=Msolve(r)
528           rhat_dot_r_new = bilinearform(rhat, r)
529           beta = rhat_dot_r_new / rhat_dot_r
530           rhat+=beta * d
531           d=rhat
532    
533           rhat_dot_r = rhat_dot_r_new
534           if rhat_dot_r<0: raise NegativeNorm,"negative norm."
535           if verbose: print "PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))
536       if verbose: print "PCG: tolerance reached after %s steps."%iter
537       return x,r
538    
539    class Defect(object):
540        """
541        Defines a non-linear defect F(x) of a variable x.
542        """
543        def __init__(self):
544            """
545            Initializes defect.
546            """
547            self.setDerivativeIncrementLength()
548    
549        def bilinearform(self, x0, x1):
550            """
551            Returns the inner product of x0 and x1
552    
553            @param x0: value for x0
554            @param x1: value for x1
555            @return: the inner product of x0 and x1
556            @rtype: C{float}
557            """
558            return 0
559    
560        def norm(self,x):
561            """
562            Returns the norm of argument C{x}.
563    
564            @param x: a value
565            @return: norm of argument x
566            @rtype: C{float}
567            @note: by default C{sqrt(self.bilinearform(x,x)} is returned.
568            """
569            s=self.bilinearform(x,x)
570            if s<0: raise NegativeNorm,"negative norm."
571            return math.sqrt(s)
572    
573        def eval(self,x):
574            """
575            Returns the value F of a given C{x}.
576    
577            @param x: value for which the defect C{F} is evaluated
578            @return: value of the defect at C{x}
579            """
580            return 0
581    
582        def __call__(self,x):
583            return self.eval(x)
584    
585        def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):
586            """
587            Sets the relative length of the increment used to approximate the
588            derivative of the defect. The increment is inc*norm(x)/norm(v)*v in the
589            direction of v with x as a starting point.
590    
591            @param inc: relative increment length
592            @type inc: positive C{float}
593            """
594            if inc<=0: raise ValueError,"positive increment required."
595            self.__inc=inc
596    
597        def getDerivativeIncrementLength(self):
598            """
599            Returns the relative increment length used to approximate the
600            derivative of the defect.
601            @return: value of the defect at C{x}
602            @rtype: positive C{float}
603            """
604            return self.__inc
605    
606        def derivative(self, F0, x0, v, v_is_normalised=True):
607            """
608            Returns the directional derivative at C{x0} in the direction of C{v}.
609    
610            @param F0: value of this defect at x0
611            @param x0: value at which derivative is calculated
612            @param v: direction
613            @param v_is_normalised: True to indicate that C{v} is nomalized
614                                    (self.norm(v)=0)
615            @return: derivative of this defect at x0 in the direction of C{v}
616            @note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is
617                   used but this method maybe overwritten to use exact evaluation.
618            """
619            normx=self.norm(x0)
620            if normx>0:
621                 epsnew = self.getDerivativeIncrementLength() * normx
622            else:
623                 epsnew = self.getDerivativeIncrementLength()
624            if not v_is_normalised:
625                normv=self.norm(v)
626                if normv<=0:
627                   return F0*0
628                else:
629                   epsnew /= normv
630            F1=self.eval(x0 + epsnew * v)
631            return (F1-F0)/epsnew
632    
633    ######################################
634    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):
635       """
636       Solves a non-linear problem M{F(x)=0} for unknown M{x} using the stopping
637       criterion:
638    
639       M{norm(F(x) <= atol + rtol * norm(F(x0)}
640    
641       where M{x0} is the initial guess.
642    
643       @param defect: object defining the function M{F}. C{defect.norm} defines the
644                      M{norm} used in the stopping criterion.
645       @type defect: L{Defect}
646       @param x: initial guess for the solution, C{x} is altered.
647       @type x: any object type allowing basic operations such as
648                C{numarray.NumArray}, L{Data}
649       @param iter_max: maximum number of iteration steps
650       @type iter_max: positive C{int}
651       @param sub_iter_max: maximum number of inner iteration steps
652       @type sub_iter_max: positive C{int}
653       @param atol: absolute tolerance for the solution
654       @type atol: positive C{float}
655       @param rtol: relative tolerance for the solution
656       @type rtol: positive C{float}
657       @param gamma: tolerance safety factor for inner iteration
658       @type gamma: positive C{float}, less than 1
659       @param sub_tol_max: upper bound for inner tolerance
660       @type sub_tol_max: positive C{float}, less than 1
661       @return: an approximation of the solution with the desired accuracy
662       @rtype: same type as the initial guess
663       """
664       lmaxit=iter_max
665       if atol<0: raise ValueError,"atol needs to be non-negative."
666       if rtol<0: raise ValueError,"rtol needs to be non-negative."
667       if rtol+atol<=0: raise ValueError,"rtol or atol needs to be non-negative."
668       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
669       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
670    
671       F=defect(x)
672       fnrm=defect.norm(F)
673       stop_tol=atol + rtol*fnrm
674       sub_tol=sub_tol_max
675       if verbose: print "NewtonGMRES: initial residual = %e."%fnrm
676       if verbose: print "             tolerance = %e."%sub_tol
677       iter=1
678       #
679       # main iteration loop
680       #
681       while not fnrm<=stop_tol:
682                if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
683                #
684            #   adjust sub_tol_
685            #
686                if iter > 1:
687               rat=fnrm/fnrmo
688                   sub_tol_old=sub_tol
689               sub_tol=gamma*rat**2
690               if gamma*sub_tol_old**2 > .1: sub_tol=max(sub_tol,gamma*sub_tol_old**2)
691               sub_tol=max(min(sub_tol,sub_tol_max), .5*stop_tol/fnrm)
692            #
693            # calculate newton increment xc
694                #     if iter_max in __FDGMRES is reached MaxIterReached is thrown
695                #     if iter_restart -1 is returned as sub_iter
696                #     if  atol is reached sub_iter returns the numer of steps performed to get there
697                #
698                #
699                if verbose: print "             subiteration (GMRES) is called with relative tolerance %e."%sub_tol
700                try:
701                   xc, sub_iter=__FDGMRES(F, defect, x, sub_tol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)
702                except MaxIterReached:
703                   raise MaxIterReached,"maximum number of %s steps reached."%iter_max
704                if sub_iter<0:
705                   iter+=sub_iter_max
706                else:
707                   iter+=sub_iter
708                # ====
709            x+=xc
710                F=defect(x)
711            iter+=1
712                fnrmo, fnrm=fnrm, defect.norm(F)
713                if verbose: print "             step %s: residual %e."%(iter,fnrm)
714       if verbose: print "NewtonGMRES: completed after %s steps."%iter
715       return x
716    
717    def __givapp(c,s,vin):
718        """
719        Applies a sequence of Givens rotations (c,s) recursively to the vector
720        C{vin}
721    
722        @warning: C{vin} is altered.
723        """
724        vrot=vin
725        if isinstance(c,float):
726            vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
727        else:
728            for i in range(len(c)):
729                w1=c[i]*vrot[i]-s[i]*vrot[i+1]
730            w2=s[i]*vrot[i]+c[i]*vrot[i+1]
731                vrot[i:i+2]=w1,w2
732        return vrot
733    
734    def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
735       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
736       c=numarray.zeros(iter_restart,numarray.Float64)
737       s=numarray.zeros(iter_restart,numarray.Float64)
738       g=numarray.zeros(iter_restart,numarray.Float64)
739       v=[]
740    
741       rho=defect.norm(F0)
742       if rho<=0.: return x0*0
743    
744       v.append(-F0/rho)
745       g[0]=rho
746       iter=0
747       while rho > atol and iter<iter_restart-1:
748            if iter  >= iter_max:
749                raise MaxIterReached,"maximum number of %s steps reached."%iter_max
750    
751            p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
752            v.append(p)
753    
754            v_norm1=defect.norm(v[iter+1])
755    
756            # Modified Gram-Schmidt
757            for j in range(iter+1):
758                h[j,iter]=defect.bilinearform(v[j],v[iter+1])
759                v[iter+1]-=h[j,iter]*v[j]
760    
761            h[iter+1,iter]=defect.norm(v[iter+1])
762            v_norm2=h[iter+1,iter]
763    
764            # Reorthogonalize if needed
765            if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
766                for j in range(iter+1):
767                    hr=defect.bilinearform(v[j],v[iter+1])
768                    h[j,iter]=h[j,iter]+hr
769                    v[iter+1] -= hr*v[j]
770    
771                v_norm2=defect.norm(v[iter+1])
772                h[iter+1,iter]=v_norm2
773            #   watch out for happy breakdown
774            if not v_norm2 == 0:
775                v[iter+1]=v[iter+1]/h[iter+1,iter]
776    
777            #   Form and store the information for the new Givens rotation
778            if iter > 0 :
779                hhat=numarray.zeros(iter+1,numarray.Float64)
780                for i in range(iter+1) : hhat[i]=h[i,iter]
781                hhat=__givapp(c[0:iter],s[0:iter],hhat);
782                for i in range(iter+1) : h[i,iter]=hhat[i]
783    
784            mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
785    
786            if mu!=0 :
787                c[iter]=h[iter,iter]/mu
788                s[iter]=-h[iter+1,iter]/mu
789                h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
790                h[iter+1,iter]=0.0
791                g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])
792    
793            # Update the residual norm
794            rho=abs(g[iter+1])
795            iter+=1
796    
797       # At this point either iter > iter_max or rho < tol.
798       # It's time to compute x and leave.
799       if iter > 0 :
800         y=numarray.zeros(iter,numarray.Float64)
801         y[iter-1] = g[iter-1] / h[iter-1,iter-1]
802         if iter > 1 :
803            i=iter-2
804            while i>=0 :
805              y[i] = ( g[i] - numarray.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
806              i=i-1
807         xhat=v[iter-1]*y[iter-1]
808         for i in range(iter-1):
809        xhat += v[i]*y[i]
810       else :
811          xhat=v[0] * 0
812    
813       if iter<iter_restart-1:
814          stopped=iter
815       else:
816          stopped=-1
817    
818       return xhat,stopped
819    
820    def GMRES(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100, iter_restart=20, verbose=False):
821       """
822       Solver for
823    
824       M{Ax=b}
825    
826       with a general operator A (more details required!).
827       It uses the generalized minimum residual method (GMRES).
828    
829       The iteration is terminated if
830    
831       M{|r| <= atol+rtol*|r0|}
832    
833       where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
834    
835       M{|r| = sqrt( bilinearform(r,r))}
836    
837       @param r: initial residual M{r=b-Ax}. C{r} is altered.
838       @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
839       @param x: an initial guess for the solution
840       @type x: same like C{r}
841       @param Aprod: returns the value Ax
842       @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
843                    argument C{x}. The returned object needs to be of the same
844                    type like argument C{r}.
845       @param bilinearform: inner product C{<x,r>}
846       @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
847                           type like argument C{x} and C{r}. The returned value is
848                           a C{float}.
849       @param atol: absolute tolerance
850       @type atol: non-negative C{float}
851       @param rtol: relative tolerance
852       @type rtol: non-negative C{float}
853       @param iter_max: maximum number of iteration steps
854       @type iter_max: C{int}
855       @param iter_restart: in order to save memory the orthogonalization process
856                            is terminated after C{iter_restart} steps and the
857                            iteration is restarted.
858       @type iter_restart: C{int}
859       @return: the solution approximation and the corresponding residual
860       @rtype: C{tuple}
861       @warning: C{r} and C{x} are altered.
862       """
863       m=iter_restart
864       restarted=False
865       iter=0
866       if rtol>0:
867          r_dot_r = bilinearform(r, r)
868          if r_dot_r<0: raise NegativeNorm,"negative norm."
869          atol2=atol+rtol*math.sqrt(r_dot_r)
870          if verbose: print "GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)
871       else:
872          atol2=atol
873          if verbose: print "GMRES: absolute tolerance = %e"%atol2
874       if atol2<=0:
875          raise ValueError,"Non-positive tolarance."
876    
877       while True:
878          if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
879          if restarted:
880             r2 = r-Aprod(x-x2)
881          else:
882             r2=1*r
883          x2=x*1.
884          x,stopped=_GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose)
885          iter+=iter_restart
886          if stopped: break
887          if verbose: print "GMRES: restart."
888          restarted=True
889       if verbose: print "GMRES: tolerance has reached."
890       return x
891    
892    def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False):
893       iter=0
894    
895       h=numarray.zeros((iter_restart+1,iter_restart),numarray.Float64)
896       c=numarray.zeros(iter_restart,numarray.Float64)
897       s=numarray.zeros(iter_restart,numarray.Float64)
898       g=numarray.zeros(iter_restart+1,numarray.Float64)
899       v=[]
900    
901       r_dot_r = bilinearform(r, r)
902       if r_dot_r<0: raise NegativeNorm,"negative norm."
903       rho=math.sqrt(r_dot_r)
904    
905       v.append(r/rho)
906       g[0]=rho
907    
908       if verbose: print "GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)
909       while not (rho<=atol or iter==iter_restart):
910    
911        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
912    
913        p=Aprod(v[iter])
914        v.append(p)
915    
916        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
917    
918    # Modified Gram-Schmidt
919        for j in range(iter+1):
920          h[j,iter]=bilinearform(v[j],v[iter+1])
921          v[iter+1]-=h[j,iter]*v[j]
922    
923        h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
924        v_norm2=h[iter+1,iter]
925    
926    # Reorthogonalize if needed
927        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
928         for j in range(iter+1):
929            hr=bilinearform(v[j],v[iter+1])
930                h[j,iter]=h[j,iter]+hr
931                v[iter+1] -= hr*v[j]
932    
933         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
934         h[iter+1,iter]=v_norm2
935    
936    #   watch out for happy breakdown
937            if not v_norm2 == 0:
938             v[iter+1]=v[iter+1]/h[iter+1,iter]
939    
940    #   Form and store the information for the new Givens rotation
941        if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])
942        mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
943    
944        if mu!=0 :
945            c[iter]=h[iter,iter]/mu
946            s[iter]=-h[iter+1,iter]/mu
947            h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
948            h[iter+1,iter]=0.0
949            g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])
950    # Update the residual norm
951    
952            rho=abs(g[iter+1])
953            if verbose: print "GMRES: iteration step %s: residual %e"%(iter,rho)
954        iter+=1
955    
956    # At this point either iter > iter_max or rho < tol.
957    # It's time to compute x and leave.
958    
959       if verbose: print "GMRES: iteration stopped after %s step."%iter
960       if iter > 0 :
961         y=numarray.zeros(iter,numarray.Float64)
962         y[iter-1] = g[iter-1] / h[iter-1,iter-1]
963         if iter > 1 :
964            i=iter-2
965            while i>=0 :
966              y[i] = ( g[i] - numarray.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
967              i=i-1
968         xhat=v[iter-1]*y[iter-1]
969         for i in range(iter-1):
970        xhat += v[i]*y[i]
971       else:
972         xhat=v[0] * 0
973    
974       x += xhat
975       if iter<iter_restart-1:
976          stopped=True
977       else:
978          stopped=False
979    
980       return x,stopped
981    
982    def MINRES(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
983        """
984        Solver for
985    
986        M{Ax=b}
987    
988        with a symmetric and positive definite operator A (more details required!).
989        It uses the minimum residual method (MINRES) with preconditioner M
990        providing an approximation of A.
991    
992        The iteration is terminated if
993    
994        M{|r| <= atol+rtol*|r0|}
995    
996        where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
997    
998        M{|r| = sqrt( bilinearform(Msolve(r),r))}
999    
1000        For details on the preconditioned conjugate gradient method see the book:
1001    
1002        I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
1003        T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
1004        C. Romine, and H. van der Vorst}.
1005    
1006        @param r: initial residual M{r=b-Ax}. C{r} is altered.
1007        @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1008        @param x: an initial guess for the solution
1009        @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1010        @param Aprod: returns the value Ax
1011        @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
1012                     argument C{x}. The returned object needs to be of the same
1013                     type like argument C{r}.
1014        @param Msolve: solves Mx=r
1015        @type Msolve: function C{Msolve(r)} where C{r} is of the same type like
1016                      argument C{r}. The returned object needs to be of the same
1017                      type like argument C{x}.
1018        @param bilinearform: inner product C{<x,r>}
1019        @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
1020                            type like argument C{x} and C{r} is. The returned value
1021                            is a C{float}.
1022        @param atol: absolute tolerance
1023        @type atol: non-negative C{float}
1024        @param rtol: relative tolerance
1025        @type rtol: non-negative C{float}
1026        @param iter_max: maximum number of iteration steps
1027        @type iter_max: C{int}
1028        @return: the solution approximation and the corresponding residual
1029        @rtype: C{tuple}
1030        @warning: C{r} and C{x} are altered.
1031        """
1032        #------------------------------------------------------------------
1033        # Set up y and v for the first Lanczos vector v1.
1034        # y  =  beta1 P' v1,  where  P = C**(-1).
1035        # v is really P' v1.
1036        #------------------------------------------------------------------
1037        r1    = r
1038        y = Msolve(r)
1039        beta1 = bilinearform(y,r)
1040    
1041        if beta1< 0: raise NegativeNorm,"negative norm."
1042    
1043        #  If r = 0 exactly, stop with x
1044        if beta1==0: return x
1045    
1046        if beta1> 0: beta1  = math.sqrt(beta1)
1047    
1048        #------------------------------------------------------------------
1049        # Initialize quantities.
1050        # ------------------------------------------------------------------
1051        iter   = 0
1052        Anorm = 0
1053        ynorm = 0
1054        oldb   = 0
1055        beta   = beta1
1056        dbar   = 0
1057        epsln  = 0
1058        phibar = beta1
1059        rhs1   = beta1
1060        rhs2   = 0
1061        rnorm  = phibar
1062        tnorm2 = 0
1063        ynorm2 = 0
1064        cs     = -1
1065        sn     = 0
1066        w      = r*0.
1067        w2     = r*0.
1068        r2     = r1
1069        eps    = 0.0001
1070    
1071        #---------------------------------------------------------------------
1072        # Main iteration loop.
1073        # --------------------------------------------------------------------
1074        while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL
1075    
1076        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1077            iter    = iter  +  1
1078    
1079            #-----------------------------------------------------------------
1080            # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
1081            # The general iteration is similar to the case k = 1 with v0 = 0:
1082            #
1083            #   p1      = Operator * v1  -  beta1 * v0,
1084            #   alpha1  = v1'p1,
1085            #   q2      = p2  -  alpha1 * v1,
1086            #   beta2^2 = q2'q2,
1087            #   v2      = (1/beta2) q2.
1088            #
1089            # Again, y = betak P vk,  where  P = C**(-1).
1090            #-----------------------------------------------------------------
1091            s = 1/beta                 # Normalize previous vector (in y).
1092            v = s*y                    # v = vk if P = I
1093    
1094            y      = Aprod(v)
1095    
1096            if iter >= 2:
1097              y = y - (beta/oldb)*r1
1098    
1099            alfa   = bilinearform(v,y)              # alphak
1100            y      += (- alfa/beta)*r2
1101            r1     = r2
1102            r2     = y
1103            y = Msolve(r2)
1104            oldb   = beta                           # oldb = betak
1105            beta   = bilinearform(y,r2)             # beta = betak+1^2
1106            if beta < 0: raise NegativeNorm,"negative norm."
1107    
1108            beta   = math.sqrt( beta )
1109            tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
1110    
1111            if iter==1:                 # Initialize a few things.
1112              gmax   = abs( alfa )      # alpha1
1113              gmin   = gmax             # alpha1
1114    
1115            # Apply previous rotation Qk-1 to get
1116            #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
1117            #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
1118    
1119            oldeps = epsln
1120            delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
1121            gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
1122            epsln  =               sn * beta  # epsln2 = 0         epslnk+1
1123            dbar   =            -  cs * beta  # dbar 2 = beta2     dbar k+1
1124    
1125            # Compute the next plane rotation Qk
1126    
1127            gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
1128            gamma  = max(gamma,eps)
1129            cs     = gbar / gamma             # ck
1130            sn     = beta / gamma             # sk
1131            phi    = cs * phibar              # phik
1132            phibar = sn * phibar              # phibark+1
1133    
1134            # Update  x.
1135    
1136            denom = 1/gamma
1137            w1    = w2
1138            w2    = w
1139            w     = (v - oldeps*w1 - delta*w2) * denom
1140            x     +=  phi*w
1141    
1142            # Go round again.
1143    
1144            gmax   = max(gmax,gamma)
1145            gmin   = min(gmin,gamma)
1146            z      = rhs1 / gamma
1147            ynorm2 = z*z  +  ynorm2
1148            rhs1   = rhs2 -  delta*z
1149            rhs2   =      -  epsln*z
1150    
1151            # Estimate various norms and test for convergence.
1152    
1153            Anorm  = math.sqrt( tnorm2 )
1154            ynorm  = math.sqrt( ynorm2 )
1155    
1156            rnorm  = phibar
1157    
1158        return x
1159    
1160    def TFQMR(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1161      """
1162      Solver for
1163    
1164      M{Ax=b}
1165    
1166      with a general operator A (more details required!).
1167      It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).
1168    
1169      The iteration is terminated if
1170    
1171      M{|r| <= atol+rtol*|r0|}
1172    
1173      where M{r0} is the initial residual and M{|.|} is the energy norm. In fact
1174    
1175      M{|r| = sqrt( bilinearform(r,r))}
1176    
1177      @param r: initial residual M{r=b-Ax}. C{r} is altered.
1178      @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1179      @param x: an initial guess for the solution
1180      @type x: same like C{r}
1181      @param Aprod: returns the value Ax
1182      @type Aprod: function C{Aprod(x)} where C{x} is of the same object like
1183                   argument C{x}. The returned object needs to be of the same type
1184                   like argument C{r}.
1185      @param bilinearform: inner product C{<x,r>}
1186      @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same
1187                          type like argument C{x} and C{r}. The returned value is
1188                          a C{float}.
1189      @param atol: absolute tolerance
1190      @type atol: non-negative C{float}
1191      @param rtol: relative tolerance
1192      @type rtol: non-negative C{float}
1193      @param iter_max: maximum number of iteration steps
1194      @type iter_max: C{int}
1195      @rtype: C{tuple}
1196      @warning: C{r} and C{x} are altered.
1197      """
1198      u1=0
1199      u2=0
1200      y1=0
1201      y2=0
1202    
1203      w = r
1204      y1 = r
1205      iter = 0
1206      d = 0
1207      v = Aprod(y1)
1208      u1 = v
1209    
1210      theta = 0.0;
1211      eta = 0.0;
1212      rho=bilinearform(r,r)
1213      if rho < 0: raise NegativeNorm,"negative norm."
1214      tau = math.sqrt(rho)
1215      norm_r0=tau
1216      while tau>atol+rtol*norm_r0:
1217        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1218    
1219        sigma = bilinearform(r,v)
1220        if sigma == 0.0: raise IterationBreakDown,'TFQMR breakdown, sigma=0'
1221    
1222        alpha = rho / sigma
1223    
1224        for j in range(2):
1225    #
1226    #   Compute y2 and u2 only if you have to
1227    #
1228          if ( j == 1 ):
1229            y2 = y1 - alpha * v
1230            u2 = Aprod(y2)
1231    
1232          m = 2 * (iter+1) - 2 + (j+1)
1233          if j==0:
1234             w = w - alpha * u1
1235             d = y1 + ( theta * theta * eta / alpha ) * d
1236          if j==1:
1237             w = w - alpha * u2
1238             d = y2 + ( theta * theta * eta / alpha ) * d
1239    
1240          theta = math.sqrt(bilinearform(w,w))/ tau
1241          c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1242          tau = tau * theta * c
1243          eta = c * c * alpha
1244          x = x + eta * d
1245    #
1246    #  Try to terminate the iteration at each pass through the loop
1247    #
1248        if rho == 0.0: raise IterationBreakDown,'TFQMR breakdown, rho=0'
1249    
1250        rhon = bilinearform(r,w)
1251        beta = rhon / rho;
1252        rho = rhon;
1253        y1 = w + beta * y2;
1254        u1 = Aprod(y1)
1255        v = u1 + beta * ( u2 + beta * v )
1256    
1257        iter += 1
1258    
1259      return x
1260    
1261    
1262    #############################################
1263    
1264    class ArithmeticTuple(object):
1265       """
1266       Tuple supporting inplace update x+=y and scaling x=a*y where C{x,y} is an
1267       ArithmeticTuple and C{a} is a float.
1268    
1269       Example of usage::
1270    
1271           from esys.escript import Data
1272           from numarray import array
1273           a=Data(...)
1274           b=array([1.,4.])
1275           x=ArithmeticTuple(a,b)
1276           y=5.*x
1277    
1278       """
1279       def __init__(self,*args):
1280           """
1281           Initializes object with elements C{args}.
1282    
1283           @param args: tuple of objects that support inplace add (x+=y) and
1284                        scaling (x=a*y)
1285           """
1286           self.__items=list(args)
1287    
1288       def __len__(self):
1289           """
1290           Returns the number of items.
1291    
1292           @return: number of items
1293           @rtype: C{int}
1294           """
1295           return len(self.__items)
1296    
1297       def __getitem__(self,index):
1298           """
1299           Returns item at specified position.
1300    
1301           @param index: index of item to be returned
1302           @type index: C{int}
1303           @return: item with index C{index}
1304           """
1305           return self.__items.__getitem__(index)
1306    
1307       def __mul__(self,other):
1308           """
1309           Scales by C{other} from the right.
1310    
1311           @param other: scaling factor
1312           @type other: C{float}
1313           @return: itemwise self*other
1314           @rtype: L{ArithmeticTuple}
1315           """
1316           out=[]
1317           try:
1318               l=len(other)
1319               if l!=len(self):
1320                   raise ValueError,"length of arguments don't match."
1321               for i in range(l): out.append(self[i]*other[i])
1322           except TypeError:
1323               for i in range(len(self)): out.append(self[i]*other)
1324           return ArithmeticTuple(*tuple(out))
1325    
1326       def __rmul__(self,other):
1327           """
1328           Scales by C{other} from the left.
1329    
1330           @param other: scaling factor
1331           @type other: C{float}
1332           @return: itemwise other*self
1333           @rtype: L{ArithmeticTuple}
1334           """
1335           out=[]
1336           try:
1337               l=len(other)
1338               if l!=len(self):
1339                   raise ValueError,"length of arguments don't match."
1340               for i in range(l): out.append(other[i]*self[i])
1341           except TypeError:
1342               for i in range(len(self)): out.append(other*self[i])
1343           return ArithmeticTuple(*tuple(out))
1344    
1345       def __div__(self,other):
1346           """
1347           Scales by (1/C{other}) from the right.
1348    
1349           @param other: scaling factor
1350           @type other: C{float}
1351           @return: itemwise self/other
1352           @rtype: L{ArithmeticTuple}
1353           """
1354           return self*(1/other)
1355    
1356       def __rdiv__(self,other):
1357           """
1358           Scales by (1/C{other}) from the left.
1359    
1360           @param other: scaling factor
1361           @type other: C{float}
1362           @return: itemwise other/self
1363           @rtype: L{ArithmeticTuple}
1364           """
1365           out=[]
1366           try:
1367               l=len(other)
1368               if l!=len(self):
1369                   raise ValueError,"length of arguments don't match."
1370               for i in range(l): out.append(other[i]/self[i])
1371           except TypeError:
1372               for i in range(len(self)): out.append(other/self[i])
1373           return ArithmeticTuple(*tuple(out))
1374    
1375       def __iadd__(self,other):
1376           """
1377           Inplace addition of C{other} to self.
1378    
1379           @param other: increment
1380           @type other: C{ArithmeticTuple}
1381           """
1382           if len(self) != len(other):
1383               raise ValueError,"tuple lengths must match."
1384           for i in range(len(self)):
1385               self.__items[i]+=other[i]
1386           return self
1387    
1388       def __add__(self,other):
1389           """
1390           Adds C{other} to self.
1391    
1392           @param other: increment
1393           @type other: C{ArithmeticTuple}
1394           """
1395           out=[]
1396           try:
1397               l=len(other)
1398               if l!=len(self):
1399                   raise ValueError,"length of arguments don't match."
1400               for i in range(l): out.append(self[i]+other[i])
1401           except TypeError:
1402               for i in range(len(self)): out.append(self[i]+other)
1403           return ArithmeticTuple(*tuple(out))
1404    
1405       def __sub__(self,other):
1406           """
1407           Subtracts C{other} from self.
1408    
1409           @param other: decrement
1410           @type other: C{ArithmeticTuple}
1411           """
1412           out=[]
1413           try:
1414               l=len(other)
1415               if l!=len(self):
1416                   raise ValueError,"length of arguments don't match."
1417               for i in range(l): out.append(self[i]-other[i])
1418           except TypeError:
1419               for i in range(len(self)): out.append(self[i]-other)
1420           return ArithmeticTuple(*tuple(out))
1421    
1422       def __isub__(self,other):
1423           """
1424           Inplace subtraction of C{other} from self.
1425    
1426           @param other: decrement
1427           @type other: C{ArithmeticTuple}
1428           """
1429           if len(self) != len(other):
1430               raise ValueError,"tuple length must match."
1431           for i in range(len(self)):
1432               self.__items[i]-=other[i]
1433           return self
1434    
1435       def __neg__(self):
1436           """
1437           Negates values.
1438           """
1439           out=[]
1440           for i in range(len(self)):
1441               out.append(-self[i])
1442           return ArithmeticTuple(*tuple(out))
1443    
1444    
1445    class HomogeneousSaddlePointProblem(object):
1446          """
1447          This class provides a framework for solving linear homogeneous saddle
1448          point problems of the form::
1449    
1450              M{Av+B^*p=f}
1451              M{Bv     =0}
1452    
1453          for the unknowns M{v} and M{p} and given operators M{A} and M{B} and
1454          given right hand side M{f}. M{B^*} is the adjoint operator of M{B}.
1455          """
1456          def __init__(self,**kwargs):
1457            self.setTolerance()
1458            self.setAbsoluteTolerance()
1459            self.setSubProblemTolerance()
1460    
1461          #=============================================================
1462          def initialize(self):
1463            """
1464            Initializes the problem (overwrite).
1465            """
1466            pass
1467    
1468          def B(self,v):
1469            """
1470            Returns Bv (overwrite).
1471    
1472            @rtype: equal to the type of p
1473            @note: boundary conditions on p should be zero!
1474            """
1475            raise NotImplementedError, "no operator B implemented."
1476    
1477          def inner_pBv(self,p,Bv):
1478             """
1479             Returns inner product of element p and Bv (overwrite).
1480    
1481             @type p: equal to the type of p
1482             @type Bv: equal to the type of result of operator B
1483             @rtype: equal to the type of p
1484             """
1485             raise NotImplementedError,"no inner product for p implemented."
1486    
1487          def inner_p(self,p0,p1):
1488             """
1489             Returns inner product of element p0 and p1 (overwrite).
1490    
1491             @type p0: equal to the type of p
1492             @type p1: equal to the type of p
1493             @rtype: equal to the type of p
1494             """
1495             raise NotImplementedError,"no inner product for p implemented."
1496    
1497          def inner_v(self,v0,v1):
1498             """
1499             Returns inner product of two elements v0 and v1 (overwrite).
1500    
1501             @type v0: equal to the type of v
1502             @type v1: equal to the type of v
1503             @rtype: equal to the type of v
1504             """
1505             raise NotImplementedError,"no inner product for v implemented."
1506             pass
1507    
1508          def solve_A(self,u,p):
1509             """
1510             Solves M{Av=f-Au-B^*p} with accuracy L{self.getSubProblemTolerance()}
1511             (overwrite).
1512    
1513             @rtype: equal to the type of v
1514             @note: boundary conditions on v should be zero!
1515             """
1516             raise NotImplementedError,"no operator A implemented."
1517    
1518          def solve_prec(self,p):
1519             """
1520             Provides a preconditioner for M{BA^{-1}B^*} with accuracy
1521             L{self.getSubProblemTolerance()} (overwrite).
1522    
1523             @rtype: equal to the type of p
1524             @note: boundary conditions on p should be zero!
1525             """
1526             raise NotImplementedError,"no preconditioner for Schur complement implemented."
1527          #=============================================================
1528          def __Aprod_PCG(self,p):
1529              # return (v,Bv) with v=A^-1B*p
1530              #solve Av =B^*p as Av =f-Az-B^*(-p)
1531              v=self.solve_A(self.__z,-p)
1532              return ArithmeticTuple(v, self.B(v))
1533    
1534          def __inner_PCG(self,p,r):
1535             return self.inner_pBv(p,r[1])
1536    
1537          def __Msolve_PCG(self,r):
1538              return self.solve_prec(r[1])
1539          #=============================================================
1540          def __Aprod_GMRES(self,x):
1541              # return w,q  from v, p
1542              # solve Aw =Av+B^*p as Aw =f-A(z-v)-B^*(-p)
1543              #  and  Sq=B(v-w)
1544              v=x[0]
1545              p=x[1]
1546              w=self.solve_A(self.__z-v,-p)
1547              Bw=self.B(w-v)
1548              q=self.solve_prec(Bw)
1549              return ArithmeticTuple(w,q)
1550    
1551          def __inner_GMRES(self,a1,a2):
1552             return self.inner_p(a1[1],a2[1])+self.inner_v(a1[0],a2[0])
1553    
1554          #=============================================================
1555          def norm(self,v,p):
1556            f=self.inner_p(p,p)+self.inner_v(v,v)
1557            if f<0:
1558                raise ValueError,"negative norm."
1559            return math.sqrt(f)
1560    
1561          def solve(self,v,p,max_iter=20, verbose=False, show_details=False, useUzawa=True, iter_restart=20,max_correction_steps=3):
1562             """
1563             Solves the saddle point problem using initial guesses v and p.
1564    
1565             @param v: initial guess for velocity
1566             @param p: initial guess for pressure
1567             @type v: L{Data}
1568             @type p: L{Data}
1569             @param useUzawa: indicates the usage of the Uzawa scheme. Otherwise
1570                              the problem is solved in coupled form.
1571             @param max_iter: maximum number of iteration steps per correction
1572                              attempt
1573             @param verbose: if True, shows information on the progress of the
1574                             saddlepoint problem solver.
1575             @param show_details: if True, shows details of the sub problem solver
1576             @param iter_restart: restart the iteration after C{iter_restart} steps
1577                                  (only used if useUzaw=False)
1578             @param max_correction_steps: maximum number of iteration steps in the
1579                                          attempt to get |Bv| to zero.
1580             @return: new approximations for velocity and pressure
1581             @type useUzawa: C{bool}
1582             @type max_iter: C{int}
1583             @type verbose: C{bool}
1584             @type show_details: C{bool}
1585             @type iter_restart: C{int}
1586             @type max_correction_steps: C{int}
1587             @rtype: C{tuple} of L{Data} objects
1588             """
1589             self.verbose=verbose
1590             self.show_details=show_details and self.verbose
1591             #=====================================================================
1592             # Az=f is solved as A(z-v)=f-Av-B^*0 (typically with z-v=0 on boundary)
1593             self.__z=v+self.solve_A(v,p*0)
1594             # tolerances:
1595             rtol=self.getTolerance()
1596             atol=self.getAbsoluteTolerance()
1597             if useUzawa:
1598                p0=self.solve_prec(self.B(self.__z))
1599                f0=self.norm(self.__z,p0)
1600             else:
1601                f0=util.sqrt(self.inner_v(self.__z,self.__z))
1602             if not f0>0: return self.__z, p*0
1603             ATOL=rtol*f0+atol
1604             if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."
1605             if self.verbose: print "saddle point solver: initial residual %e, tolerance = %e."%(f0,ATOL)
1606             # initialization
1607             self.iter=0
1608             correction_step=0
1609             converged=False
1610             # initial guess:
1611             q=p*1
1612             u=v*1
1613             while not converged :
1614                if useUzawa:
1615                   # assume p is known: then v=z-A^-1B^*p
1616                   #
1617                   # which leads to BA^-1B^*p = Bz
1618                   #
1619                   # note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1620                   # we use the (v,Bv) to represent the residual
1621                   #
1622                   # the norm of the right hand side Bv = f0
1623                   #
1624                   #                  and the initial residual
1625                   #
1626                   #     r=Bz-BA^-1B^*q = B(z-A^{-1}B^*q)=Bw with A(w-z)=Az-Az-B^*q = f -A*0 - B^{*}q
1627                   #
1628                   w=self.solve_A(self.__z,q)+self.__z
1629                   if self.verbose: print "enter PCG method (iter_max=%s, atol=%e, subtolerance=%e)"%(max_iter,ATOL, self.getSubProblemTolerance())
1630                   q,r=PCG(ArithmeticTuple(w,self.B(w)),self.__Aprod_PCG,q,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
1631                   u=r[0]
1632                   Bu=r[1]
1633                else:
1634                   #
1635                   #  with v=dv+z
1636                   #
1637                   #   A v + B p = f
1638                   #   B v       = 0
1639                   #
1640                   # apply the preconditioner [[A^{-1} 0][(S^{-1} B A^{-1})  -S^{-1}]]
1641                   #
1642                   w=self.solve_A(u,q)
1643                   if self.verbose: print "enter GMRES (iter_max=%s, atol=%e, subtolerance=%e)"%(max_iter,ATOL,self.getSubProblemTolerance())
1644                   x=GMRES(ArithmeticTuple(w,self.solve_prec(self.B(u+w))),self.__Aprod_GMRES, ArithmeticTuple(u,q), \
1645                             self.__inner_GMRES,atol=ATOL, rtol=0.,iter_max=max_iter, iter_restart=iter_restart, verbose=self.verbose)
1646                   u=x[0]
1647                   q=x[1]
1648                   Bu=self.B(u)
1649                # now we check if |Bu| ~ 0 or more precise |Bu|_p  <= rtol * |v|_v
1650                nu=self.inner_v(u,u)
1651                p2=self.solve_prec(Bu)
1652                nBu=self.inner_p(p2,p2)
1653                if not nu>=0 and not nBu>=0: raise NegativeNorm,"negative norm."
1654                nu=math.sqrt(nu)
1655                nBu=math.sqrt(nBu)
1656                if self.verbose: print "saddle point solver: norm v= %e (Bv = %e)"%(nu,nBu)
1657                QTOL=atol+nu*rtol
1658                if nBu <= QTOL:
1659                    converged=True
1660                else:
1661                    ATOL=QTOL/nBu*ATOL*0.3
1662                    if self.verbose: print "correction step %s: tolerance reduced to %e."%(correction_step,ATOL)
1663                    converged=False
1664                correction_step+=1
1665                if correction_step>max_correction_steps:
1666                   raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1667             if self.verbose: print "saddle point solver: tolerance reached."
1668         return u,q
1669    
1670          #========================================================================
1671          def setTolerance(self,tolerance=1.e-4):
1672             """
1673             Sets the relative tolerance for (v,p).
1674    
1675             @param tolerance: tolerance to be used
1676             @type tolerance: non-negative C{float}
1677             """
1678             if tolerance<0:
1679                 raise ValueError,"tolerance must be positive."
1680             self.__rtol=tolerance
1681             self.setSubProblemTolerance()
1682    
1683          def getTolerance(self):
1684             """
1685             Returns the relative tolerance.
1686    
1687             @return: relative tolerance
1688             @rtype: C{float}
1689             """
1690             return self.__rtol
1691    
1692          def setAbsoluteTolerance(self,tolerance=0.):
1693             """
1694             Sets the absolute tolerance.
1695    
1696             @param tolerance: tolerance to be used
1697             @type tolerance: non-negative C{float}
1698             """
1699             if tolerance<0:
1700                 raise ValueError,"tolerance must be non-negative."
1701             self.__atol=tolerance
1702    
1703          def getAbsoluteTolerance(self):
1704             """
1705             Returns the absolute tolerance.
1706    
1707             @return: absolute tolerance
1708             @rtype: C{float}
1709             """
1710             return self.__atol
1711    
1712          def setSubProblemTolerance(self,rtol=None):
1713             """
1714             Sets the relative tolerance to solve the subproblem(s).
1715    
1716             @param rtol: relative tolerence
1717             @type rtol: positive C{float}
1718             """
1719             if rtol == None:
1720                  rtol=max(200.*util.EPSILON,self.getTolerance()**2)
1721             if rtol<=0:
1722                 raise ValueError,"tolerance must be positive."
1723             self.__sub_tol=rtol
1724    
1725          def getSubProblemTolerance(self):
1726             """
1727             Returns the subproblem reduction factor.
1728    
1729             @return: subproblem reduction factor
1730             @rtype: C{float}
1731             """
1732             return self.__sub_tol
1733    
1734    def MaskFromBoundaryTag(domain,*tags):
1735       """
1736       Creates a mask on the Solution(domain) function space where the value is
1737       one for samples that touch the boundary tagged by tags.
1738    
1739       Usage: m=MaskFromBoundaryTag(domain, "left", "right")
1740    
1741       @param domain: domain to be used
1742       @type domain: L{escript.Domain}
1743       @param tags: boundary tags
1744       @type tags: C{str}
1745       @return: a mask which marks samples that are touching the boundary tagged
1746                by any of the given tags
1747       @rtype: L{escript.Data} of rank 0
1748       """
1749       pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1750       d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
1751       for t in tags: d.setTaggedValue(t,1.)
1752       pde.setValue(y=d)
1753       return util.whereNonZero(pde.getRightHandSide())
1754    
1755    #==============================================================================
1756    class SaddlePointProblem(object):
1757       """
1758       This implements a solver for a saddle point problem
1759    
1760       M{f(u,p)=0}
1761       M{g(u)=0}
1762    
1763       for u and p. The problem is solved with an inexact Uszawa scheme for p:
1764    
1765       M{Q_f (u^{k+1}-u^{k}) = - f(u^{k},p^{k})}
1766       M{Q_g (p^{k+1}-p^{k}) =   g(u^{k+1})}
1767    
1768       where Q_f is an approximation of the Jacobian A_f of f with respect to u and
1769       Q_f is an approximation of A_g A_f^{-1} A_g with A_g is the Jacobian of g
1770       with respect to p. As a the construction of a 'proper' Q_g can be difficult,
1771       non-linear conjugate gradient method is applied to solve for p, so Q_g plays
1772       in fact the role of a preconditioner.
1773       """
1774       def __init__(self,verbose=False,*args):
1775           """
1776           Initializes the problem.
1777    
1778           @param verbose: if True, some information is printed in the process
1779           @type verbose: C{bool}
1780           @note: this method may be overwritten by a particular saddle point
1781                  problem
1782           """
1783           print "SaddlePointProblem should not be used anymore!"
1784           if not isinstance(verbose,bool):
1785                raise TypeError("verbose needs to be of type bool.")
1786           self.__verbose=verbose
1787           self.relaxation=1.
1788           DeprecationWarning("SaddlePointProblem should not be used anymore.")
1789    
1790       def trace(self,text):
1791           """
1792           Prints C{text} if verbose has been set.
1793    
1794           @param text: a text message
1795           @type text: C{str}
1796           """
1797           if self.__verbose: print "%s: %s"%(str(self),text)
1798    
1799       def solve_f(self,u,p,tol=1.e-8):
1800           """
1801           Solves
1802    
1803           A_f du = f(u,p)
1804    
1805           with tolerance C{tol} and returns du. A_f is Jacobian of f with respect
1806           to u.
1807    
1808           @param u: current approximation of u
1809           @type u: L{escript.Data}
1810           @param p: current approximation of p
1811           @type p: L{escript.Data}
1812           @param tol: tolerance expected for du
1813           @type tol: C{float}
1814           @return: increment du
1815           @rtype: L{escript.Data}
1816           @note: this method has to be overwritten by a particular saddle point
1817                  problem
1818           """
1819           pass
1820    
1821       def solve_g(self,u,tol=1.e-8):
1822           """
1823           Solves
1824    
1825           Q_g dp = g(u)
1826    
1827           where Q_g is a preconditioner for A_g A_f^{-1} A_g with A_g is the
1828           Jacobian of g with respect to p.
1829    
1830           @param u: current approximation of u
1831           @type u: L{escript.Data}
1832           @param tol: tolerance expected for dp
1833           @type tol: C{float}
1834           @return: increment dp
1835           @rtype: L{escript.Data}
1836           @note: this method has to be overwritten by a particular saddle point
1837                  problem
1838           """
1839           pass
1840    
1841       def inner(self,p0,p1):
1842           """
1843           Inner product of p0 and p1 approximating p. Typically this returns
1844           C{integrate(p0*p1)}.
1845           @return: inner product of p0 and p1
1846           @rtype: C{float}
1847           """
1848           pass
1849    
1850       subiter_max=3
1851       def solve(self,u0,p0,tolerance=1.e-6,tolerance_u=None,iter_max=100,accepted_reduction=0.995,relaxation=None):
1852            """
1853            Runs the solver.
1854    
1855            @param u0: initial guess for C{u}
1856            @type u0: L{esys.escript.Data}
1857            @param p0: initial guess for C{p}
1858            @type p0: L{esys.escript.Data}
1859            @param tolerance: tolerance for relative error in C{u} and C{p}
1860            @type tolerance: positive C{float}
1861            @param tolerance_u: tolerance for relative error in C{u} if different
1862                                from C{tolerance}
1863            @type tolerance_u: positive C{float}
1864            @param iter_max: maximum number of iteration steps
1865            @type iter_max: C{int}
1866            @param accepted_reduction: if the norm g cannot be reduced by
1867                                       C{accepted_reduction} backtracking to adapt
1868                                       the relaxation factor. If
1869                                       C{accepted_reduction=None} no backtracking
1870                                       is used.
1871            @type accepted_reduction: positive C{float} or C{None}
1872            @param relaxation: initial relaxation factor. If C{relaxation==None},
1873                               the last relaxation factor is used.
1874            @type relaxation: C{float} or C{None}
1875            """
1876            tol=1.e-2
1877            if tolerance_u==None: tolerance_u=tolerance
1878            if not relaxation==None: self.relaxation=relaxation
1879            if accepted_reduction ==None:
1880                  angle_limit=0.
1881            elif accepted_reduction>=1.:
1882                  angle_limit=0.
1883            else:
1884                  angle_limit=util.sqrt(1-accepted_reduction**2)
1885            self.iter=0
1886            u=u0
1887            p=p0
1888            #
1889            #   initialize things:
1890            #
1891            converged=False
1892            #
1893            #  start loop:
1894            #
1895            #  initial search direction is g
1896            #
1897            while not converged :
1898                if self.iter>iter_max:
1899                    raise ArithmeticError("no convergence after %s steps."%self.iter)
1900                f_new=self.solve_f(u,p,tol)
1901                norm_f_new = util.Lsup(f_new)
1902                u_new=u-f_new
1903                g_new=self.solve_g(u_new,tol)
1904                self.iter+=1
1905                norm_g_new = util.sqrt(self.inner(g_new,g_new))
1906                if norm_f_new==0. and norm_g_new==0.: return u, p
1907                if self.iter>1 and not accepted_reduction==None:
1908                   #
1909                   #   did we manage to reduce the norm of G? I
1910                   #   if not we start a backtracking procedure
1911                   #
1912                   # print "new/old norm = ",norm_g_new, norm_g, norm_g_new/norm_g
1913                   if norm_g_new > accepted_reduction * norm_g:
1914                      sub_iter=0
1915                      s=self.relaxation
1916                      d=g
1917                      g_last=g
1918                      self.trace("    start substepping: f = %s, g = %s, relaxation = %s."%(norm_f_new, norm_g_new, s))
1919                      while sub_iter < self.subiter_max and  norm_g_new > accepted_reduction * norm_g:
1920                         dg= g_new-g_last
1921                         norm_dg=abs(util.sqrt(self.inner(dg,dg))/self.relaxation)
1922                         rad=self.inner(g_new,dg)/self.relaxation
1923                         # print "   ",sub_iter,": rad, norm_dg:",abs(rad), norm_dg*norm_g_new * angle_limit
1924                         # print "   ",sub_iter,": rad, norm_dg:",rad, norm_dg, norm_g_new, norm_g
1925                         if abs(rad) < norm_dg*norm_g_new * angle_limit:
1926                             if sub_iter>0: self.trace("    no further improvements expected from backtracking.")
1927                             break
1928                         r=self.relaxation
1929                         self.relaxation= - rad/norm_dg**2
1930                         s+=self.relaxation
1931                         #####
1932                         # a=g_new+self.relaxation*dg/r
1933                         # print "predicted new norm = ",util.sqrt(self.inner(a,a)),util.sqrt(self.inner(g_new,g_new)), self.relaxation
1934                         #####
1935                         g_last=g_new
1936                         p+=self.relaxation*d
1937                         f_new=self.solve_f(u,p,tol)
1938                         u_new=u-f_new
1939                         g_new=self.solve_g(u_new,tol)
1940                         self.iter+=1
1941                         norm_f_new = util.Lsup(f_new)
1942                         norm_g_new = util.sqrt(self.inner(g_new,g_new))
1943                         # print "   ",sub_iter," new g norm",norm_g_new
1944                         self.trace("    %s th sub-step: f = %s, g = %s, relaxation = %s."%(sub_iter, norm_f_new, norm_g_new, s))
1945                         #
1946                         #   can we expect reduction of g?
1947                         #
1948                         # u_last=u_new
1949                         sub_iter+=1
1950                      self.relaxation=s
1951                #
1952                #  check for convergence:
1953                #
1954                norm_u_new = util.Lsup(u_new)
1955                p_new=p+self.relaxation*g_new
1956                norm_p_new = util.sqrt(self.inner(p_new,p_new))
1957                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))
1958    
1959                if self.iter>1:
1960                   dg2=g_new-g
1961                   df2=f_new-f
1962                   norm_dg2=util.sqrt(self.inner(dg2,dg2))
1963                   norm_df2=util.Lsup(df2)
1964                   # print norm_g_new, norm_g, norm_dg, norm_p, tolerance
1965                   tol_eq_g=tolerance*norm_dg2/(norm_g*abs(self.relaxation))*norm_p_new
1966                   tol_eq_f=tolerance_u*norm_df2/norm_f*norm_u_new
1967                   if norm_g_new <= tol_eq_g and norm_f_new <= tol_eq_f:
1968                       converged=True
1969                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
1970            self.trace("convergence after %s steps."%self.iter)
1971            return u,p
1972    

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

  ViewVC Help
Powered by ViewVC 1.1.26