/[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 1331 by gross, Tue Oct 23 00:42:15 2007 UTC revision 1482 by artak, Wed Apr 9 02:29:47 2008 UTC
# Line 48  import numarray Line 48  import numarray
48  import util  import util
49  import math  import math
50    
51    ##### Added by Artak
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 472  class IterationHistory(object): Line 477  class IterationHistory(object):
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):
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           self.history.append(norm_r)
494           if self.verbose: print "iter: %s:  norm(r) = %e"%(len(self.history)-1, self.history[-1])
495           return self.history[-1]<=self.tolerance * norm_b
496    
497  def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
498     """     """
499     Solver for     Solver for
# Line 539  type like argument C{x}. Line 561  type like argument C{x}.
561    
562     return x,r     return x,r
563    
564    
565    ############################
566    # Added by Artak
567    #################################3
568    
569    #Apply a sequence of k Givens rotations, used within gmres codes
570    # vrot=givapp(c, s, vin, k)
571    def givapp(c,s,vin):
572        vrot=vin # warning: vin is altered!!!!
573        if isinstance(c,float):
574            vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
575        else:
576            for i in range(len(c)):
577                w1=c[i]*vrot[i]-s[i]*vrot[i+1]
578            w2=s[i]*vrot[i]+c[i]*vrot[i+1]
579                vrot[i:i+2]=w1,w2
580        return vrot
581    
582    def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=10):
583       m=iter_restart
584       iter=0
585       while True:
586          if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
587          x,stopped=GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=x, iter_max=iter_max-iter, iter_restart=m)
588          iter+=iter_restart
589          if stopped: break
590       return x
591    
592    def GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=10):
593       iter=0
594       r=Msolve(b)
595       r_dot_r = bilinearform(r, r)
596       if r_dot_r<0: raise NegativeNorm,"negative norm."
597       norm_b=math.sqrt(r_dot_r)
598    
599       if x==None:
600          x=0*b
601       else:
602          r=Msolve(b-Aprod(x))
603          r_dot_r = bilinearform(r, r)
604          if r_dot_r<0: raise NegativeNorm,"negative norm."
605      
606       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
607       c=numarray.zeros(iter_restart,numarray.Float64)
608       s=numarray.zeros(iter_restart,numarray.Float64)
609       g=numarray.zeros(iter_restart,numarray.Float64)
610       v=[]
611    
612       rho=math.sqrt(r_dot_r)
613       v.append(r/rho)
614       g[0]=rho
615    
616       while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
617    
618        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
619    
620        
621        p=Msolve(Aprod(v[iter]))
622    
623        v.append(p)
624    
625        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
626    
627    # Modified Gram-Schmidt
628        for j in range(iter+1):
629          h[j][iter]=bilinearform(v[j],v[iter+1])  
630          v[iter+1]+=(-1.)*h[j][iter]*v[j]
631          
632        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
633        v_norm2=h[iter+1][iter]
634    
635    
636    # Reorthogonalize if needed
637        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
638         for j in range(iter+1):
639            hr=bilinearform(v[j],v[iter+1])
640                h[j][iter]=h[j][iter]+hr #vhat
641                v[iter+1] +=(-1.)*hr*v[j]
642    
643         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
644         h[iter+1][iter]=v_norm2
645    
646    #   watch out for happy breakdown
647            if v_norm2 != 0:
648             v[iter+1]=v[iter+1]/h[iter+1][iter]
649    
650    #   Form and store the information for the new Givens rotation
651        if iter > 0 :
652            hhat=[]
653            for i in range(iter+1) : hhat.append(h[i][iter])
654            hhat=givapp(c[0:iter],s[0:iter],hhat);
655                for i in range(iter+1) : h[i][iter]=hhat[i]
656    
657        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
658        if mu!=0 :
659            c[iter]=h[iter][iter]/mu
660            s[iter]=-h[iter+1][iter]/mu
661            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
662            h[iter+1][iter]=0.0
663            g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
664    
665    # Update the residual norm
666            rho=abs(g[iter+1])
667        iter+=1
668    
669    # At this point either iter > iter_max or rho < tol.
670    # It's time to compute x and leave.        
671    
672       if iter > 0 :
673         y=numarray.zeros(iter,numarray.Float64)    
674         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
675         if iter > 1 :  
676            i=iter-2  
677            while i>=0 :
678              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
679              i=i-1
680         xhat=v[iter-1]*y[iter-1]
681         for i in range(iter-1):
682        xhat += v[i]*y[i]
683       else : xhat=v[0]
684        
685       x += xhat
686       if iter!=iter_restart-1:
687          stopped=True
688       else:
689          stopped=False
690    
691       return x,stopped
692    
693    def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
694    
695        #
696        #  minres solves the system of linear equations Ax = b
697        #  where A is a symmetric matrix (possibly indefinite or singular)
698        #  and b is a given vector.
699        #  
700        #  "A" may be a dense or sparse matrix (preferably sparse!)
701        #  or the name of a function such that
702        #               y = A(x)
703        #  returns the product y = Ax for any given vector x.
704        #
705        #  "M" defines a positive-definite preconditioner M = C C'.
706        #  "M" may be a dense or sparse matrix (preferably sparse!)
707        #  or the name of a function such that
708        #  solves the system My = x for any given vector x.
709        #
710        #
711        
712        #------------------------------------------------------------------
713        # Set up y and v for the first Lanczos vector v1.
714        # y  =  beta1 P' v1,  where  P = C**(-1).
715        # v is really P' v1.
716        #------------------------------------------------------------------
717        if x==None:
718          x=0*b
719        else:
720          b += (-1)*Aprod(x)
721    
722        r1    = b
723        y = Msolve(b)
724        beta1 = bilinearform(b,y)
725    
726        if beta1< 0: raise NegativeNorm,"negative norm."
727    
728        #  If b = 0 exactly, stop with x = 0.
729        if beta1==0: return x*0.
730    
731        if beta1> 0:
732          beta1  = math.sqrt(beta1)       # Normalize y to get v1 later.
733    
734        #------------------------------------------------------------------
735        # Initialize other quantities.
736        # ------------------------------------------------------------------
737    #  Initialize                              
738    
739        iter   = 0
740        Anorm = 0
741        ynorm = 0
742    #    x=x*0
743    
744        oldb   = 0
745        beta   = beta1
746        dbar   = 0
747        epsln  = 0
748        phibar = beta1
749        rhs1   = beta1
750        rhs2   = 0
751        rnorm  = phibar
752        tnorm2 = 0
753        ynorm2 = 0
754        cs     = -1
755        sn     = 0
756        w      = b*0.
757        w2     = b*0.
758        r2     = r1
759        eps    = 0.0001
760    
761        #---------------------------------------------------------------------
762        # Main iteration loop.
763        # --------------------------------------------------------------------
764        while not stoppingcriterium(rnorm,Anorm*ynorm):    #  ||r|| / (||A|| ||x||)
765    
766        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
767            iter    = iter  +  1
768    
769            #-----------------------------------------------------------------
770            # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
771            # The general iteration is similar to the case k = 1 with v0 = 0:
772            #
773            #   p1      = Operator * v1  -  beta1 * v0,
774            #   alpha1  = v1'p1,
775            #   q2      = p2  -  alpha1 * v1,
776            #   beta2^2 = q2'q2,
777            #   v2      = (1/beta2) q2.
778            #
779            # Again, y = betak P vk,  where  P = C**(-1).
780            #-----------------------------------------------------------------
781            s = 1/beta                 # Normalize previous vector (in y).
782            v = s*y                    # v = vk if P = I
783        
784            y      = Aprod(v)
785        
786            if iter >= 2:
787              y = y - (beta/oldb)*r1
788    
789            alfa   = bilinearform(v,y)              # alphak
790            y      = (- alfa/beta)*r2 + y
791            r1     = r2
792            r2     = y
793            y = Msolve(r2)
794            oldb   = beta                           # oldb = betak
795            beta   = bilinearform(r2,y)             # beta = betak+1^2
796            if beta < 0: raise NegativeNorm,"negative norm."
797    
798            beta   = math.sqrt( beta )
799            tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
800            
801            if iter==1:                 # Initialize a few things.
802              gmax   = abs( alfa )      # alpha1
803              gmin   = gmax             # alpha1
804    
805            # Apply previous rotation Qk-1 to get
806            #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
807            #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
808        
809            oldeps = epsln
810            delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
811            gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
812            epsln  =               sn * beta  # epsln2 = 0         epslnk+1
813            dbar   =            -  cs * beta  # dbar 2 = beta2     dbar k+1
814    
815            # Compute the next plane rotation Qk
816    
817            gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
818            gamma  = max(gamma,eps)
819            cs     = gbar / gamma             # ck
820            sn     = beta / gamma             # sk
821            phi    = cs * phibar              # phik
822            phibar = sn * phibar              # phibark+1
823    
824            # Update  x.
825    
826            denom = 1/gamma
827            w1    = w2
828            w2    = w
829            w     = (v - oldeps*w1 - delta*w2) * denom
830            x     = x  +  phi*w
831    
832            # Go round again.
833    
834            gmax   = max(gmax,gamma)
835            gmin   = min(gmin,gamma)
836            z      = rhs1 / gamma
837            ynorm2 = z*z  +  ynorm2
838            rhs1   = rhs2 -  delta*z
839            rhs2   =      -  epsln*z
840    
841            # Estimate various norms and test for convergence.
842    
843            Anorm  = math.sqrt( tnorm2 )
844            ynorm  = math.sqrt( ynorm2 )
845    
846            rnorm  = phibar
847    
848        # Return final answer.
849        print iter
850        return x
851        
852    #############################################
853    
854  class ArithmeticTuple(object):  class ArithmeticTuple(object):
855     """     """
856     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 920  class ArithmeticTuple(object):
920             out.append(other*self[i])             out.append(other*self[i])
921         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
922    
923    #########################
924    # Added by Artak
925    #########################
926       def __div__(self,other):
927           """
928           dividing from the right
929    
930           @param other: scaling factor
931           @type other: C{float}
932           @return: itemwise self/other
933           @rtype: L{ArithmeticTuple}
934           """
935           out=[]
936           for i in range(len(self)):
937               out.append(self[i]/other)
938           return ArithmeticTuple(*tuple(out))
939    
940       def __rdiv__(self,other):
941           """
942           dividing from the left
943    
944           @param other: scaling factor
945           @type other: C{float}
946           @return: itemwise other/self
947           @rtype: L{ArithmeticTuple}
948           """
949           out=[]
950           for i in range(len(self)):
951               out.append(other/self[i])
952           return ArithmeticTuple(*tuple(out))
953      
954    ##########################################33
955    
956     def __iadd__(self,other):     def __iadd__(self,other):
957         """         """
958         in-place add of other to self         in-place add of other to self
# Line 621  class ArithmeticTuple(object): Line 966  class ArithmeticTuple(object):
966             self.__items[i]+=other[i]             self.__items[i]+=other[i]
967         return self         return self
968    
969    class HomogeneousSaddlePointProblem(object):
970          """
971          This provides a framwork for solving homogeneous saddle point problem of the form
972    
973                 Av+B^*p=f
974                 Bv    =0
975    
976          for the unknowns v and p and given operators A and B and given right hand side f.
977          B^* is the adjoint operator of B is the given inner product.
978    
979          """
980          def __init__(self,**kwargs):
981            self.setTolerance()
982            self.setToleranceReductionFactor()
983    
984          def initialize(self):
985            """
986            initialize the problem (overwrite)
987            """
988            pass
989          def B(self,v):
990             """
991             returns Bv (overwrite)
992             @rtype: equal to the type of p
993    
994             @note: boundary conditions on p should be zero!
995             """
996             pass
997    
998          def inner(self,p0,p1):
999             """
1000             returns inner product of two element p0 and p1  (overwrite)
1001            
1002             @type p0: equal to the type of p
1003             @type p1: equal to the type of p
1004             @rtype: C{float}
1005    
1006             @rtype: equal to the type of p
1007             """
1008             pass
1009    
1010          def solve_A(self,u,p):
1011             """
1012             solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
1013    
1014             @rtype: equal to the type of v
1015             @note: boundary conditions on v should be zero!
1016             """
1017             pass
1018    
1019          def solve_prec(self,p):
1020             """
1021             provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
1022    
1023             @rtype: equal to the type of p
1024             """
1025             pass
1026    
1027          def stoppingcriterium(self,Bv,v,p):
1028             """
1029             returns a True if iteration is terminated. (overwrite)
1030    
1031             @rtype: C{bool}
1032             """
1033             pass
1034                
1035          def __inner(self,p,r):
1036             return self.inner(p,r[1])
1037    
1038          def __inner_p(self,p1,p2):
1039             return self.inner(p1,p2)
1040    
1041          def __stoppingcriterium(self,norm_r,r,p):
1042              return self.stoppingcriterium(r[1],r[0],p)
1043    
1044          def __stoppingcriterium_GMRES(self,norm_r,norm_b):
1045              return self.stoppingcriterium_GMRES(norm_r,norm_b)
1046    
1047          def __stoppingcriterium_MINRES(self,norm_r,norm_Ax):
1048              return self.stoppingcriterium_MINRES(norm_r,norm_Ax)
1049    
1050    
1051          def setTolerance(self,tolerance=1.e-8):
1052                  self.__tol=tolerance
1053          def getTolerance(self):
1054                  return self.__tol
1055          def setToleranceReductionFactor(self,reduction=0.01):
1056                  self.__reduction=reduction
1057          def getSubProblemTolerance(self):
1058                  return self.__reduction*self.getTolerance()
1059    
1060          def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG'):
1061                  """
1062                  solves the saddle point problem using initial guesses v and p.
1063    
1064                  @param max_iter: maximum number of iteration steps.
1065                  """
1066                  self.verbose=verbose
1067                  self.show_details=show_details and self.verbose
1068    
1069                  # assume p is known: then v=A^-1(f-B^*p)
1070                  # which leads to BA^-1B^*p = BA^-1f  
1071    
1072              # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)      
1073    
1074              
1075              self.__z=v+self.solve_A(v,p*0)
1076    
1077                  Bz=self.B(self.__z)
1078                  #
1079              #   solve BA^-1B^*p = Bz
1080                  #
1081                  #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1082                  #
1083                  #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1084                  #                           A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)
1085                  #
1086                  self.iter=0
1087              if solver=='GMRES':      
1088                    if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1089                    p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1.)
1090                    # solve Au=f-B^*p
1091                    #       A(u-v)=f-B^*p-Av
1092                    #       u=v+(u-v)
1093            u=v+self.solve_A(v,p)
1094    
1095              if solver=='MINRES':      
1096                    if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1097                    p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_MINRES,iter_max=max_iter, x=p*1.)
1098                    # solve Au=f-B^*p
1099                    #       A(u-v)=f-B^*p-Av
1100                    #       u=v+(u-v)
1101            u=v+self.solve_A(v,p)
1102        
1103                  if solver=='PCG':
1104                    if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1105                    p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1106                u=r[0]  
1107    
1108                  print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1109    
1110              return u,p
1111    
1112          def __Msolve(self,r):
1113              return self.solve_prec(r[1])
1114    
1115          def __Msolve_GMRES(self,r):
1116              return self.solve_prec(r)
1117    
1118    
1119          def __Aprod(self,p):
1120              # return BA^-1B*p
1121              #solve Av =-B^*p as Av =f-Az-B^*p
1122              v=self.solve_A(self.__z,-p)
1123              return ArithmeticTuple(v, self.B(v))
1124    
1125          def __Aprod_GMRES(self,p):
1126              # return BA^-1B*p
1127              #solve Av =-B^*p as Av =f-Az-B^*p
1128          v=self.solve_A(self.__z,-p)
1129              return self.B(v)
1130    
1131  class SaddlePointProblem(object):  class SaddlePointProblem(object):
1132     """     """
1133     This implements a solver for a saddlepoint problem     This implements a solver for a saddlepoint problem
# Line 875  def MaskFromBoundaryTag(function_space,* Line 1382  def MaskFromBoundaryTag(function_space,*
1382     else:     else:
1383        return util.whereNonZero(util.interpolate(out,function_space))        return util.whereNonZero(util.interpolate(out,function_space))
1384    
1385    
1386    

Legend:
Removed from v.1331  
changed lines
  Added in v.1482

  ViewVC Help
Powered by ViewVC 1.1.26