/[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 1464 by gross, Thu Feb 14 10:01:43 2008 UTC revision 1465 by artak, Wed Apr 2 03:28:25 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,Float32,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 539  type like argument C{x}. Line 544  type like argument C{x}.
544    
545     return x,r     return x,r
546    
547    
548    ############################
549    # Added by Artak
550    #################################3
551    
552    #Apply a sequence of k Givens rotations, used within gmres codes
553    # vrot=givapp(c, s, vin, k)
554    def givapp(c,s,vin,k):
555        vrot=vin
556        for i in range(k+1):
557            w1=c[i]*vrot[i]-s[i]*vrot[i+1]
558        w2=s[i]*vrot[i]+c[i]*vrot[i+1]
559            vrot[i:i+2]=w1,w2
560        return vrot
561    
562    def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
563    
564       from numarray import dot
565      
566       v0=b
567       iter=0
568       if x==None:
569          x=0*b
570       else:
571          b += (-1.)*Aprod(x)
572       r=b
573    
574       rhat=Msolve(r)
575      
576       rhat_dot_r = bilinearform(rhat, rhat)
577       norm_r=math.sqrt(rhat_dot_r)
578    
579       if rhat_dot_r<0: raise NegativeNorm,"negative norm."
580      
581       h=zeros((iter_max,iter_max),Float32)
582       c=zeros(iter_max,Float32)
583       s=zeros(iter_max,Float32)
584       g=zeros(iter_max,Float32)
585       v=[]
586    
587       v.append(rhat/norm_r)
588       rho=norm_r
589       g[0]=rho
590    
591       while not stoppingcriterium(rho,rho,norm_r):
592    
593        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
594    
595        
596        vhat=Aprod(v[iter])
597        p=Msolve(vhat)  
598    
599        v.append(p)
600    
601        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
602    
603    # Modified Gram-Schmidt
604        for j in range(iter+1):
605          h[j][iter]=bilinearform(v[j],v[iter+1])  
606          v[iter+1]+=(-1.)*h[j][iter]*v[j]
607          
608        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
609        v_norm2=h[iter+1][iter]
610    
611    
612    # Reorthogonalize if needed
613        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
614         for j in range(iter+1):
615            hr=bilinearform(v[j],v[iter+1])
616                h[j][iter]=h[j][iter]+hr #vhat
617                v[iter+1] +=(-1.)*hr*v[j]
618    
619         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
620         h[iter+1][iter]=v_norm2
621    
622    #   watch out for happy breakdown
623            if v_norm2 != 0:
624             v[iter+1]=v[iter+1]/h[iter+1][iter]
625    
626    #   Form and store the information for the new Givens rotation
627        if iter > 0 :
628            hhat=[]
629            for i in range(iter+1) : hhat.append(h[i][iter])
630            hhat=givapp(c[0:iter],s[0:iter],hhat,iter-1);
631                for i in range(iter+1) : h[i][iter]=hhat[i]
632    
633        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
634        if mu!=0 :
635            c[iter]=h[iter][iter]/mu
636            s[iter]=-h[iter+1][iter]/mu
637            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
638            h[iter+1][iter]=0.0
639            g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2],0)
640    
641    # Update the residual norm
642            rho=abs(g[iter+1])
643        iter+=1
644    
645    # At this point either iter > iter_max or rho < tol.
646    # It's time to compute x and leave.        
647    
648       if iter > 0 :
649         y=zeros(iter,Float32)  
650         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
651         if iter > 1 :  
652            i=iter-2  
653            while i>=0 :
654              y[i] = ( g[i] - dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
655              i=i-1
656         xhat=v[iter-1]*y[iter-1]
657         for i in range(iter-1):
658        xhat += v[i]*y[i]
659       else : xhat=v[0]
660        
661       x += xhat
662    
663       return x
664        
665    #############################################
666    
667  class ArithmeticTuple(object):  class ArithmeticTuple(object):
668     """     """
669     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 733  class ArithmeticTuple(object):
733             out.append(other*self[i])             out.append(other*self[i])
734         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
735    
736    #########################
737    # Added by Artak
738    #########################
739       def __div__(self,other):
740           """
741           dividing from the right
742    
743           @param other: scaling factor
744           @type other: C{float}
745           @return: itemwise self/other
746           @rtype: L{ArithmeticTuple}
747           """
748           out=[]
749           for i in range(len(self)):
750               out.append(self[i]/other)
751           return ArithmeticTuple(*tuple(out))
752    
753       def __rdiv__(self,other):
754           """
755           dividing from the left
756    
757           @param other: scaling factor
758           @type other: C{float}
759           @return: itemwise other/self
760           @rtype: L{ArithmeticTuple}
761           """
762           out=[]
763           for i in range(len(self)):
764               out.append(other/self[i])
765           return ArithmeticTuple(*tuple(out))
766      
767    ##########################################33
768    
769     def __iadd__(self,other):     def __iadd__(self,other):
770         """         """
771         in-place add of other to self         in-place add of other to self
# Line 690  class HomogeneousSaddlePointProblem(obje Line 848  class HomogeneousSaddlePointProblem(obje
848        def __inner(self,p,r):        def __inner(self,p,r):
849           return self.inner(p,r[1])           return self.inner(p,r[1])
850    
851          def __inner_p(self,p1,p2):
852             return self.inner(p1,p2)
853    
854        def __stoppingcriterium(self,norm_r,r,p):        def __stoppingcriterium(self,norm_r,r,p):
855            return self.stoppingcriterium(r[1],r[0],p)            return self.stoppingcriterium(r[1],r[0],p)
856    
857          def __stoppingcriterium_GMRES(self,norm_r,rho,r):
858              return self.stoppingcriterium_GMRES(rho,r)
859    
860        def setTolerance(self,tolerance=1.e-8):        def setTolerance(self,tolerance=1.e-8):
861                self.__tol=tolerance                self.__tol=tolerance
862        def getTolerance(self):        def getTolerance(self):
# Line 702  class HomogeneousSaddlePointProblem(obje Line 866  class HomogeneousSaddlePointProblem(obje
866        def getSubProblemTolerance(self):        def getSubProblemTolerance(self):
867                return self.__reduction*self.getTolerance()                return self.__reduction*self.getTolerance()
868    
869        def solve(self,v,p,max_iter=20, verbose=False, show_details=False):        def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='GMRES'):
870                """                """
871                solves the saddle point problem using initial guesses v and p.                solves the saddle point problem using initial guesses v and p.
872    
# Line 713  class HomogeneousSaddlePointProblem(obje Line 877  class HomogeneousSaddlePointProblem(obje
877    
878            # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)                  # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)      
879    
880              
881            self.__z=v+self.solve_A(v,p*0)            self.__z=v+self.solve_A(v,p*0)
882    
883                Bz=self.B(self.__z)                Bz=self.B(self.__z)
884                #                #
885            #   solve BA^-1B^*p = Bz            #   solve BA^-1B^*p = Bz
# Line 724  class HomogeneousSaddlePointProblem(obje Line 890  class HomogeneousSaddlePointProblem(obje
890                #                           A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)                #                           A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)
891                #                #
892                self.iter=0                self.iter=0
893                if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter            if solver=='GMRES':      
894                p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)                  if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
895            return r[0],p                  p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1)
896            u=v+self.solve_A(v,p)
897        
898                  else:
899                    if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
900                    p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p*1)
901                u=r[0]  
902                    u=v+self.solve_A(v,p)     # Lutz to check !!!!!
903    
904              return u,p
905    
906        def __Msolve(self,r):        def __Msolve(self,r):
907            return self.solve_prec(r[1])            return self.solve_prec(r[1])
908    
909          def __Msolve_GMRES(self,r):
910              return self.solve_prec(r)
911    
912    
913        def __Aprod(self,p):        def __Aprod(self,p):
914            # return BA^-1B*p            # return BA^-1B*p
915            #solve Av =-B^*p as Av =f-Az-B^*p            #solve Av =-B^*p as Av =f-Az-B^*p
916            v=self.solve_A(self.__z,p)            v=self.solve_A(self.__z,p)
917            return ArithmeticTuple(v, self.B(v))            return ArithmeticTuple(v, self.B(v))
918    
919          def __Aprod_GMRES(self,p):
920              # return BA^-1B*p
921              #solve Av =-B^*p as Av =f-Az-B^*p
922          v=self.solve_A(self.__z,p)
923              return self.B(v)
924    
925  class SaddlePointProblem(object):  class SaddlePointProblem(object):
926     """     """
# Line 993  def MaskFromBoundaryTag(function_space,* Line 1177  def MaskFromBoundaryTag(function_space,*
1177        return util.whereNonZero(util.interpolate(out,function_space))        return util.whereNonZero(util.interpolate(out,function_space))
1178    
1179    
1180    
Legend:
Removed from v.1464  
changed lines
  Added in v.1465
  ViewVC Help
Powered by ViewVC 1.1.26