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

Legend:
Removed from v.637  
changed lines
  Added in v.1878

  ViewVC Help
Powered by ViewVC 1.1.26