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

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