/[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 1312 by ksteube, Mon Sep 24 06:18:44 2007 UTC revision 1469 by gross, Thu Apr 3 05:16:56 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 435  class NegativeNorm(SolverSchemeException Line 440  class NegativeNorm(SolverSchemeException
440     """     """
441     pass     pass
442    
443  def PCG(b,x,Aprod,Msolve,bilinearform, norm, verbose=True, iter_max=100, tolerance=math.sqrt(util.EPSILON)):  class IterationHistory(object):
444     """     """
445     Solver for     The IterationHistory class is used to define a stopping criterium. It keeps track of the
446       residual norms. The stoppingcriterium indicates termination if the residual norm has been reduced by
447       a given tolerance.
448       """
449       def __init__(self,tolerance=math.sqrt(util.EPSILON),verbose=False):
450          """
451          Initialization
452    
453     M{Ax=b}        @param tolerance: tolerance
454          @type tolerance: positive C{float}
455          @param verbose: switches on the printing out some information
456          @type verbose: C{bool}
457          """
458          if not tolerance>0.:
459              raise ValueError,"tolerance needs to be positive."
460          self.tolerance=tolerance
461          self.verbose=verbose
462          self.history=[]
463       def stoppingcriterium(self,norm_r,r,x):
464           """
465           returns True if the C{norm_r} is C{tolerance}*C{norm_r[0]} where C{norm_r[0]}  is the residual norm at the first call.
466    
467     with a symmetric and positive definite operator A (more details required!).        
468     It uses the conjugate gradient method with preconditioner M providing an approximation of A.         @param norm_r: current residual norm
469           @type norm_r: non-negative C{float}
470           @param r: current residual (not used)
471           @param x: current solution approximation (not used)
472           @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
473           @rtype: C{bool}
474    
475           """
476           self.history.append(norm_r)
477           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]
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     The iteration is terminated if        
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     M{norm(r) <= tolerance * norm(b)}         """
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     where C{norm()} defines a norm and  def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
498       """
499       Solver for
500    
501     M{r = b-Ax}     M{Ax=b}
502    
503     the residual.     with a symmetric and positive definite operator A (more details required!).
504       It uses the conjugate gradient method with preconditioner M providing an approximation of A.
505    
506       The iteration is terminated if the C{stoppingcriterium} function return C{True}.
507    
508     For details on the preconditioned conjugate gradient method see the book:     For details on the preconditioned conjugate gradient method see the book:
509    
# Line 461  def PCG(b,x,Aprod,Msolve,bilinearform, n Line 512  def PCG(b,x,Aprod,Msolve,bilinearform, n
512     C. Romine, and H. van der Vorst.     C. Romine, and H. van der Vorst.
513    
514     @param b: the right hand side of the liner system. C{b} is altered.     @param b: the right hand side of the liner system. C{b} is altered.
515     @type b: any object type R supporting inplace add (x+=y) and scaling (x=scalar*y)     @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
    @param x: an initial guess for the solution  
    @type x: any object type S supporting inplace add (x+=y), scaling (x=scalar*y)  
516     @param Aprod: returns the value Ax     @param Aprod: returns the value Ax
517     @type Aprod: function C{Aprod(x)} where C{x} is of object type S. The returned object needs to be of type R.     @type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}. The returned object needs to be of the same type like argument C{b}.
518     @param Msolve: solves Mx=r     @param Msolve: solves Mx=r
519     @type Msolve: function C{Msolve(r)} where C{r} is of object type R. The returned object needs to be of tupe S.     @type Msolve: function C{Msolve(r)} where C{r} is of the same type like argument C{b}. The returned object needs to be of the same
520    type like argument C{x}.
521     @param bilinearform: inner product C{<x,r>}     @param bilinearform: inner product C{<x,r>}
522     @type bilinearform: function C{bilinearform(x,r)} where C{x} is of object type S and C{r} is of object type R. The returned value is a C{float}.     @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is . The returned value is a C{float}.
523     @param norm: norm calculation for the residual C{r=b-Ax}.     @param stoppingcriterium: function which returns True if a stopping criterium is meet. C{stoppingcriterium} has the arguments C{norm_r}, C{r} and C{x} giving the current norm of the residual (=C{sqrt(bilinearform(Msolve(r),r)}), the current residual and the current solution approximation. C{stoppingcriterium} is called in each iteration step.
524     @type norm: function C{norm(r)} where C{r} is of object type R. The returned value is a C{float}.     @type stoppingcriterium: function that returns C{True} or C{False}
525     @param verbose: switches on the printing out some information     @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
526     @type verbose: C{bool}     @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
527     @param iter_max: maximum number of iteration steps.     @param iter_max: maximum number of iteration steps.
528     @type iter_max: C{int}     @type iter_max: C{int}
529     @param tolerance: tolerance     @return: the solution approximation and the corresponding residual
530     @type tolerance: positive C{float}     @rtype: C{tuple}
531     @return: the solution apprximation and the corresponding residual     @warning: C{b} and C{x} are altered.
    @rtype: C{tuple} of an S type and and an R type object.A  
    @warning: C{b} ans C{x} are altered.  
532     """     """
    if verbose:  
         print "Enter PCG for solving Ax=b\n\t iter_max =%s\t tolerance   =%e"%(iter_max,tolerance)  
533     iter=0     iter=0
534     normb = norm(b)     if x==None:
535     if normb<0: raise NegativeNorm        x=0*b
536       else:
537     b += (-1)*Aprod(x)        b += (-1)*Aprod(x)
538     r=b     r=b
539     rhat=Msolve(r)     rhat=Msolve(r)
540     d = rhat;     d = rhat
541     rhat_dot_r = bilinearform(rhat, r)     rhat_dot_r = bilinearform(rhat, r)
542       if rhat_dot_r<0: raise NegativeNorm,"negative norm."
543    
544     while True:     while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
545         normr=norm(r)         iter+=1
        if normr<0: raise NegativeNorm  
        if verbose: print "iter: %s: norm(r) = %e, tolerance*norm(b) = %e"%(iter, normr,tolerance * normb)  
        if normr <= tolerance * normb: return x,r  
   
        iter+=1 # k = iter = 1 first time through  
546         if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max         if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
547    
548         q=Aprod(d)         q=Aprod(d)
# Line 515  def PCG(b,x,Aprod,Msolve,bilinearform, n Line 557  def PCG(b,x,Aprod,Msolve,bilinearform, n
557         d=rhat         d=rhat
558    
559         rhat_dot_r = rhat_dot_r_new         rhat_dot_r = rhat_dot_r_new
560           if rhat_dot_r<0: raise NegativeNorm,"negative norm."
561    
562       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):
583       iter=0
584       r=Msolve(b)
585       r_dot_r = bilinearform(r, r)
586       if r_dot_r<0: raise NegativeNorm,"negative norm."
587       norm_b=math.sqrt(r_dot_r)
588    
589       if x==None:
590          x=0*b
591       else:
592          r=Msolve(b-Aprod(x))
593          r_dot_r = bilinearform(r, r)
594          if r_dot_r<0: raise NegativeNorm,"negative norm."
595      
596       h=numarray.zeros((iter_max,iter_max),numarray.Float64)
597       c=numarray.zeros(iter_max,numarray.Float64)
598       s=numarray.zeros(iter_max,numarray.Float64)
599       g=numarray.zeros(iter_max,numarray.Float64)
600       v=[]
601    
602       rho=math.sqrt(r_dot_r)
603       v.append(r/rho)
604       g[0]=rho
605    
606       while not stoppingcriterium(rho,norm_b):
607    
608        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
609    
610        
611        p=Msolve(Aprod(v[iter]))
612    
613        v.append(p)
614    
615        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
616    
617    # Modified Gram-Schmidt
618        for j in range(iter+1):
619          h[j][iter]=bilinearform(v[j],v[iter+1])  
620          v[iter+1]+=(-1.)*h[j][iter]*v[j]
621          
622        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
623        v_norm2=h[iter+1][iter]
624    
625    
626    # Reorthogonalize if needed
627        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
628         for j in range(iter+1):
629            hr=bilinearform(v[j],v[iter+1])
630                h[j][iter]=h[j][iter]+hr #vhat
631                v[iter+1] +=(-1.)*hr*v[j]
632    
633         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
634         h[iter+1][iter]=v_norm2
635    
636    #   watch out for happy breakdown
637            if v_norm2 != 0:
638             v[iter+1]=v[iter+1]/h[iter+1][iter]
639    
640    #   Form and store the information for the new Givens rotation
641        if iter > 0 :
642            hhat=[]
643            for i in range(iter+1) : hhat.append(h[i][iter])
644            hhat=givapp(c[0:iter],s[0:iter],hhat);
645                for i in range(iter+1) : h[i][iter]=hhat[i]
646    
647        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
648        if mu!=0 :
649            c[iter]=h[iter][iter]/mu
650            s[iter]=-h[iter+1][iter]/mu
651            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
652            h[iter+1][iter]=0.0
653            g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
654    
655    # Update the residual norm
656            rho=abs(g[iter+1])
657        iter+=1
658    
659    # At this point either iter > iter_max or rho < tol.
660    # It's time to compute x and leave.        
661    
662       if iter > 0 :
663         y=numarray.zeros(iter,numarray.Float64)    
664         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
665         if iter > 1 :  
666            i=iter-2  
667            while i>=0 :
668              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
669              i=i-1
670         xhat=v[iter-1]*y[iter-1]
671         for i in range(iter-1):
672        xhat += v[i]*y[i]
673       else : xhat=v[0]
674        
675       x += xhat
676    
677       return x
678        
679    #############################################
680    
681    class ArithmeticTuple(object):
682       """
683       tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
684    
685       example of usage:
686    
687       from esys.escript import Data
688       from numarray import array
689       a=Data(...)
690       b=array([1.,4.])
691       x=ArithmeticTuple(a,b)
692       y=5.*x
693    
694       """
695       def __init__(self,*args):
696           """
697           initialize object with elements args.
698    
699           @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
700           """
701           self.__items=list(args)
702    
703       def __len__(self):
704           """
705           number of items
706    
707           @return: number of items
708           @rtype: C{int}
709           """
710           return len(self.__items)
711    
712       def __getitem__(self,index):
713           """
714           get an item
715    
716           @param index: item to be returned
717           @type index: C{int}
718           @return: item with index C{index}
719           """
720           return self.__items.__getitem__(index)
721    
722       def __mul__(self,other):
723           """
724           scaling from the right
725    
726           @param other: scaling factor
727           @type other: C{float}
728           @return: itemwise self*other
729           @rtype: L{ArithmeticTuple}
730           """
731           out=[]
732           for i in range(len(self)):
733               out.append(self[i]*other)
734           return ArithmeticTuple(*tuple(out))
735    
736       def __rmul__(self,other):
737           """
738           scaling from the left
739    
740           @param other: scaling factor
741           @type other: C{float}
742           @return: itemwise other*self
743           @rtype: L{ArithmeticTuple}
744           """
745           out=[]
746           for i in range(len(self)):
747               out.append(other*self[i])
748           return ArithmeticTuple(*tuple(out))
749    
750    #########################
751    # Added by Artak
752    #########################
753       def __div__(self,other):
754           """
755           dividing from the right
756    
757           @param other: scaling factor
758           @type other: C{float}
759           @return: itemwise self/other
760           @rtype: L{ArithmeticTuple}
761           """
762           out=[]
763           for i in range(len(self)):
764               out.append(self[i]/other)
765           return ArithmeticTuple(*tuple(out))
766    
767       def __rdiv__(self,other):
768           """
769           dividing from the left
770    
771           @param other: scaling factor
772           @type other: C{float}
773           @return: itemwise other/self
774           @rtype: L{ArithmeticTuple}
775           """
776           out=[]
777           for i in range(len(self)):
778               out.append(other/self[i])
779           return ArithmeticTuple(*tuple(out))
780      
781    ##########################################33
782    
783       def __iadd__(self,other):
784           """
785           in-place add of other to self
786    
787           @param other: increment
788           @type other: C{ArithmeticTuple}
789           """
790           if len(self) != len(other):
791               raise ValueError,"tuple length must match."
792           for i in range(len(self)):
793               self.__items[i]+=other[i]
794           return self
795    
796    class HomogeneousSaddlePointProblem(object):
797          """
798          This provides a framwork for solving homogeneous saddle point problem of the form
799    
800                 Av+B^*p=f
801                 Bv    =0
802    
803          for the unknowns v and p and given operators A and B and given right hand side f.
804          B^* is the adjoint operator of B is the given inner product.
805    
806          """
807          def __init__(self,**kwargs):
808            self.setTolerance()
809            self.setToleranceReductionFactor()
810    
811          def initialize(self):
812            """
813            initialize the problem (overwrite)
814            """
815            pass
816          def B(self,v):
817             """
818             returns Bv (overwrite)
819             @rtype: equal to the type of p
820    
821             @note: boundary conditions on p should be zero!
822             """
823             pass
824    
825          def inner(self,p0,p1):
826             """
827             returns inner product of two element p0 and p1  (overwrite)
828            
829             @type p0: equal to the type of p
830             @type p1: equal to the type of p
831             @rtype: C{float}
832    
833             @rtype: equal to the type of p
834             """
835             pass
836    
837          def solve_A(self,u,p):
838             """
839             solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
840    
841             @rtype: equal to the type of v
842             @note: boundary conditions on v should be zero!
843             """
844             pass
845    
846          def solve_prec(self,p):
847             """
848             provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
849    
850             @rtype: equal to the type of p
851             """
852             pass
853    
854          def stoppingcriterium(self,Bv,v,p):
855             """
856             returns a True if iteration is terminated. (overwrite)
857    
858             @rtype: C{bool}
859             """
860             pass
861                
862          def __inner(self,p,r):
863             return self.inner(p,r[1])
864    
865          def __inner_p(self,p1,p2):
866             return self.inner(p1,p2)
867    
868          def __stoppingcriterium(self,norm_r,r,p):
869              return self.stoppingcriterium(r[1],r[0],p)
870    
871          def __stoppingcriterium_GMRES(self,norm_r,norm_b):
872              return self.stoppingcriterium_GMRES(norm_r,norm_b)
873    
874          def setTolerance(self,tolerance=1.e-8):
875                  self.__tol=tolerance
876          def getTolerance(self):
877                  return self.__tol
878          def setToleranceReductionFactor(self,reduction=0.01):
879                  self.__reduction=reduction
880          def getSubProblemTolerance(self):
881                  return self.__reduction*self.getTolerance()
882    
883          def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='GMRES'):
884                  """
885                  solves the saddle point problem using initial guesses v and p.
886    
887                  @param max_iter: maximum number of iteration steps.
888                  """
889                  self.verbose=verbose
890                  self.show_details=show_details and self.verbose
891    
892                  # assume p is known: then v=A^-1(f-B^*p)
893                  # which leads to BA^-1B^*p = BA^-1f  
894    
895              # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)      
896    
897              
898              self.__z=v+self.solve_A(v,p*0)
899    
900                  Bz=self.B(self.__z)
901                  #
902              #   solve BA^-1B^*p = Bz
903                  #
904                  #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
905                  #
906                  #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
907                  #                           A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)
908                  #
909                  self.iter=0
910              if solver=='GMRES':      
911                    if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
912                    p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1.)
913                    # solve Au=f-B^*p
914                    #       A(u-v)=f-B^*p-Av
915                    #       u=v+(u-v)
916            u=v+self.solve_A(v,p)
917        
918                  else:
919                    if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
920                    p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
921                u=r[0]  
922                  print "div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
923    
924              return u,p
925    
926          def __Msolve(self,r):
927              return self.solve_prec(r[1])
928    
929          def __Msolve_GMRES(self,r):
930              return self.solve_prec(r)
931    
932    
933          def __Aprod(self,p):
934              # return BA^-1B*p
935              #solve Av =-B^*p as Av =f-Az-B^*p
936              v=self.solve_A(self.__z,-p)
937              return ArithmeticTuple(v, self.B(v))
938    
939          def __Aprod_GMRES(self,p):
940              # return BA^-1B*p
941              #solve Av =-B^*p as Av =f-Az-B^*p
942          v=self.solve_A(self.__z,-p)
943              return self.B(v)
944    
945  class SaddlePointProblem(object):  class SaddlePointProblem(object):
946     """     """
# Line 771  def MaskFromBoundaryTag(function_space,* Line 1196  def MaskFromBoundaryTag(function_space,*
1196     else:     else:
1197        return util.whereNonZero(util.interpolate(out,function_space))        return util.whereNonZero(util.interpolate(out,function_space))
1198    
1199    
1200    

Legend:
Removed from v.1312  
changed lines
  Added in v.1469

  ViewVC Help
Powered by ViewVC 1.1.26