/[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 1105 by gross, Thu Apr 19 01:10:49 2007 UTC revision 1552 by gross, Thu May 8 08:52:41 2008 UTC
# Line 1  Line 1 
1    #
2  # $Id$  # $Id$
3    #
4    #######################################################
5    #
6    #           Copyright 2003-2007 by ACceSS MNRF
7    #       Copyright 2007 by University of Queensland
8    #
9    #                http://esscc.uq.edu.au
10    #        Primary Business: Queensland, Australia
11    #  Licensed under the Open Software License version 3.0
12    #     http://www.opensource.org/licenses/osl-3.0.php
13    #
14    #######################################################
15    #
16    
17  """  """
18  Provides some tools related to PDEs.  Provides some tools related to PDEs.
# Line 32  import escript Line 46  import escript
46  import linearPDEs  import linearPDEs
47  import numarray  import numarray
48  import util  import util
49    import math
50    
51    ##### Added by Artak
52    # from Numeric import zeros,Int,Float64
53    ###################################
54    
55    
56  class TimeIntegrationManager:  class TimeIntegrationManager:
57    """    """
# Line 132  class Projector: Line 152  class Projector:
152      self.__pde.setValue(D = 1.)      self.__pde.setValue(D = 1.)
153      return      return
154    
   def __del__(self):  
     return  
   
155    def __call__(self, input_data):    def __call__(self, input_data):
156      """      """
157      Projects input_data onto a continuous function      Projects input_data onto a continuous function
# Line 142  class Projector: Line 159  class Projector:
159      @param input_data: The input_data to be projected.      @param input_data: The input_data to be projected.
160      """      """
161      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
162        self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())
163      if input_data.getRank()==0:      if input_data.getRank()==0:
164          self.__pde.setValue(Y = input_data)          self.__pde.setValue(Y = input_data)
165          out=self.__pde.getSolution()          out=self.__pde.getSolution()
# Line 395  class Locator: Line 413  class Locator:
413          else:          else:
414             return data             return data
415    
416    class SolverSchemeException(Exception):
417       """
418       exceptions thrown by solvers
419       """
420       pass
421    
422    class IndefinitePreconditioner(SolverSchemeException):
423       """
424       the preconditioner is not positive definite.
425       """
426       pass
427    class MaxIterReached(SolverSchemeException):
428       """
429       maxium number of iteration steps is reached.
430       """
431       pass
432    class IterationBreakDown(SolverSchemeException):
433       """
434       iteration scheme econouters an incurable breakdown.
435       """
436       pass
437    class NegativeNorm(SolverSchemeException):
438       """
439       a norm calculation returns a negative norm.
440       """
441       pass
442    
443    class IterationHistory(object):
444       """
445       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          @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          
468           @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,solver="GMRES",TOL=None):
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           if TOL==None:
494              TOL=self.tolerance
495           self.history.append(norm_r)
496           if self.verbose: print "iter: #s:  norm(r) = #e"#(len(self.history)-1, self.history[-1])
497           return self.history[-1]<=TOL * norm_b
498    
499    def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
500       """
501       Solver for
502    
503       M{Ax=b}
504    
505       with a symmetric and positive definite operator A (more details required!).
506       It uses the conjugate gradient method with preconditioner M providing an approximation of A.
507    
508       The iteration is terminated if the C{stoppingcriterium} function return C{True}.
509    
510       For details on the preconditioned conjugate gradient method see the book:
511    
512       Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
513       T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
514       C. Romine, and H. van der Vorst.
515    
516       @param b: the right hand side of the liner system. C{b} is altered.
517       @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
518       @param Aprod: returns the value Ax
519       @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}.
520       @param Msolve: solves Mx=r
521       @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
522    type like argument C{x}.
523       @param bilinearform: inner product C{<x,r>}
524       @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}.
525       @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.
526       @type stoppingcriterium: function that returns C{True} or C{False}
527       @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
528       @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
529       @param iter_max: maximum number of iteration steps.
530       @type iter_max: C{int}
531       @return: the solution approximation and the corresponding residual
532       @rtype: C{tuple}
533       @warning: C{b} and C{x} are altered.
534       """
535       iter=0
536       if x==None:
537          x=0*b
538       else:
539          b += (-1)*Aprod(x)
540       r=b
541       rhat=Msolve(r)
542       d = rhat
543       rhat_dot_r = bilinearform(rhat, r)
544       if rhat_dot_r<0: raise NegativeNorm,"negative norm."
545    
546       while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
547           iter+=1
548           if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
549    
550           q=Aprod(d)
551           alpha = rhat_dot_r / bilinearform(d, q)
552           x += alpha * d
553           r += (-alpha) * q
554    
555           rhat=Msolve(r)
556           rhat_dot_r_new = bilinearform(rhat, r)
557           beta = rhat_dot_r_new / rhat_dot_r
558           rhat+=beta * d
559           d=rhat
560    
561           rhat_dot_r = rhat_dot_r_new
562           if rhat_dot_r<0: raise NegativeNorm,"negative norm."
563    
564       return x,r
565    
566    
567    ############################
568    # Added by Artak
569    #################################3
570    
571    #Apply a sequence of k Givens rotations, used within gmres codes
572    # vrot=givapp(c, s, vin, k)
573    def givapp(c,s,vin):
574        vrot=vin # warning: vin is altered!!!!
575        if isinstance(c,float):
576            vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
577        else:
578            for i in range(len(c)):
579                w1=c[i]*vrot[i]-s[i]*vrot[i+1]
580            w2=s[i]*vrot[i]+c[i]*vrot[i+1]
581                vrot[i:i+2]=w1,w2
582        return vrot
583    
584    ##############################################
585    def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
586    ################################################
587       m=iter_restart
588       iter=0
589       xc=x
590       while True:
591          if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
592          xc,stopped=GMRESm(b*1, Aprod, Msolve, bilinearform, stoppingcriterium, x=xc*1, iter_max=iter_max-iter, iter_restart=m)
593          if stopped: break
594          iter+=iter_restart    
595       return xc
596    
597    def GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
598       iter=0
599       r=Msolve(b)
600       r_dot_r = bilinearform(r, r)
601       if r_dot_r<0: raise NegativeNorm,"negative norm."
602       norm_b=math.sqrt(r_dot_r)
603    
604       if x==None:
605          x=0*b
606       else:
607          r=Msolve(b-Aprod(x))
608          r_dot_r = bilinearform(r, r)
609          if r_dot_r<0: raise NegativeNorm,"negative norm."
610      
611       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
612       c=numarray.zeros(iter_restart,numarray.Float64)
613       s=numarray.zeros(iter_restart,numarray.Float64)
614       g=numarray.zeros(iter_restart,numarray.Float64)
615       v=[]
616    
617       rho=math.sqrt(r_dot_r)
618      
619       v.append(r/rho)
620       g[0]=rho
621       while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
622    
623        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
624    
625        
626        p=Msolve(Aprod(v[iter]))
627    
628        v.append(p)
629    
630        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
631    
632    # Modified Gram-Schmidt
633        for j in range(iter+1):
634          h[j][iter]=bilinearform(v[j],v[iter+1])  
635          v[iter+1]-=h[j][iter]*v[j]
636          
637        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
638        v_norm2=h[iter+1][iter]
639    
640    
641    # Reorthogonalize if needed
642        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
643         for j in range(iter+1):
644            hr=bilinearform(v[j],v[iter+1])
645                h[j][iter]=h[j][iter]+hr #vhat
646                v[iter+1] -= hr*v[j]
647    
648         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
649         h[iter+1][iter]=v_norm2
650    
651    #   watch out for happy breakdown
652            if not v_norm2 == 0:
653             v[iter+1]=v[iter+1]/h[iter+1][iter]
654    
655    #   Form and store the information for the new Givens rotation
656        if iter > 0 :
657            hhat=numarray.zeros(iter+1,numarray.Float64)
658            for i in range(iter+1) : hhat[i]=h[i][iter]
659            hhat=givapp(c[0:iter],s[0:iter],hhat);
660                for i in range(iter+1) : h[i][iter]=hhat[i]
661    
662        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
663        if mu!=0 :
664            c[iter]=h[iter][iter]/mu
665            s[iter]=-h[iter+1][iter]/mu
666            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
667            h[iter+1][iter]=0.0
668            g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
669    
670    # Update the residual norm
671            rho=abs(g[iter+1])
672        iter+=1
673    
674    # At this point either iter > iter_max or rho < tol.
675    # It's time to compute x and leave.        
676    
677       if iter > 0 :
678         y=numarray.zeros(iter,numarray.Float64)    
679         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
680         if iter > 1 :  
681            i=iter-2  
682            while i>=0 :
683              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
684              i=i-1
685         xhat=v[iter-1]*y[iter-1]
686         for i in range(iter-1):
687        xhat += v[i]*y[i]
688       else : xhat=v[0]
689        
690       x += xhat
691       if iter<iter_restart-1:
692          stopped=True
693       else:
694          stopped=False
695    
696       return x,stopped
697    
698    
699    ######################################################
700    def dirder(x, w, bilinearform, Aprod, Msolve, f0, b ):
701    ######################################################
702    
703    # DIRDER estimates the directional derivative of a function F.
704    
705    
706    # Hardwired difference increment.
707    #
708      epsnew = 1.0e-07
709    #
710    #  Scale the step.
711    #
712      norm_w=math.sqrt(bilinearform(w,w))
713      if norm_w== 0.0:
714        return x/x
715    
716      epsnew = epsnew / norm_w
717    
718      if norm_w > 0.0:
719        epsnew = epsnew * math.sqrt(bilinearform(x,x))
720    #
721    #  DEL and F1 could share the same space if storage
722    #  is more important than clarity.
723    #
724    
725      DEL = x + epsnew * w
726      f1 = -Msolve(Aprod(DEL))
727      z = ( f1 - f0 ) / epsnew
728      return z
729    
730    ######################################################
731    def FDGMRES(f0, Aprod, Msolve, bilinearform, stoppingcriterium, xc=None, x=None, iter_max=100, iter_restart=20,TOL=None):
732    ######################################################
733       b=-f0
734       b_dot_b = bilinearform(b, b)
735       if b_dot_b<0: raise NegativeNorm,"negative norm."
736       norm_b=math.sqrt(b_dot_b)
737    
738       r=b
739    
740       if x==None:
741          x=0*f0
742       else:
743          r=-dirder(xc,x,bilinearform,Aprod,Msolve,f0,b)-f0  
744          
745       r_dot_r = bilinearform(r, r)
746       if r_dot_r<0: raise NegativeNorm,"negative norm."
747      
748       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
749       c=numarray.zeros(iter_restart,numarray.Float64)
750       s=numarray.zeros(iter_restart,numarray.Float64)
751       g=numarray.zeros(iter_restart,numarray.Float64)
752       v=[]
753    
754       rho=math.sqrt(r_dot_r)
755      
756       v.append(r/rho)
757       g[0]=rho
758       iter=0
759    
760       while not (stoppingcriterium(rho,norm_b,solver="FDGMRES",TOL=TOL) or iter==iter_restart-1):
761    
762        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
763    
764        
765            p=dirder(xc, v[iter], bilinearform,Aprod,Msolve,f0,b)
766    
767        v.append(p)
768    
769        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
770    
771    # Modified Gram-Schmidt
772        for j in range(iter+1):
773          h[j][iter]=bilinearform(v[j],v[iter+1])  
774          v[iter+1]+=(-1.)*h[j][iter]*v[j]
775          
776        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
777        v_norm2=h[iter+1][iter]
778    
779    
780    # Reorthogonalize if needed
781        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
782         for j in range(iter+1):
783            hr=bilinearform(v[j],v[iter+1])
784                h[j][iter]=h[j][iter]+hr #vhat
785                v[iter+1] +=(-1.)*hr*v[j]
786    
787         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
788         h[iter+1][iter]=v_norm2
789    
790    #   watch out for happy breakdown
791            if v_norm2 != 0:
792             v[iter+1]=v[iter+1]/h[iter+1][iter]
793    
794    #   Form and store the information for the new Givens rotation
795        if iter > 0 :
796            hhat=[]
797            for i in range(iter+1) : hhat.append(h[i][iter])
798            hhat=givapp(c[0:iter],s[0:iter],hhat);
799                for i in range(iter+1) : h[i][iter]=hhat[i]
800    
801        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
802        if mu!=0 :
803            c[iter]=h[iter][iter]/mu
804            s[iter]=-h[iter+1][iter]/mu
805            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
806            h[iter+1][iter]=0.0
807            g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
808    
809    # Update the residual norm
810            rho=abs(g[iter+1])
811        iter+=1
812    
813       if iter > 0 :
814         y=numarray.zeros(iter,numarray.Float64)    
815         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
816         if iter > 1 :  
817            i=iter-2  
818            while i>=0 :
819              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
820              i=i-1
821         xhat=v[iter-1]*y[iter-1]
822         for i in range(iter-1):
823        xhat += v[i]*y[i]
824       else : xhat=v[0]
825        
826       x += xhat
827       if iter<iter_restart-1:
828          stopped=True
829       else:
830          stopped=False
831    
832       return x,stopped
833    
834    #################################################
835    def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
836    #################################################
837        #
838        #  minres solves the system of linear equations Ax = b
839        #  where A is a symmetric matrix (possibly indefinite or singular)
840        #  and b is a given vector.
841        #  
842        #  "A" may be a dense or sparse matrix (preferably sparse!)
843        #  or the name of a function such that
844        #               y = A(x)
845        #  returns the product y = Ax for any given vector x.
846        #
847        #  "M" defines a positive-definite preconditioner M = C C'.
848        #  "M" may be a dense or sparse matrix (preferably sparse!)
849        #  or the name of a function such that
850        #  solves the system My = x for any given vector x.
851        #
852        #
853        
854        #------------------------------------------------------------------
855        # Set up y and v for the first Lanczos vector v1.
856        # y  =  beta1 P' v1,  where  P = C**(-1).
857        # v is really P' v1.
858        #------------------------------------------------------------------
859        if x==None:
860          x=0*b
861        else:
862          b += (-1)*Aprod(x)
863    
864        r1    = b
865        y = Msolve(b)
866        beta1 = bilinearform(y,b)
867    
868        if beta1< 0: raise NegativeNorm,"negative norm."
869    
870        #  If b = 0 exactly, stop with x = 0.
871        if beta1==0: return x*0.
872    
873        if beta1> 0:
874          beta1  = math.sqrt(beta1)      
875    
876        #------------------------------------------------------------------
877        # Initialize quantities.
878        # ------------------------------------------------------------------
879        iter   = 0
880        Anorm = 0
881        ynorm = 0
882        oldb   = 0
883        beta   = beta1
884        dbar   = 0
885        epsln  = 0
886        phibar = beta1
887        rhs1   = beta1
888        rhs2   = 0
889        rnorm  = phibar
890        tnorm2 = 0
891        ynorm2 = 0
892        cs     = -1
893        sn     = 0
894        w      = b*0.
895        w2     = b*0.
896        r2     = r1
897        eps    = 0.0001
898    
899        #---------------------------------------------------------------------
900        # Main iteration loop.
901        # --------------------------------------------------------------------
902        while not stoppingcriterium(rnorm,Anorm*ynorm,'MINRES'):    #  checks ||r|| < (||A|| ||x||) * TOL
903    
904        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
905            iter    = iter  +  1
906    
907            #-----------------------------------------------------------------
908            # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
909            # The general iteration is similar to the case k = 1 with v0 = 0:
910            #
911            #   p1      = Operator * v1  -  beta1 * v0,
912            #   alpha1  = v1'p1,
913            #   q2      = p2  -  alpha1 * v1,
914            #   beta2^2 = q2'q2,
915            #   v2      = (1/beta2) q2.
916            #
917            # Again, y = betak P vk,  where  P = C**(-1).
918            #-----------------------------------------------------------------
919            s = 1/beta                 # Normalize previous vector (in y).
920            v = s*y                    # v = vk if P = I
921        
922            y      = Aprod(v)
923        
924            if iter >= 2:
925              y = y - (beta/oldb)*r1
926    
927            alfa   = bilinearform(v,y)              # alphak
928            y      += (- alfa/beta)*r2
929            r1     = r2
930            r2     = y
931            y = Msolve(r2)
932            oldb   = beta                           # oldb = betak
933            beta   = bilinearform(y,r2)             # beta = betak+1^2
934            if beta < 0: raise NegativeNorm,"negative norm."
935    
936            beta   = math.sqrt( beta )
937            tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
938            
939            if iter==1:                 # Initialize a few things.
940              gmax   = abs( alfa )      # alpha1
941              gmin   = gmax             # alpha1
942    
943            # Apply previous rotation Qk-1 to get
944            #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
945            #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
946        
947            oldeps = epsln
948            delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
949            gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
950            epsln  =               sn * beta  # epsln2 = 0         epslnk+1
951            dbar   =            -  cs * beta  # dbar 2 = beta2     dbar k+1
952    
953            # Compute the next plane rotation Qk
954    
955            gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
956            gamma  = max(gamma,eps)
957            cs     = gbar / gamma             # ck
958            sn     = beta / gamma             # sk
959            phi    = cs * phibar              # phik
960            phibar = sn * phibar              # phibark+1
961    
962            # Update  x.
963    
964            denom = 1/gamma
965            w1    = w2
966            w2    = w
967            w     = (v - oldeps*w1 - delta*w2) * denom
968            x     +=  phi*w
969    
970            # Go round again.
971    
972            gmax   = max(gmax,gamma)
973            gmin   = min(gmin,gamma)
974            z      = rhs1 / gamma
975            ynorm2 = z*z  +  ynorm2
976            rhs1   = rhs2 -  delta*z
977            rhs2   =      -  epsln*z
978    
979            # Estimate various norms and test for convergence.
980    
981            Anorm  = math.sqrt( tnorm2 )
982            ynorm  = math.sqrt( ynorm2 )
983    
984            rnorm  = phibar
985    
986        return x
987    
988    ######################################    
989    def NewtonGMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium,x=None, iter_max=100,iter_restart=20,atol=1.e-2,rtol=1.e-4):
990    #####################################
991        gamma=.9
992        lmaxit=100
993        etamax=.5
994    
995        n = 1 #len(x)
996        iter=0
997        
998        # evaluate f at the initial iterate
999        # compute the stop tolerance
1000        #
1001        r=b
1002        if x==None:
1003          x=0*b
1004        else:
1005          r += (-1)*Aprod(x)
1006    
1007        f0=-Msolve(r)
1008        fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
1009        fnrmo=1
1010        stop_tol=atol + rtol*fnrm
1011        #
1012        # main iteration loop
1013        #
1014        while not stoppingcriterium(fnrm*1,stop_tol,'NewtonGMRES',TOL=1.):
1015    
1016                if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1017            #
1018            # keep track of the ratio (rat = fnrm/frnmo)
1019            # of successive residual norms and
1020            # the iteration counter (iter)
1021            #
1022            #rat=fnrm/fnrmo
1023            fnrmo=fnrm
1024            iter+=1
1025            #
1026                # compute the step using a GMRES(m) routine especially designed for this purpose
1027            #
1028                initer=0
1029                while True:
1030                   xc,stopped=FDGMRES(f0*1, Aprod, Msolve, bilinearform, stoppingcriterium, xc=x*1, iter_max=lmaxit-initer, iter_restart=iter_restart, TOL=etamax)
1031                   if stopped: break
1032                   initer+=iter_restart
1033            x+=xc
1034            f0=-Msolve(Aprod(x))
1035            fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
1036            rat=fnrm/fnrmo
1037    
1038    
1039        #   adjust eta
1040        #
1041            if etamax > 0:
1042                etaold=etamax
1043                etanew=gamma*rat*rat
1044                if gamma*etaold*etaold > .1 :
1045                    etanew=max(etanew,gamma*etaold*etaold)
1046                
1047                etamax=min(etanew,etamax)
1048                etamax=max(etamax,.5*stop_tol/fnrm)
1049    
1050        return x
1051    
1052    def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
1053    
1054    # TFQMR solver for linear systems
1055    #
1056    #
1057    # initialization
1058    #
1059      errtol = math.sqrt(bilinearform(b,b))
1060      norm_b=errtol
1061      kmax  = iter_max
1062      error = []
1063    
1064      if math.sqrt(bilinearform(x,x)) != 0.0:
1065        r = b - Aprod(x)
1066      else:
1067        r = b
1068    
1069      r=Msolve(r)
1070    
1071      u1=0
1072      u2=0
1073      y1=0
1074      y2=0
1075    
1076      w = r
1077      y1 = r
1078      iter = 0
1079      d = 0
1080      
1081      v = Msolve(Aprod(y1))
1082      u1 = v
1083      
1084      theta = 0.0;
1085      eta = 0.0;
1086      tau = math.sqrt(bilinearform(r,r))
1087      error = [ error, tau ]
1088      rho = tau * tau
1089      m=1
1090    #
1091    #  TFQMR iteration
1092    #
1093    #  while ( iter < kmax-1 ):
1094      
1095      while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):
1096        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1097    
1098        sigma = bilinearform(r,v)
1099    
1100        if ( sigma == 0.0 ):
1101          raise 'TFQMR breakdown, sigma=0'
1102        
1103    
1104        alpha = rho / sigma
1105    
1106        for j in range(2):
1107    #
1108    #   Compute y2 and u2 only if you have to
1109    #
1110          if ( j == 1 ):
1111            y2 = y1 - alpha * v
1112            u2 = Msolve(Aprod(y2))
1113    
1114          m = 2 * (iter+1) - 2 + (j+1)
1115          if j==0:
1116             w = w - alpha * u1
1117             d = y1 + ( theta * theta * eta / alpha ) * d
1118          if j==1:
1119             w = w - alpha * u2
1120             d = y2 + ( theta * theta * eta / alpha ) * d
1121    
1122          theta = math.sqrt(bilinearform(w,w))/ tau
1123          c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1124          tau = tau * theta * c
1125          eta = c * c * alpha
1126          x = x + eta * d
1127    #
1128    #  Try to terminate the iteration at each pass through the loop
1129    #
1130         # if ( tau * math.sqrt ( m + 1 ) <= errtol ):
1131         #   error = [ error, tau ]
1132         #   total_iters = iter
1133         #   break
1134          
1135    
1136        if ( rho == 0.0 ):
1137          raise 'TFQMR breakdown, rho=0'
1138        
1139    
1140        rhon = bilinearform(r,w)
1141        beta = rhon / rho;
1142        rho = rhon;
1143        y1 = w + beta * y2;
1144        u1 = Msolve(Aprod(y1))
1145        v = u1 + beta * ( u2 + beta * v )
1146        error = [ error, tau ]
1147        total_iters = iter
1148        
1149        iter = iter + 1
1150    
1151      return x
1152    
1153    
1154    #############################################
1155    
1156    class ArithmeticTuple(object):
1157       """
1158       tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
1159    
1160       example of usage:
1161    
1162       from esys.escript import Data
1163       from numarray import array
1164       a=Data(...)
1165       b=array([1.,4.])
1166       x=ArithmeticTuple(a,b)
1167       y=5.*x
1168    
1169       """
1170       def __init__(self,*args):
1171           """
1172           initialize object with elements args.
1173    
1174           @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
1175           """
1176           self.__items=list(args)
1177    
1178       def __len__(self):
1179           """
1180           number of items
1181    
1182           @return: number of items
1183           @rtype: C{int}
1184           """
1185           return len(self.__items)
1186    
1187       def __getitem__(self,index):
1188           """
1189           get an item
1190    
1191           @param index: item to be returned
1192           @type index: C{int}
1193           @return: item with index C{index}
1194           """
1195           return self.__items.__getitem__(index)
1196    
1197       def __mul__(self,other):
1198           """
1199           scaling from the right
1200    
1201           @param other: scaling factor
1202           @type other: C{float}
1203           @return: itemwise self*other
1204           @rtype: L{ArithmeticTuple}
1205           """
1206           out=[]
1207           for i in range(len(self)):
1208               out.append(self[i]*other)
1209           return ArithmeticTuple(*tuple(out))
1210    
1211       def __rmul__(self,other):
1212           """
1213           scaling from the left
1214    
1215           @param other: scaling factor
1216           @type other: C{float}
1217           @return: itemwise other*self
1218           @rtype: L{ArithmeticTuple}
1219           """
1220           out=[]
1221           for i in range(len(self)):
1222               out.append(other*self[i])
1223           return ArithmeticTuple(*tuple(out))
1224    
1225    #########################
1226    # Added by Artak
1227    #########################
1228       def __div__(self,other):
1229           """
1230           dividing from the right
1231    
1232           @param other: scaling factor
1233           @type other: C{float}
1234           @return: itemwise self/other
1235           @rtype: L{ArithmeticTuple}
1236           """
1237           out=[]
1238           for i in range(len(self)):
1239               out.append(self[i]*(1./other))
1240           return ArithmeticTuple(*tuple(out))
1241    
1242       def __rdiv__(self,other):
1243           """
1244           dividing from the left
1245    
1246           @param other: scaling factor
1247           @type other: C{float}
1248           @return: itemwise other/self
1249           @rtype: L{ArithmeticTuple}
1250           """
1251           out=[]
1252           for i in range(len(self)):
1253               out.append(other/self[i])
1254           return ArithmeticTuple(*tuple(out))
1255      
1256    ##########################################33
1257    
1258       def __iadd__(self,other):
1259           """
1260           in-place add of other to self
1261    
1262           @param other: increment
1263           @type other: C{ArithmeticTuple}
1264           """
1265           if len(self) != len(other):
1266               raise ValueError,"tuple length must match."
1267           for i in range(len(self)):
1268               self.__items[i]+=other[i]
1269           return self
1270    
1271       def __add__(self,other):
1272           """
1273           add other to self
1274    
1275           @param other: increment
1276           @type other: C{ArithmeticTuple}
1277           """
1278           if not isinstance(other,float):
1279            if len(self) != len(other):
1280               raise ValueError,"tuple length must match."
1281            for i in range(len(self)):
1282               self.__items[i]+=other[i]
1283           else:
1284            for i in range(len(self)):
1285               self.__items[i]+=other
1286    
1287           return self
1288    
1289       def __sub__(self,other):
1290           """
1291           subtract other from self
1292    
1293           @param other: increment
1294           @type other: C{ArithmeticTuple}
1295           """
1296           if len(self) != len(other):
1297               raise ValueError,"tuple length must match."
1298           for i in range(len(self)):
1299               self.__items[i]-=other[i]
1300           return self
1301    
1302       def __neg__(self):
1303           """
1304           negate
1305    
1306           """
1307           for i in range(len(self)):
1308               self.__items[i]=-self.__items[i]
1309           return self
1310    
1311    
1312    class HomogeneousSaddlePointProblem(object):
1313          """
1314          This provides a framwork for solving homogeneous saddle point problem of the form
1315    
1316                 Av+B^*p=f
1317                 Bv    =0
1318    
1319          for the unknowns v and p and given operators A and B and given right hand side f.
1320          B^* is the adjoint operator of B is the given inner product.
1321    
1322          """
1323          def __init__(self,**kwargs):
1324            self.setTolerance()
1325            self.setToleranceReductionFactor()
1326    
1327          def initialize(self):
1328            """
1329            initialize the problem (overwrite)
1330            """
1331            pass
1332          def B(self,v):
1333             """
1334             returns Bv (overwrite)
1335             @rtype: equal to the type of p
1336    
1337             @note: boundary conditions on p should be zero!
1338             """
1339             pass
1340    
1341          def inner(self,p0,p1):
1342             """
1343             returns inner product of two element p0 and p1  (overwrite)
1344            
1345             @type p0: equal to the type of p
1346             @type p1: equal to the type of p
1347             @rtype: C{float}
1348    
1349             @rtype: equal to the type of p
1350             """
1351             pass
1352    
1353          def solve_A(self,u,p):
1354             """
1355             solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
1356    
1357             @rtype: equal to the type of v
1358             @note: boundary conditions on v should be zero!
1359             """
1360             pass
1361    
1362          def solve_prec(self,p):
1363             """
1364             provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
1365    
1366             @rtype: equal to the type of p
1367             """
1368             pass
1369    
1370          def stoppingcriterium(self,Bv,v,p):
1371             """
1372             returns a True if iteration is terminated. (overwrite)
1373    
1374             @rtype: C{bool}
1375             """
1376             pass
1377                
1378          def __inner(self,p,r):
1379             return self.inner(p,r[1])
1380    
1381          def __inner_p(self,p1,p2):
1382             return self.inner(p1,p2)
1383          
1384          def __inner_a(self,a1,a2):
1385             return self.inner_a(a1,a2)
1386    
1387          def __inner_a1(self,a1,a2):
1388             return self.inner(a1[1],a2[1])
1389    
1390          def __stoppingcriterium(self,norm_r,r,p):
1391              return self.stoppingcriterium(r[1],r[0],p)
1392    
1393          def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
1394              return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)
1395    
1396          def setTolerance(self,tolerance=1.e-8):
1397                  self.__tol=tolerance
1398          def getTolerance(self):
1399                  return self.__tol
1400          def setToleranceReductionFactor(self,reduction=0.01):
1401                  self.__reduction=reduction
1402          def getSubProblemTolerance(self):
1403                  return self.__reduction*self.getTolerance()
1404    
1405          def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):
1406                  """
1407                  solves the saddle point problem using initial guesses v and p.
1408    
1409                  @param max_iter: maximum number of iteration steps.
1410                  """
1411                  self.verbose=verbose
1412                  self.show_details=show_details and self.verbose
1413    
1414                  # assume p is known: then v=A^-1(f-B^*p)
1415                  # which leads to BA^-1B^*p = BA^-1f  
1416    
1417              # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)      
1418              self.__z=v+self.solve_A(v,p*0)
1419                  Bz=self.B(self.__z)
1420                  #
1421              #   solve BA^-1B^*p = Bz
1422                  #
1423                  #
1424                  #
1425                  self.iter=0
1426              if solver=='GMRES':      
1427                    if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1428                    p=GMRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)
1429                    # solve Au=f-B^*p
1430                    #       A(u-v)=f-B^*p-Av
1431                    #       u=v+(u-v)
1432            u=v+self.solve_A(v,p)
1433    
1434              if solver=='NewtonGMRES':    
1435                    if self.verbose: print "enter NewtonGMRES method (iter_max=%s)"%max_iter
1436                    p=NewtonGMRES(Bz,self.__Aprod_Newton,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,atol=0,rtol=self.getTolerance())
1437                    # solve Au=f-B^*p
1438                    #       A(u-v)=f-B^*p-Av
1439                    #       u=v+(u-v)
1440            u=v+self.solve_A(v,p)
1441                    
1442    
1443              if solver=='TFQMR':      
1444                    if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter
1445                    p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1446                    # solve Au=f-B^*p
1447                    #       A(u-v)=f-B^*p-Av
1448                    #       u=v+(u-v)
1449            u=v+self.solve_A(v,p)
1450    
1451              if solver=='MINRES':      
1452                    if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1453                    p=MINRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1454                    # solve Au=f-B^*p
1455                    #       A(u-v)=f-B^*p-Av
1456                    #       u=v+(u-v)
1457            u=v+self.solve_A(v,p)
1458                  
1459              if solver=='GMRESC':      
1460                    if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter
1461                    p0=self.solve_prec(Bz)
1462               #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))
1463                   #p-=alfa
1464                    x=GMRES(ArithmeticTuple(self.__z*1.,p0*1),self.__Anext,self.__Mempty,self.__inner_a,self.__stoppingcriterium2,iter_max=max_iter, x=ArithmeticTuple(v*1,p*1),iter_restart=20)
1465                    #x=NewtonGMRES(ArithmeticTuple(self.__z*1.,p0*1),self.__Aprod_Newton2,self.__Mempty,self.__inner_a,self.__stoppingcriterium2,iter_max=max_iter, x=ArithmeticTuple(v*1,p*1),atol=0,rtol=self.getTolerance())
1466                    # solve Au=f-B^*p
1467                    #       A(u-v)=f-B^*p-Av
1468                    #       u=v+(u-v)
1469                p=x[1]
1470            #u=v+self.solve_A(v,p)      
1471            #p=x[1]
1472            #u=x[0]
1473    
1474                  if solver=='PCG':
1475                    #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1476                    #
1477                    #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1478                    #                           A(v-z)= f -Az - B^*p (v-z=0 on fixed_u_mask)
1479                    if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1480                    p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1481                u=r[0]  
1482                    print "DIFF=",util.integrate(p)
1483    
1484                  print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1485    
1486              return u,p
1487    
1488          def __Msolve(self,r):
1489              return self.solve_prec(r[1])
1490    
1491          def __Msolve2(self,r):
1492              return self.solve_prec(r*1)
1493    
1494          def __Mempty(self,r):
1495              return r
1496    
1497    
1498          def __Aprod(self,p):
1499              # return BA^-1B*p
1500              #solve Av =B^*p as Av =f-Az-B^*(-p)
1501              v=self.solve_A(self.__z,-p)
1502              return ArithmeticTuple(v, self.B(v))
1503    
1504          def __Anext(self,x):
1505              # return next v,p
1506              #solve Av =-B^*p as Av =f-Az-B^*p
1507    
1508          pc=x[1]
1509              v=self.solve_A(self.__z,-pc)
1510          p=self.solve_prec(self.B(v))
1511    
1512              return ArithmeticTuple(v,p)
1513    
1514    
1515          def __Aprod2(self,p):
1516              # return BA^-1B*p
1517              #solve Av =B^*p as Av =f-Az-B^*(-p)
1518          v=self.solve_A(self.__z,-p)
1519              return self.B(v)
1520    
1521          def __Aprod_Newton(self,p):
1522              # return BA^-1B*p - Bz
1523              #solve Av =-B^*p as Av =f-Az-B^*p
1524          v=self.solve_A(self.__z,-p)
1525              return self.B(v-self.__z)
1526    
1527          def __Aprod_Newton2(self,x):
1528              # return BA^-1B*p - Bz
1529              #solve Av =-B^*p as Av =f-Az-B^*p
1530              pc=x[1]
1531          v=self.solve_A(self.__z,-pc)
1532              p=self.solve_prec(self.B(v-self.__z))
1533              return ArithmeticTuple(v,p)
1534    
1535  class SaddlePointProblem(object):  class SaddlePointProblem(object):
1536     """     """
1537     This implements a solver for a saddlepoint problem     This implements a solver for a saddlepoint problem
# Line 420  class SaddlePointProblem(object): Line 1557  class SaddlePointProblem(object):
1557         @type verbose: C{bool}         @type verbose: C{bool}
1558         @note: this method may be overwritten by a particular saddle point problem         @note: this method may be overwritten by a particular saddle point problem
1559         """         """
1560         if verbose:         if not isinstance(verbose,bool):
1561            self.__verbose=True              raise TypeError("verbose needs to be of type bool.")
1562         else:         self.__verbose=verbose
           self.__verbose=False  
1563         self.relaxation=1.         self.relaxation=1.
1564    
1565     def trace(self,text):     def trace(self,text):
# Line 433  class SaddlePointProblem(object): Line 1569  class SaddlePointProblem(object):
1569         @param text: a text message         @param text: a text message
1570         @type text: C{str}         @type text: C{str}
1571         """         """
1572         if self.__verbose: print "%s: %s"%(str(self),text)         if self.__verbose: print "#s: #s"%(str(self),text)
1573    
1574     def solve_f(self,u,p,tol=1.e-8):     def solve_f(self,u,p,tol=1.e-8):
1575         """         """
# Line 583  class SaddlePointProblem(object): Line 1719  class SaddlePointProblem(object):
1719              norm_u_new = util.Lsup(u_new)              norm_u_new = util.Lsup(u_new)
1720              p_new=p+self.relaxation*g_new              p_new=p+self.relaxation*g_new
1721              norm_p_new = util.sqrt(self.inner(p_new,p_new))              norm_p_new = util.sqrt(self.inner(p_new,p_new))
1722              self.trace("%s th step: f/u = %s, g/p = %s, relaxation = %s."%(self.iter,norm_f_new/norm_u_new, norm_g_new/norm_p_new, self.relaxation))              self.trace("%s th step: f = %s, f/u = %s, g = %s, g/p = %s, relaxation = %s."%(self.iter, norm_f_new ,norm_f_new/norm_u_new, norm_g_new, norm_g_new/norm_p_new, self.relaxation))
1723    
1724              if self.iter>1:              if self.iter>1:
1725                 dg2=g_new-g                 dg2=g_new-g
# Line 625  class SaddlePointProblem(object): Line 1761  class SaddlePointProblem(object):
1761  #  #
1762  #      return u,p  #      return u,p
1763                        
1764  # vim: expandtab shiftwidth=4:  def MaskFromBoundaryTag(function_space,*tags):
1765       """
1766       create a mask on the given function space which one for samples
1767       that touch the boundary tagged by tags.
1768    
1769       usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
1770    
1771       @param function_space: a given function space
1772       @type function_space: L{escript.FunctionSpace}
1773       @param tags: boundray tags
1774       @type tags: C{str}
1775       @return: a mask which marks samples used by C{function_space} that are touching the
1776                boundary tagged by any of the given tags.
1777       @rtype: L{escript.Data} of rank 0
1778       """
1779       pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
1780       d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
1781       for t in tags: d.setTaggedValue(t,1.)
1782       pde.setValue(y=d)
1783       out=util.whereNonZero(pde.getRightHandSide())
1784       if out.getFunctionSpace() == function_space:
1785          return out
1786       else:
1787          return util.whereNonZero(util.interpolate(out,function_space))
1788    
1789    
1790    

Legend:
Removed from v.1105  
changed lines
  Added in v.1552

  ViewVC Help
Powered by ViewVC 1.1.26