/[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 1122 by gross, Tue May 1 03:21:04 2007 UTC revision 1809 by ksteube, Thu Sep 25 06:43:44 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 18  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    """    """
# Line 393  class Locator: Line 410  class Locator:
410          else:          else:
411             return data             return data
412    
413    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    
564    ############################
565    # Added by Artak
566    #################################3
567    
568    #Apply a sequence of k Givens rotations, used within gmres codes
569    # vrot=givapp(c, s, vin, k)
570    def givapp(c,s,vin):
571        vrot=vin # warning: vin is altered!!!!
572        if isinstance(c,float):
573            vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
574        else:
575            for i in range(len(c)):
576                w1=c[i]*vrot[i]-s[i]*vrot[i+1]
577            w2=s[i]*vrot[i]+c[i]*vrot[i+1]
578                vrot[i:i+2]=w1,w2
579        return vrot
580    
581    ##############################################
582    def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
583    ################################################
584       m=iter_restart
585       iter=0
586       xc=x
587       while True:
588          if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
589          xc,stopped=GMRESm(b*1, Aprod, Msolve, bilinearform, stoppingcriterium, x=xc*1, iter_max=iter_max-iter, iter_restart=m)
590          if stopped: break
591          iter+=iter_restart    
592       return xc
593    
594    def GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
595       iter=0
596       r=Msolve(b)
597       r_dot_r = bilinearform(r, r)
598       if r_dot_r<0: raise NegativeNorm,"negative norm."
599       norm_b=math.sqrt(r_dot_r)
600    
601       if x==None:
602          x=0*b
603       else:
604          r=Msolve(b-Aprod(x))
605          r_dot_r = bilinearform(r, r)
606          if r_dot_r<0: raise NegativeNorm,"negative norm."
607      
608       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
609       c=numarray.zeros(iter_restart,numarray.Float64)
610       s=numarray.zeros(iter_restart,numarray.Float64)
611       g=numarray.zeros(iter_restart,numarray.Float64)
612       v=[]
613    
614       rho=math.sqrt(r_dot_r)
615      
616       v.append(r/rho)
617       g[0]=rho
618    
619       while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
620    
621        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
622    
623        p=Msolve(Aprod(v[iter]))
624    
625        v.append(p)
626    
627        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
628    
629    # Modified Gram-Schmidt
630        for j in range(iter+1):
631          h[j][iter]=bilinearform(v[j],v[iter+1])  
632          v[iter+1]-=h[j][iter]*v[j]
633          
634        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
635        v_norm2=h[iter+1][iter]
636    
637    # Reorthogonalize if needed
638        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
639         for j in range(iter+1):  
640            hr=bilinearform(v[j],v[iter+1])
641                h[j][iter]=h[j][iter]+hr
642                v[iter+1] -= hr*v[j]
643    
644         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
645         h[iter+1][iter]=v_norm2
646    
647    #   watch out for happy breakdown
648            if not v_norm2 == 0:
649             v[iter+1]=v[iter+1]/h[iter+1][iter]
650    
651    #   Form and store the information for the new Givens rotation
652        if iter > 0 :
653            hhat=numarray.zeros(iter+1,numarray.Float64)
654            for i in range(iter+1) : hhat[i]=h[i][iter]
655            hhat=givapp(c[0:iter],s[0:iter],hhat);
656                for i in range(iter+1) : h[i][iter]=hhat[i]
657    
658        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
659    
660        if mu!=0 :
661            c[iter]=h[iter][iter]/mu
662            s[iter]=-h[iter+1][iter]/mu
663            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
664            h[iter+1][iter]=0.0
665            g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
666    
667    # Update the residual norm
668                  
669            rho=abs(g[iter+1])
670        iter+=1
671    
672    # At this point either iter > iter_max or rho < tol.
673    # It's time to compute x and leave.        
674    
675       if iter > 0 :
676         y=numarray.zeros(iter,numarray.Float64)    
677         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
678         if iter > 1 :  
679            i=iter-2  
680            while i>=0 :
681              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
682              i=i-1
683         xhat=v[iter-1]*y[iter-1]
684         for i in range(iter-1):
685        xhat += v[i]*y[i]
686       else : xhat=v[0]
687    
688       x += xhat
689       if iter<iter_restart-1:
690          stopped=True
691       else:
692          stopped=False
693    
694       return x,stopped
695    
696    
697    ######################################################
698    def dirder(x, w, bilinearform, Aprod, Msolve, f0, b ):
699    ######################################################
700    
701    # DIRDER estimates the directional derivative of a function F.
702    
703    
704    # Hardwired difference increment.
705    #
706      epsnew = 1.0e-07
707    #
708    #  Scale the step.
709    #
710      norm_w=math.sqrt(bilinearform(w,w))
711      if norm_w== 0.0:
712        return x/x
713    
714      epsnew = epsnew / norm_w
715    
716      if norm_w > 0.0:
717        epsnew = epsnew * math.sqrt(bilinearform(x,x))
718    #
719    #  DEL and F1 could share the same space if storage
720    #  is more important than clarity.
721    #
722    
723      DEL = x + epsnew * w
724      f1 = -Msolve(Aprod(DEL))
725      z = ( f1 - f0 ) / epsnew
726      return z
727    
728    ######################################################
729    def FDGMRES(f0, Aprod, Msolve, bilinearform, stoppingcriterium, xc=None, x=None, iter_max=100, iter_restart=20,TOL=None):
730    ######################################################
731       b=-f0
732       b_dot_b = bilinearform(b, b)
733       if b_dot_b<0: raise NegativeNorm,"negative norm."
734       norm_b=math.sqrt(b_dot_b)
735    
736       r=b
737    
738       if x==None:
739          x=0*f0
740       else:
741          r=-dirder(xc,x,bilinearform,Aprod,Msolve,f0,b)-f0  
742          
743       r_dot_r = bilinearform(r, r)
744       if r_dot_r<0: raise NegativeNorm,"negative norm."
745      
746       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
747       c=numarray.zeros(iter_restart,numarray.Float64)
748       s=numarray.zeros(iter_restart,numarray.Float64)
749       g=numarray.zeros(iter_restart,numarray.Float64)
750       v=[]
751    
752       rho=math.sqrt(r_dot_r)
753      
754       v.append(r/rho)
755       g[0]=rho
756       iter=0
757    
758       while not (stoppingcriterium(rho,norm_b,solver="FDGMRES",TOL=TOL) or iter==iter_restart-1):
759    
760        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
761    
762        
763            p=dirder(xc, v[iter], bilinearform,Aprod,Msolve,f0,b)
764    
765        v.append(p)
766    
767        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
768    
769    # Modified Gram-Schmidt
770        for j in range(iter+1):
771          h[j][iter]=bilinearform(v[j],v[iter+1])  
772          v[iter+1]+=(-1.)*h[j][iter]*v[j]
773          
774        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
775        v_norm2=h[iter+1][iter]
776    
777    
778    # Reorthogonalize if needed
779        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
780         for j in range(iter+1):
781            hr=bilinearform(v[j],v[iter+1])
782                h[j][iter]=h[j][iter]+hr #vhat
783                v[iter+1] +=(-1.)*hr*v[j]
784    
785         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
786         h[iter+1][iter]=v_norm2
787    
788    #   watch out for happy breakdown
789            if v_norm2 != 0:
790             v[iter+1]=v[iter+1]/h[iter+1][iter]
791    
792    #   Form and store the information for the new Givens rotation
793        if iter > 0 :
794            hhat=[]
795            for i in range(iter+1) : hhat.append(h[i][iter])
796            hhat=givapp(c[0:iter],s[0:iter],hhat);
797                for i in range(iter+1) : h[i][iter]=hhat[i]
798    
799        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
800        if mu!=0 :
801            c[iter]=h[iter][iter]/mu
802            s[iter]=-h[iter+1][iter]/mu
803            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
804            h[iter+1][iter]=0.0
805            g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
806    
807    # Update the residual norm
808            rho=abs(g[iter+1])
809        iter+=1
810    
811       if iter > 0 :
812         y=numarray.zeros(iter,numarray.Float64)    
813         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
814         if iter > 1 :  
815            i=iter-2  
816            while i>=0 :
817              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
818              i=i-1
819         xhat=v[iter-1]*y[iter-1]
820         for i in range(iter-1):
821        xhat += v[i]*y[i]
822       else : xhat=v[0]
823        
824       x += xhat
825       if iter<iter_restart-1:
826          stopped=True
827       else:
828          stopped=False
829    
830       return x,stopped
831    
832    #################################################
833    def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
834    #################################################
835        #
836        #  minres solves the system of linear equations Ax = b
837        #  where A is a symmetric matrix (possibly indefinite or singular)
838        #  and b is a given vector.
839        #  
840        #  "A" may be a dense or sparse matrix (preferably sparse!)
841        #  or the name of a function such that
842        #               y = A(x)
843        #  returns the product y = Ax for any given vector x.
844        #
845        #  "M" defines a positive-definite preconditioner M = C C'.
846        #  "M" may be a dense or sparse matrix (preferably sparse!)
847        #  or the name of a function such that
848        #  solves the system My = x for any given vector x.
849        #
850        #
851        
852        #------------------------------------------------------------------
853        # Set up y and v for the first Lanczos vector v1.
854        # y  =  beta1 P' v1,  where  P = C**(-1).
855        # v is really P' v1.
856        #------------------------------------------------------------------
857        if x==None:
858          x=0*b
859        else:
860          b += (-1)*Aprod(x)
861    
862        r1    = b
863        y = Msolve(b)
864        beta1 = bilinearform(y,b)
865    
866        if beta1< 0: raise NegativeNorm,"negative norm."
867    
868        #  If b = 0 exactly, stop with x = 0.
869        if beta1==0: return x*0.
870    
871        if beta1> 0:
872          beta1  = math.sqrt(beta1)      
873    
874        #------------------------------------------------------------------
875        # Initialize quantities.
876        # ------------------------------------------------------------------
877        iter   = 0
878        Anorm = 0
879        ynorm = 0
880        oldb   = 0
881        beta   = beta1
882        dbar   = 0
883        epsln  = 0
884        phibar = beta1
885        rhs1   = beta1
886        rhs2   = 0
887        rnorm  = phibar
888        tnorm2 = 0
889        ynorm2 = 0
890        cs     = -1
891        sn     = 0
892        w      = b*0.
893        w2     = b*0.
894        r2     = r1
895        eps    = 0.0001
896    
897        #---------------------------------------------------------------------
898        # Main iteration loop.
899        # --------------------------------------------------------------------
900        while not stoppingcriterium(rnorm,Anorm*ynorm,'MINRES'):    #  checks ||r|| < (||A|| ||x||) * TOL
901    
902        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
903            iter    = iter  +  1
904    
905            #-----------------------------------------------------------------
906            # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
907            # The general iteration is similar to the case k = 1 with v0 = 0:
908            #
909            #   p1      = Operator * v1  -  beta1 * v0,
910            #   alpha1  = v1'p1,
911            #   q2      = p2  -  alpha1 * v1,
912            #   beta2^2 = q2'q2,
913            #   v2      = (1/beta2) q2.
914            #
915            # Again, y = betak P vk,  where  P = C**(-1).
916            #-----------------------------------------------------------------
917            s = 1/beta                 # Normalize previous vector (in y).
918            v = s*y                    # v = vk if P = I
919        
920            y      = Aprod(v)
921        
922            if iter >= 2:
923              y = y - (beta/oldb)*r1
924    
925            alfa   = bilinearform(v,y)              # alphak
926            y      += (- alfa/beta)*r2
927            r1     = r2
928            r2     = y
929            y = Msolve(r2)
930            oldb   = beta                           # oldb = betak
931            beta   = bilinearform(y,r2)             # beta = betak+1^2
932            if beta < 0: raise NegativeNorm,"negative norm."
933    
934            beta   = math.sqrt( beta )
935            tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
936            
937            if iter==1:                 # Initialize a few things.
938              gmax   = abs( alfa )      # alpha1
939              gmin   = gmax             # alpha1
940    
941            # Apply previous rotation Qk-1 to get
942            #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
943            #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
944        
945            oldeps = epsln
946            delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
947            gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
948            epsln  =               sn * beta  # epsln2 = 0         epslnk+1
949            dbar   =            -  cs * beta  # dbar 2 = beta2     dbar k+1
950    
951            # Compute the next plane rotation Qk
952    
953            gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
954            gamma  = max(gamma,eps)
955            cs     = gbar / gamma             # ck
956            sn     = beta / gamma             # sk
957            phi    = cs * phibar              # phik
958            phibar = sn * phibar              # phibark+1
959    
960            # Update  x.
961    
962            denom = 1/gamma
963            w1    = w2
964            w2    = w
965            w     = (v - oldeps*w1 - delta*w2) * denom
966            x     +=  phi*w
967    
968            # Go round again.
969    
970            gmax   = max(gmax,gamma)
971            gmin   = min(gmin,gamma)
972            z      = rhs1 / gamma
973            ynorm2 = z*z  +  ynorm2
974            rhs1   = rhs2 -  delta*z
975            rhs2   =      -  epsln*z
976    
977            # Estimate various norms and test for convergence.
978    
979            Anorm  = math.sqrt( tnorm2 )
980            ynorm  = math.sqrt( ynorm2 )
981    
982            rnorm  = phibar
983    
984        return x
985    
986    ######################################    
987    def NewtonGMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium,x=None, iter_max=100,iter_restart=20,atol=1.e-2,rtol=1.e-4):
988    #####################################
989        gamma=.9
990        lmaxit=100
991        etamax=.5
992    
993        n = 1 #len(x)
994        iter=0
995        
996        # evaluate f at the initial iterate
997        # compute the stop tolerance
998        #
999        r=b
1000        if x==None:
1001          x=0*b
1002        else:
1003          r += (-1)*Aprod(x)
1004    
1005        f0=-Msolve(r)
1006        fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
1007        fnrmo=1
1008        stop_tol=atol + rtol*fnrm
1009        #
1010        # main iteration loop
1011        #
1012        while not stoppingcriterium(fnrm*1,stop_tol,'NewtonGMRES',TOL=1.):
1013    
1014                if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1015            #
1016            # keep track of the ratio (rat = fnrm/frnmo)
1017            # of successive residual norms and
1018            # the iteration counter (iter)
1019            #
1020            #rat=fnrm/fnrmo
1021            fnrmo=fnrm
1022            iter+=1
1023            #
1024                # compute the step using a GMRES(m) routine especially designed for this purpose
1025            #
1026                initer=0
1027                while True:
1028                   xc,stopped=FDGMRES(f0*1, Aprod, Msolve, bilinearform, stoppingcriterium, xc=x*1, iter_max=lmaxit-initer, iter_restart=iter_restart, TOL=etamax)
1029                   if stopped: break
1030                   initer+=iter_restart
1031            x+=xc
1032            f0=-Msolve(Aprod(x))
1033            fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
1034            rat=fnrm/fnrmo
1035    
1036    
1037        #   adjust eta
1038        #
1039            if etamax > 0:
1040                etaold=etamax
1041                etanew=gamma*rat*rat
1042                if gamma*etaold*etaold > .1 :
1043                    etanew=max(etanew,gamma*etaold*etaold)
1044                etamax=min(etanew,etamax)
1045                etamax=max(etamax,.5*stop_tol/fnrm)
1046        return x
1047    
1048    def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
1049    
1050    # TFQMR solver for linear systems
1051    #
1052    #
1053    # initialization
1054    #
1055      errtol = math.sqrt(bilinearform(b,b))
1056      norm_b=errtol
1057      kmax  = iter_max
1058      error = []
1059    
1060      if math.sqrt(bilinearform(x,x)) != 0.0:
1061        r = b - Aprod(x)
1062      else:
1063        r = b
1064    
1065      r=Msolve(r)
1066    
1067      u1=0
1068      u2=0
1069      y1=0
1070      y2=0
1071    
1072      w = r
1073      y1 = r
1074      iter = 0
1075      d = 0
1076      
1077      v = Msolve(Aprod(y1))
1078      u1 = v
1079      
1080      theta = 0.0;
1081      eta = 0.0;
1082      tau = math.sqrt(bilinearform(r,r))
1083      error = [ error, tau ]
1084      rho = tau * tau
1085      m=1
1086    #
1087    #  TFQMR iteration
1088    #
1089    #  while ( iter < kmax-1 ):
1090      
1091      while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):
1092        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1093    
1094        sigma = bilinearform(r,v)
1095    
1096        if ( sigma == 0.0 ):
1097          raise 'TFQMR breakdown, sigma=0'
1098        
1099    
1100        alpha = rho / sigma
1101    
1102        for j in range(2):
1103    #
1104    #   Compute y2 and u2 only if you have to
1105    #
1106          if ( j == 1 ):
1107            y2 = y1 - alpha * v
1108            u2 = Msolve(Aprod(y2))
1109    
1110          m = 2 * (iter+1) - 2 + (j+1)
1111          if j==0:
1112             w = w - alpha * u1
1113             d = y1 + ( theta * theta * eta / alpha ) * d
1114          if j==1:
1115             w = w - alpha * u2
1116             d = y2 + ( theta * theta * eta / alpha ) * d
1117    
1118          theta = math.sqrt(bilinearform(w,w))/ tau
1119          c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1120          tau = tau * theta * c
1121          eta = c * c * alpha
1122          x = x + eta * d
1123    #
1124    #  Try to terminate the iteration at each pass through the loop
1125    #
1126         # if ( tau * math.sqrt ( m + 1 ) <= errtol ):
1127         #   error = [ error, tau ]
1128         #   total_iters = iter
1129         #   break
1130          
1131    
1132        if ( rho == 0.0 ):
1133          raise 'TFQMR breakdown, rho=0'
1134        
1135    
1136        rhon = bilinearform(r,w)
1137        beta = rhon / rho;
1138        rho = rhon;
1139        y1 = w + beta * y2;
1140        u1 = Msolve(Aprod(y1))
1141        v = u1 + beta * ( u2 + beta * v )
1142        error = [ error, tau ]
1143        total_iters = iter
1144        
1145        iter = iter + 1
1146    
1147      return x
1148    
1149    
1150    #############################################
1151    
1152    class ArithmeticTuple(object):
1153       """
1154       tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
1155    
1156       example of usage:
1157    
1158       from esys.escript import Data
1159       from numarray import array
1160       a=Data(...)
1161       b=array([1.,4.])
1162       x=ArithmeticTuple(a,b)
1163       y=5.*x
1164    
1165       """
1166       def __init__(self,*args):
1167           """
1168           initialize object with elements args.
1169    
1170           @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
1171           """
1172           self.__items=list(args)
1173    
1174       def __len__(self):
1175           """
1176           number of items
1177    
1178           @return: number of items
1179           @rtype: C{int}
1180           """
1181           return len(self.__items)
1182    
1183       def __getitem__(self,index):
1184           """
1185           get an item
1186    
1187           @param index: item to be returned
1188           @type index: C{int}
1189           @return: item with index C{index}
1190           """
1191           return self.__items.__getitem__(index)
1192    
1193       def __mul__(self,other):
1194           """
1195           scaling from the right
1196    
1197           @param other: scaling factor
1198           @type other: C{float}
1199           @return: itemwise self*other
1200           @rtype: L{ArithmeticTuple}
1201           """
1202           out=[]
1203           other=1.*other
1204           if isinstance(other,float):
1205        for i in range(len(self)):
1206               out.append(self[i]*other)
1207           else:
1208            for i in range(len(self)):
1209               out.append(self[i]*other[i])
1210           return ArithmeticTuple(*tuple(out))
1211    
1212       def __rmul__(self,other):
1213           """
1214           scaling from the left
1215    
1216           @param other: scaling factor
1217           @type other: C{float}
1218           @return: itemwise other*self
1219           @rtype: L{ArithmeticTuple}
1220           """
1221           out=[]
1222           other=1.*other
1223           if isinstance(other,float):
1224        for i in range(len(self)):
1225               out.append(other*self[i])
1226           else:
1227            for i in range(len(self)):
1228               out.append(other[i]*self[i])
1229           return ArithmeticTuple(*tuple(out))
1230    
1231    #########################
1232    # Added by Artak
1233    #########################
1234       def __div__(self,other):
1235           """
1236           dividing from the right
1237    
1238           @param other: scaling factor
1239           @type other: C{float}
1240           @return: itemwise self/other
1241           @rtype: L{ArithmeticTuple}
1242           """
1243           out=[]
1244           other=1.*other
1245           if isinstance(other,float):
1246        for i in range(len(self)):
1247               out.append(self[i]*(1./other))
1248           else:
1249            for i in range(len(self)):
1250               out.append(self[i]*(1./other[i]))
1251           return ArithmeticTuple(*tuple(out))
1252    
1253       def __rdiv__(self,other):
1254           """
1255           dividing from the left
1256    
1257           @param other: scaling factor
1258           @type other: C{float}
1259           @return: itemwise other/self
1260           @rtype: L{ArithmeticTuple}
1261           """
1262           out=[]
1263           other=1.*other
1264           if isinstance(other,float):
1265            for i in range(len(self)):
1266               out.append(other*(1./self[i]))
1267           else:
1268            for i in range(len(self)):
1269               out.append(other[i]*(1./self[i]))
1270           return ArithmeticTuple(*tuple(out))
1271      
1272    ##########################################33
1273    
1274       def __iadd__(self,other):
1275           """
1276           in-place add of other to self
1277    
1278           @param other: increment
1279           @type other: C{ArithmeticTuple}
1280           """
1281           if len(self) != len(other):
1282               raise ValueError,"tuple length must match."
1283           for i in range(len(self)):
1284               self.__items[i]+=other[i]
1285           return self
1286    
1287       def __add__(self,other):
1288           """
1289           add other to self
1290    
1291           @param other: increment
1292           @type other: C{ArithmeticTuple}
1293           """
1294    #      if not isinstance(other,float):
1295           if len(self) != len(other):
1296              raise ValueError,"tuple length must match."
1297           for i in range(len(self)):
1298              self.__items[i]+=other[i]
1299    #       else:
1300    #        for i in range(len(self)):
1301    #           self.__items[i]+=other
1302    
1303           return self
1304    
1305       def __sub__(self,other):
1306           """
1307           subtract other from self
1308    
1309           @param other: increment
1310           @type other: C{ArithmeticTuple}
1311           """
1312           if len(self) != len(other):
1313               raise ValueError,"tuple length must match."
1314           for i in range(len(self)):
1315               self.__items[i]-=other[i]
1316           return self
1317      
1318       def __isub__(self,other):
1319           """
1320           subtract other from self
1321    
1322           @param other: increment
1323           @type other: C{ArithmeticTuple}
1324           """
1325           if len(self) != len(other):
1326               raise ValueError,"tuple length must match."
1327           for i in range(len(self)):
1328               self.__items[i]-=other[i]
1329           return self
1330    
1331       def __neg__(self):
1332           """
1333           negate
1334    
1335           """
1336           for i in range(len(self)):
1337               self.__items[i]=-self.__items[i]
1338           return self
1339    
1340    
1341    class HomogeneousSaddlePointProblem(object):
1342          """
1343          This provides a framwork for solving homogeneous saddle point problem of the form
1344    
1345                 Av+B^*p=f
1346                 Bv    =0
1347    
1348          for the unknowns v and p and given operators A and B and given right hand side f.
1349          B^* is the adjoint operator of B is the given inner product.
1350    
1351          """
1352          def __init__(self,**kwargs):
1353            self.setTolerance()
1354            self.setToleranceReductionFactor()
1355    
1356          def initialize(self):
1357            """
1358            initialize the problem (overwrite)
1359            """
1360            pass
1361          def B(self,v):
1362             """
1363             returns Bv (overwrite)
1364             @rtype: equal to the type of p
1365    
1366             @note: boundary conditions on p should be zero!
1367             """
1368             pass
1369    
1370          def inner(self,p0,p1):
1371             """
1372             returns inner product of two element p0 and p1  (overwrite)
1373            
1374             @type p0: equal to the type of p
1375             @type p1: equal to the type of p
1376             @rtype: C{float}
1377    
1378             @rtype: equal to the type of p
1379             """
1380             pass
1381    
1382          def solve_A(self,u,p):
1383             """
1384             solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
1385    
1386             @rtype: equal to the type of v
1387             @note: boundary conditions on v should be zero!
1388             """
1389             pass
1390    
1391          def solve_prec(self,p):
1392             """
1393             provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
1394    
1395             @rtype: equal to the type of p
1396             """
1397             pass
1398    
1399          def stoppingcriterium(self,Bv,v,p):
1400             """
1401             returns a True if iteration is terminated. (overwrite)
1402    
1403             @rtype: C{bool}
1404             """
1405             pass
1406                
1407          def __inner(self,p,r):
1408             return self.inner(p,r[1])
1409    
1410          def __inner_p(self,p1,p2):
1411             return self.inner(p1,p2)
1412          
1413          def __inner_a(self,a1,a2):
1414             return self.inner_a(a1,a2)
1415    
1416          def __inner_a1(self,a1,a2):
1417             return self.inner(a1[1],a2[1])
1418    
1419          def __stoppingcriterium(self,norm_r,r,p):
1420              return self.stoppingcriterium(r[1],r[0],p)
1421    
1422          def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
1423              return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)
1424    
1425          def setTolerance(self,tolerance=1.e-8):
1426                  self.__tol=tolerance
1427          def getTolerance(self):
1428                  return self.__tol
1429          def setToleranceReductionFactor(self,reduction=0.01):
1430                  self.__reduction=reduction
1431          def getSubProblemTolerance(self):
1432                  return self.__reduction*self.getTolerance()
1433    
1434          def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):
1435                  """
1436                  solves the saddle point problem using initial guesses v and p.
1437    
1438                  @param max_iter: maximum number of iteration steps.
1439                  """
1440                  self.verbose=verbose
1441                  self.show_details=show_details and self.verbose
1442    
1443                  # assume p is known: then v=A^-1(f-B^*p)
1444                  # which leads to BA^-1B^*p = BA^-1f  
1445    
1446              # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)      
1447              self.__z=v+self.solve_A(v,p*0)
1448                  Bz=self.B(self.__z)
1449                  #
1450              #   solve BA^-1B^*p = Bz
1451                  #
1452                  #
1453                  #
1454                  self.iter=0
1455              if solver=='GMRES':      
1456                    if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1457                    p=GMRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)
1458                    # solve Au=f-B^*p
1459                    #       A(u-v)=f-B^*p-Av
1460                    #       u=v+(u-v)
1461            u=v+self.solve_A(v,p)
1462    
1463              if solver=='NewtonGMRES':    
1464                    if self.verbose: print "enter NewtonGMRES method (iter_max=%s)"%max_iter
1465                    p=NewtonGMRES(Bz,self.__Aprod_Newton,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,atol=0,rtol=self.getTolerance())
1466                    # solve Au=f-B^*p
1467                    #       A(u-v)=f-B^*p-Av
1468                    #       u=v+(u-v)
1469            u=v+self.solve_A(v,p)
1470                    
1471    
1472              if solver=='TFQMR':      
1473                    if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter
1474                    p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1475                    # solve Au=f-B^*p
1476                    #       A(u-v)=f-B^*p-Av
1477                    #       u=v+(u-v)
1478            u=v+self.solve_A(v,p)
1479    
1480              if solver=='MINRES':      
1481                    if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1482                    p=MINRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1483                    # solve Au=f-B^*p
1484                    #       A(u-v)=f-B^*p-Av
1485                    #       u=v+(u-v)
1486            u=v+self.solve_A(v,p)
1487                  
1488              if solver=='GMRESC':      
1489                    if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter
1490                    p0=self.solve_prec1(Bz)
1491                #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))
1492                    #p-=alfa
1493                    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)
1494                    #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())
1495    
1496                    # solve Au=f-B^*p
1497                    #       A(u-v)=f-B^*p-Av
1498                    #       u=v+(u-v)
1499                p=x[1]
1500            u=v+self.solve_A(v,p)      
1501            #p=x[1]
1502            #u=x[0]
1503    
1504                  if solver=='PCG':
1505                    #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1506                    #
1507                    #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1508                    #                           A(v-z)= f -Az - B^*p (v-z=0 on fixed_u_mask)
1509                    if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1510                    p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1511                u=r[0]  
1512                    # print "DIFF=",util.integrate(p)
1513    
1514                  # print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1515    
1516              return u,p
1517    
1518          def __Msolve(self,r):
1519              return self.solve_prec(r[1])
1520    
1521          def __Msolve2(self,r):
1522              return self.solve_prec(r*1)
1523    
1524          def __Mempty(self,r):
1525              return r
1526    
1527    
1528          def __Aprod(self,p):
1529              # return BA^-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 __Anext(self,x):
1535              # return next v,p
1536              #solve Av =-B^*p as Av =f-Az-B^*p
1537    
1538          pc=x[1]
1539              v=self.solve_A(self.__z,-pc)
1540          p=self.solve_prec1(self.B(v))
1541    
1542              return ArithmeticTuple(v,p)
1543    
1544    
1545          def __Aprod2(self,p):
1546              # return BA^-1B*p
1547              #solve Av =B^*p as Av =f-Az-B^*(-p)
1548          v=self.solve_A(self.__z,-p)
1549              return self.B(v)
1550    
1551          def __Aprod_Newton(self,p):
1552              # return BA^-1B*p - Bz
1553              #solve Av =-B^*p as Av =f-Az-B^*p
1554          v=self.solve_A(self.__z,-p)
1555              return self.B(v-self.__z)
1556    
1557          def __Aprod_Newton2(self,x):
1558              # return BA^-1B*p - Bz
1559              #solve Av =-B^*p as Av =f-Az-B^*p
1560              pc=x[1]
1561          v=self.solve_A(self.__z,-pc)
1562              p=self.solve_prec1(self.B(v-self.__z))
1563              return ArithmeticTuple(v,p)
1564    
1565  class SaddlePointProblem(object):  class SaddlePointProblem(object):
1566     """     """
1567     This implements a solver for a saddlepoint problem     This implements a solver for a saddlepoint problem
# Line 418  class SaddlePointProblem(object): Line 1587  class SaddlePointProblem(object):
1587         @type verbose: C{bool}         @type verbose: C{bool}
1588         @note: this method may be overwritten by a particular saddle point problem         @note: this method may be overwritten by a particular saddle point problem
1589         """         """
1590           print "SaddlePointProblem should not be used anymore!"
1591         if not isinstance(verbose,bool):         if not isinstance(verbose,bool):
1592              raise TypeError("verbose needs to be of type bool.")              raise TypeError("verbose needs to be of type bool.")
1593         self.__verbose=verbose         self.__verbose=verbose
# Line 580  class SaddlePointProblem(object): Line 1750  class SaddlePointProblem(object):
1750              norm_u_new = util.Lsup(u_new)              norm_u_new = util.Lsup(u_new)
1751              p_new=p+self.relaxation*g_new              p_new=p+self.relaxation*g_new
1752              norm_p_new = util.sqrt(self.inner(p_new,p_new))              norm_p_new = util.sqrt(self.inner(p_new,p_new))
1753              self.trace("%s th step: f/u = %s, g/p = %s, relaxation = %s."%(self.iter,norm_f_new/norm_u_new, norm_g_new/norm_p_new, self.relaxation))              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))
1754    
1755              if self.iter>1:              if self.iter>1:
1756                 dg2=g_new-g                 dg2=g_new-g
# Line 622  class SaddlePointProblem(object): Line 1792  class SaddlePointProblem(object):
1792  #  #
1793  #      return u,p  #      return u,p
1794                        
1795  # vim: expandtab shiftwidth=4:  def MaskFromBoundaryTag(domain,*tags):
1796       """
1797       create a mask on the Solution(domain) function space which one for samples
1798       that touch the boundary tagged by tags.
1799    
1800       usage: m=MaskFromBoundaryTag(domain,"left", "right")
1801    
1802       @param domain: a given domain
1803       @type domain: L{escript.Domain}
1804       @param tags: boundray tags
1805       @type tags: C{str}
1806       @return: a mask which marks samples that are touching the boundary tagged by any of the given tags.
1807       @rtype: L{escript.Data} of rank 0
1808       """
1809       pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1810       d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
1811       for t in tags: d.setTaggedValue(t,1.)
1812       pde.setValue(y=d)
1813       return util.whereNonZero(pde.getRightHandSide())

Legend:
Removed from v.1122  
changed lines
  Added in v.1809

  ViewVC Help
Powered by ViewVC 1.1.26