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

revision 1331 by gross, Tue Oct 23 00:42:15 2007 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
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    ############################
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    #########################
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
770         """         """
771         in-place add of other to self         in-place add of other to self
# Line 621  class ArithmeticTuple(object): Line 779  class ArithmeticTuple(object):
779             self.__items[i]+=other[i]             self.__items[i]+=other[i]
780         return self         return self
781
783          """
784          This provides a framwork for solving homogeneous saddle point problem of the form
785
786                 Av+B^*p=f
787                 Bv    =0
788
789          for the unknowns v and p and given operators A and B and given right hand side f.
790          B^* is the adjoint operator of B is the given inner product.
791
792          """
793          def __init__(self,**kwargs):
794            self.setTolerance()
795            self.setToleranceReductionFactor()
796
797          def initialize(self):
798            """
799            initialize the problem (overwrite)
800            """
801            pass
802          def B(self,v):
803             """
804             returns Bv (overwrite)
805             @rtype: equal to the type of p
806
807             @note: boundary conditions on p should be zero!
808             """
809             pass
810
811          def inner(self,p0,p1):
812             """
813             returns inner product of two element p0 and p1  (overwrite)
814
815             @type p0: equal to the type of p
816             @type p1: equal to the type of p
817             @rtype: C{float}
818
819             @rtype: equal to the type of p
820             """
821             pass
822
823          def solve_A(self,u,p):
824             """
825             solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
826
827             @rtype: equal to the type of v
828             @note: boundary conditions on v should be zero!
829             """
830             pass
831
832          def solve_prec(self,p):
833             """
834             provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
835
836             @rtype: equal to the type of p
837             """
838             pass
839
840          def stoppingcriterium(self,Bv,v,p):
841             """
842             returns a True if iteration is terminated. (overwrite)
843
844             @rtype: C{bool}
845             """
846             pass
847
848          def __inner(self,p,r):
849             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):
855              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):
861                  self.__tol=tolerance
862          def getTolerance(self):
863                  return self.__tol
864          def setToleranceReductionFactor(self,reduction=0.01):
865                  self.__reduction=reduction
866          def getSubProblemTolerance(self):
867                  return self.__reduction*self.getTolerance()
868
869          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.
872
873                  @param max_iter: maximum number of iteration steps.
874                  """
875                  self.verbose=verbose
876                  self.show_details=show_details and self.verbose
877
878              # 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)
882
883                  Bz=self.B(self.__z)
884                  #
885              #   solve BA^-1B^*p = Bz
886                  #
887                  #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
888                  #
889                  #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
890                  #                           A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)
891                  #
892                  self.iter=0
893              if solver=='GMRES':
894                    if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
895                    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):
907              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):
914              # return BA^-1B*p
915              #solve Av =-B^*p as Av =f-Az-B^*p
916              v=self.solve_A(self.__z,p)
917              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