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

revision 1331 by gross, Tue Oct 23 00:42:15 2007 UTC revision 1467 by gross, Wed Apr 2 08:10:37 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 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):
555        vrot=vin # warning: vin is altered!!!!
556        if isinstance(c,float):
557            vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
558        else:
559            for i in range(len(c)):
560                w1=c[i]*vrot[i]-s[i]*vrot[i+1]
561            w2=s[i]*vrot[i]+c[i]*vrot[i+1]
562                vrot[i:i+2]=w1,w2
563        return vrot
564
565    def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
566       iter=0
567       r=Msolve(b)
568       r_dot_r = bilinearform(r, r)
569       if r_dot_r<0: raise NegativeNorm,"negative norm."
570       norm_b=math.sqrt(r_dot_r)
571
572       if x==None:
573          x=0*b
574       else:
575          r=Msolve(b-Aprod(x))
576          r_dot_r = bilinearform(r, r)
577          if r_dot_r<0: raise NegativeNorm,"negative norm."
578
579       h=numarray.zeros((iter_max,iter_max),numarray.Float64)
580       c=numarray.zeros(iter_max,numarray.Float64)
581       s=numarray.zeros(iter_max,numarray.Float64)
582       g=numarray.zeros(iter_max,numarray.Float64)
583       v=[]
584
585       rho=math.sqrt(r_dot_r)
586       v.append(r/rho)
587       g[0]=rho
588
589       while not stoppingcriterium(rho,norm_b):
590
591        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
592
593
594        p=Msolve(Aprod(v[iter]))
595
596        v.append(p)
597
598        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
599
600    # Modified Gram-Schmidt
601        for j in range(iter+1):
602          h[j][iter]=bilinearform(v[j],v[iter+1])
603          v[iter+1]+=(-1.)*h[j][iter]*v[j]
604
605        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
606        v_norm2=h[iter+1][iter]
607
608
609    # Reorthogonalize if needed
610        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
611         for j in range(iter+1):
612            hr=bilinearform(v[j],v[iter+1])
613                h[j][iter]=h[j][iter]+hr #vhat
614                v[iter+1] +=(-1.)*hr*v[j]
615
616         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
617         h[iter+1][iter]=v_norm2
618
619    #   watch out for happy breakdown
620            if v_norm2 != 0:
621             v[iter+1]=v[iter+1]/h[iter+1][iter]
622
623    #   Form and store the information for the new Givens rotation
624        if iter > 0 :
625            hhat=[]
626            for i in range(iter+1) : hhat.append(h[i][iter])
627            hhat=givapp(c[0:iter],s[0:iter],hhat);
628                for i in range(iter+1) : h[i][iter]=hhat[i]
629
630        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
631        if mu!=0 :
632            c[iter]=h[iter][iter]/mu
633            s[iter]=-h[iter+1][iter]/mu
634            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
635            h[iter+1][iter]=0.0
636            g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
637
638    # Update the residual norm
639            rho=abs(g[iter+1])
640        iter+=1
641
642    # At this point either iter > iter_max or rho < tol.
643    # It's time to compute x and leave.
644
645       if iter > 0 :
646         y=numarray.zeros(iter,numarray.Float64)
647         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
648         if iter > 1 :
649            i=iter-2
650            while i>=0 :
651              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
652              i=i-1
653         xhat=v[iter-1]*y[iter-1]
654         for i in range(iter-1):
655        xhat += v[i]*y[i]
656       else : xhat=v[0]
657
658       x += xhat
659
660       return x
661
662    #############################################
663
664  class ArithmeticTuple(object):  class ArithmeticTuple(object):
665     """     """
666     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 730  class ArithmeticTuple(object):
730             out.append(other*self[i])             out.append(other*self[i])
731         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
732
733    #########################
735    #########################
736       def __div__(self,other):
737           """
738           dividing from the right
739
740           @param other: scaling factor
741           @type other: C{float}
742           @return: itemwise self/other
743           @rtype: L{ArithmeticTuple}
744           """
745           out=[]
746           for i in range(len(self)):
747               out.append(self[i]/other)
748           return ArithmeticTuple(*tuple(out))
749
750       def __rdiv__(self,other):
751           """
752           dividing from the left
753
754           @param other: scaling factor
755           @type other: C{float}
756           @return: itemwise other/self
757           @rtype: L{ArithmeticTuple}
758           """
759           out=[]
760           for i in range(len(self)):
761               out.append(other/self[i])
762           return ArithmeticTuple(*tuple(out))
763
764    ##########################################33
765
767         """         """
768         in-place add of other to self         in-place add of other to self
# Line 621  class ArithmeticTuple(object): Line 776  class ArithmeticTuple(object):
776             self.__items[i]+=other[i]             self.__items[i]+=other[i]
777         return self         return self
778
780          """
781          This provides a framwork for solving homogeneous saddle point problem of the form
782
783                 Av+B^*p=f
784                 Bv    =0
785
786          for the unknowns v and p and given operators A and B and given right hand side f.
787          B^* is the adjoint operator of B is the given inner product.
788
789          """
790          def __init__(self,**kwargs):
791            self.setTolerance()
792            self.setToleranceReductionFactor()
793
794          def initialize(self):
795            """
796            initialize the problem (overwrite)
797            """
798            pass
799          def B(self,v):
800             """
801             returns Bv (overwrite)
802             @rtype: equal to the type of p
803
804             @note: boundary conditions on p should be zero!
805             """
806             pass
807
808          def inner(self,p0,p1):
809             """
810             returns inner product of two element p0 and p1  (overwrite)
811
812             @type p0: equal to the type of p
813             @type p1: equal to the type of p
814             @rtype: C{float}
815
816             @rtype: equal to the type of p
817             """
818             pass
819
820          def solve_A(self,u,p):
821             """
822             solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
823
824             @rtype: equal to the type of v
825             @note: boundary conditions on v should be zero!
826             """
827             pass
828
829          def solve_prec(self,p):
830             """
831             provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
832
833             @rtype: equal to the type of p
834             """
835             pass
836
837          def stoppingcriterium(self,Bv,v,p):
838             """
839             returns a True if iteration is terminated. (overwrite)
840
841             @rtype: C{bool}
842             """
843             pass
844
845          def __inner(self,p,r):
846             return self.inner(p,r[1])
847
848          def __inner_p(self,p1,p2):
849             return self.inner(p1,p2)
850
851          def __stoppingcriterium(self,norm_r,r,p):
852              return self.stoppingcriterium(r[1],r[0],p)
853
854          def __stoppingcriterium_GMRES(self,norm_r,norm_b):
855              return self.stoppingcriterium_GMRES(norm_r,norm_b)
856
857          def setTolerance(self,tolerance=1.e-8):
858                  self.__tol=tolerance
859          def getTolerance(self):
860                  return self.__tol
861          def setToleranceReductionFactor(self,reduction=0.01):
862                  self.__reduction=reduction
863          def getSubProblemTolerance(self):
864                  return self.__reduction*self.getTolerance()
865
866          def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='GMRES'):
867                  """
868                  solves the saddle point problem using initial guesses v and p.
869
870                  @param max_iter: maximum number of iteration steps.
871                  """
872                  self.verbose=verbose
873                  self.show_details=show_details and self.verbose
874
875              # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)
876
877
878              self.__z=v+self.solve_A(v,p*0)
879
880                  Bz=self.B(self.__z)
881                  #
882              #   solve BA^-1B^*p = Bz
883                  #
884                  #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
885                  #
886                  #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
887                  #                           A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)
888                  #
889                  self.iter=0
890              if solver=='GMRES':
891                    if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
892                    p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1.)
893                    # solve Au=f-B^*p
894                    #       A(u-v)=f-B^*p-Av
895                    #       u=v+(u-v)
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)
901                u=r[0]
902
903              return u,p
904
905          def __Msolve(self,r):
906              return self.solve_prec(r[1])
907
908          def __Msolve_GMRES(self,r):
909              return self.solve_prec(r)
910
911
912          def __Aprod(self,p):
913              # return BA^-1B*p
914              #solve Av =-B^*p as Av =f-Az-B^*p
915              v=self.solve_A(self.__z,p)
916              return ArithmeticTuple(v, -self.B(v))
917
918          def __Aprod_GMRES(self,p):
919              # return BA^-1B*p
920              #solve Av =-B^*p as Av =f-Az-B^*p
921          v=self.solve_A(self.__z,p)
922              return -self.B(v)
923