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

revision 1331 by gross, Tue Oct 23 00:42:15 2007 UTC revision 1476 by gross, Mon Apr 7 23:38:50 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 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    ############################
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    #############################################
694
695  class ArithmeticTuple(object):  class ArithmeticTuple(object):
696     """     """
697     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 761  class ArithmeticTuple(object):
761             out.append(other*self[i])             out.append(other*self[i])
762         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
763
764    #########################
766    #########################
767       def __div__(self,other):
768           """
769           dividing from the right
770
771           @param other: scaling factor
772           @type other: C{float}
773           @return: itemwise self/other
774           @rtype: L{ArithmeticTuple}
775           """
776           out=[]
777           for i in range(len(self)):
778               out.append(self[i]/other)
779           return ArithmeticTuple(*tuple(out))
780
781       def __rdiv__(self,other):
782           """
783           dividing from the left
784
785           @param other: scaling factor
786           @type other: C{float}
787           @return: itemwise other/self
788           @rtype: L{ArithmeticTuple}
789           """
790           out=[]
791           for i in range(len(self)):
792               out.append(other/self[i])
793           return ArithmeticTuple(*tuple(out))
794
795    ##########################################33
796
798         """         """
799         in-place add of other to self         in-place add of other to self
# Line 621  class ArithmeticTuple(object): Line 807  class ArithmeticTuple(object):
807             self.__items[i]+=other[i]             self.__items[i]+=other[i]
808         return self         return self
809
811          """
812          This provides a framwork for solving homogeneous saddle point problem of the form
813
814                 Av+B^*p=f
815                 Bv    =0
816
817          for the unknowns v and p and given operators A and B and given right hand side f.
818          B^* is the adjoint operator of B is the given inner product.
819
820          """
821          def __init__(self,**kwargs):
822            self.setTolerance()
823            self.setToleranceReductionFactor()
824
825          def initialize(self):
826            """
827            initialize the problem (overwrite)
828            """
829            pass
830          def B(self,v):
831             """
832             returns Bv (overwrite)
833             @rtype: equal to the type of p
834
835             @note: boundary conditions on p should be zero!
836             """
837             pass
838
839          def inner(self,p0,p1):
840             """
841             returns inner product of two element p0 and p1  (overwrite)
842
843             @type p0: equal to the type of p
844             @type p1: equal to the type of p
845             @rtype: C{float}
846
847             @rtype: equal to the type of p
848             """
849             pass
850
851          def solve_A(self,u,p):
852             """
853             solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
854
855             @rtype: equal to the type of v
856             @note: boundary conditions on v should be zero!
857             """
858             pass
859
860          def solve_prec(self,p):
861             """
862             provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
863
864             @rtype: equal to the type of p
865             """
866             pass
867
868          def stoppingcriterium(self,Bv,v,p):
869             """
870             returns a True if iteration is terminated. (overwrite)
871
872             @rtype: C{bool}
873             """
874             pass
875
876          def __inner(self,p,r):
877             return self.inner(p,r[1])
878
879          def __inner_p(self,p1,p2):
880             return self.inner(p1,p2)
881
882          def __stoppingcriterium(self,norm_r,r,p):
883              return self.stoppingcriterium(r[1],r[0],p)
884
885          def __stoppingcriterium_GMRES(self,norm_r,norm_b):
886              return self.stoppingcriterium_GMRES(norm_r,norm_b)
887
888          def setTolerance(self,tolerance=1.e-8):
889                  self.__tol=tolerance
890          def getTolerance(self):
891                  return self.__tol
892          def setToleranceReductionFactor(self,reduction=0.01):
893                  self.__reduction=reduction
894          def getSubProblemTolerance(self):
895                  return self.__reduction*self.getTolerance()
896
897          def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG'):
898                  """
899                  solves the saddle point problem using initial guesses v and p.
900
901                  @param max_iter: maximum number of iteration steps.
902                  """
903                  self.verbose=verbose
904                  self.show_details=show_details and self.verbose
905
906                  # assume p is known: then v=A^-1(f-B^*p)
907                  # which leads to BA^-1B^*p = BA^-1f
908
909              # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)
910
911
912              self.__z=v+self.solve_A(v,p*0)
913
914                  Bz=self.B(self.__z)
915                  #
916              #   solve BA^-1B^*p = Bz
917                  #
918                  #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
919                  #
920                  #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
921                  #                           A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)
922                  #
923                  self.iter=0
924              if solver=='GMRES':
925                    if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
926                    p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1.)
927                    # solve Au=f-B^*p
928                    #       A(u-v)=f-B^*p-Av
929                    #       u=v+(u-v)
930            u=v+self.solve_A(v,p)
931
932                  else:
933                    if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
934                    p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
935                u=r[0]
936                  print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
937
938              return u,p
939
940          def __Msolve(self,r):
941              return self.solve_prec(r[1])
942
943          def __Msolve_GMRES(self,r):
944              return self.solve_prec(r)
945
946
947          def __Aprod(self,p):
948              # return BA^-1B*p
949              #solve Av =-B^*p as Av =f-Az-B^*p
950              v=self.solve_A(self.__z,-p)
951              return ArithmeticTuple(v, self.B(v))
952
953          def __Aprod_GMRES(self,p):
954              # return BA^-1B*p
955              #solve Av =-B^*p as Av =f-Az-B^*p
956          v=self.solve_A(self.__z,-p)
957              return self.B(v)
958