/[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 1519 by artak, Tue Apr 22 03:45:36 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    def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
585       m=iter_restart
586       iter=0
587       while True:
588          if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
589          x,stopped=GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=x, iter_max=iter_max-iter, iter_restart=m)
590          if stopped: break
591          iter+=iter_restart    
592       return x
593    
594    def GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
595       iter=0
596       r=Msolve(b)
597       r_dot_r = bilinearform(r, r)
598       if r_dot_r<0: raise NegativeNorm,"negative norm."
599       norm_b=math.sqrt(r_dot_r)
600    
601       if x==None:
602          x=0
603       else:
604          r=Msolve(b-Aprod(x))
605          r_dot_r = bilinearform(r, r)
606          if r_dot_r<0: raise NegativeNorm,"negative norm."
607      
608       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
609       c=numarray.zeros(iter_restart,numarray.Float64)
610       s=numarray.zeros(iter_restart,numarray.Float64)
611       g=numarray.zeros(iter_restart,numarray.Float64)
612       v=[]
613    
614       rho=math.sqrt(r_dot_r)
615      
616       v.append(r/rho)
617       g[0]=rho
618       while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
619    
620        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
621    
622        
623        p=Msolve(Aprod(v[iter]))
624    
625        v.append(p)
626    
627        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
628    
629    # Modified Gram-Schmidt
630        for j in range(iter+1):
631          h[j][iter]=bilinearform(v[j],v[iter+1])  
632          v[iter+1]+=(-1.)*h[j][iter]*v[j]
633          
634        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
635        v_norm2=h[iter+1][iter]
636    
637    
638    # Reorthogonalize if needed
639        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
640         for j in range(iter+1):
641            hr=bilinearform(v[j],v[iter+1])
642                h[j][iter]=h[j][iter]+hr #vhat
643                v[iter+1] +=(-1.)*hr*v[j]
644    
645         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
646         h[iter+1][iter]=v_norm2
647    
648    #   watch out for happy breakdown
649            if v_norm2 != 0:
650             v[iter+1]=v[iter+1]/h[iter+1][iter]
651    
652    #   Form and store the information for the new Givens rotation
653        if iter > 0 :
654            hhat=[]
655            for i in range(iter+1) : hhat.append(h[i][iter])
656            hhat=givapp(c[0:iter],s[0:iter],hhat);
657                for i in range(iter+1) : h[i][iter]=hhat[i]
658    
659        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
660        if mu!=0 :
661            c[iter]=h[iter][iter]/mu
662            s[iter]=-h[iter+1][iter]/mu
663            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
664            h[iter+1][iter]=0.0
665            g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
666    
667    # Update the residual norm
668            rho=abs(g[iter+1])
669        iter+=1
670    
671    # At this point either iter > iter_max or rho < tol.
672    # It's time to compute x and leave.        
673    
674       if iter > 0 :
675         y=numarray.zeros(iter,numarray.Float64)    
676         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
677         if iter > 1 :  
678            i=iter-2  
679            while i>=0 :
680              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
681              i=i-1
682         xhat=v[iter-1]*y[iter-1]
683         for i in range(iter-1):
684        xhat += v[i]*y[i]
685       else : xhat=v[0]
686        
687       x += xhat
688       if iter<iter_restart-1:
689          stopped=True
690       else:
691          stopped=False
692    
693       return x,stopped
694    
695    ######################################################
696    def dirder(x, w, bilinearform, Aprod, Msolve, f0, b ):
697    ######################################################
698    
699    # DIRDER estimates the directional derivative of a function F.
700    
701    
702    # Hardwired difference increment.
703    #
704      epsnew = 1.0e-07
705    #
706    #  Scale the step.
707    #
708      norm_w=math.sqrt(bilinearform(w,w))
709      if norm_w== 0.0:
710        return x/x
711    
712      epsnew = epsnew / norm_w
713    
714      if norm_w > 0.0:
715        epsnew = epsnew * math.sqrt(bilinearform(x,x))
716    #
717    #  DEL and F1 could share the same space if storage
718    #  is more important than clarity.
719    #
720    
721      DEL = x + epsnew * w
722      f1 = -Msolve(Aprod(DEL))
723      z = ( f1 - f0 ) / epsnew
724      return z
725    
726    ######################################################
727    def FDGMRES(f0, Aprod, Msolve, bilinearform, stoppingcriterium, xc=None, x=None, iter_max=100, iter_restart=20,TOL=None):
728    ######################################################
729       b=-f0
730       b_dot_b = bilinearform(b, b)
731       if b_dot_b<0: raise NegativeNorm,"negative norm."
732       norm_b=math.sqrt(b_dot_b)
733    
734       r=b
735    
736       if x==None:
737          x=0
738       else:
739          r=-dirder(xc,x,bilinearform,Aprod,Msolve,f0,b)-f0  
740          
741       r_dot_r = bilinearform(r, r)
742       if r_dot_r<0: raise NegativeNorm,"negative norm."
743      
744       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
745       c=numarray.zeros(iter_restart,numarray.Float64)
746       s=numarray.zeros(iter_restart,numarray.Float64)
747       g=numarray.zeros(iter_restart,numarray.Float64)
748       v=[]
749    
750       rho=math.sqrt(r_dot_r)
751      
752       v.append(r/rho)
753       g[0]=rho
754       iter=0
755    
756       while not (stoppingcriterium(rho,norm_b,solver="FDGMRES",TOL=TOL) or iter==iter_restart-1):
757    
758        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
759    
760        
761            p=dirder(xc, v[iter], bilinearform,Aprod,Msolve,f0,b)
762    
763        v.append(p)
764    
765        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
766    
767    # Modified Gram-Schmidt
768        for j in range(iter+1):
769          h[j][iter]=bilinearform(v[j],v[iter+1])  
770          v[iter+1]+=(-1.)*h[j][iter]*v[j]
771          
772        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
773        v_norm2=h[iter+1][iter]
774    
775    
776    # Reorthogonalize if needed
777        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
778         for j in range(iter+1):
779            hr=bilinearform(v[j],v[iter+1])
780                h[j][iter]=h[j][iter]+hr #vhat
781                v[iter+1] +=(-1.)*hr*v[j]
782    
783         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
784         h[iter+1][iter]=v_norm2
785    
786    #   watch out for happy breakdown
787            if v_norm2 != 0:
788             v[iter+1]=v[iter+1]/h[iter+1][iter]
789    
790    #   Form and store the information for the new Givens rotation
791        if iter > 0 :
792            hhat=[]
793            for i in range(iter+1) : hhat.append(h[i][iter])
794            hhat=givapp(c[0:iter],s[0:iter],hhat);
795                for i in range(iter+1) : h[i][iter]=hhat[i]
796    
797        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
798        if mu!=0 :
799            c[iter]=h[iter][iter]/mu
800            s[iter]=-h[iter+1][iter]/mu
801            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
802            h[iter+1][iter]=0.0
803            g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
804    
805    # Update the residual norm
806            rho=abs(g[iter+1])
807        iter+=1
808    
809       if iter > 0 :
810         y=numarray.zeros(iter,numarray.Float64)    
811         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
812         if iter > 1 :  
813            i=iter-2  
814            while i>=0 :
815              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
816              i=i-1
817         xhat=v[iter-1]*y[iter-1]
818         for i in range(iter-1):
819        xhat += v[i]*y[i]
820       else : xhat=v[0]
821        
822       x += xhat
823       if iter<iter_restart-1:
824          stopped=True
825       else:
826          stopped=False
827    
828       return x,stopped
829    
830    def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
831    
832        #
833        #  minres solves the system of linear equations Ax = b
834        #  where A is a symmetric matrix (possibly indefinite or singular)
835        #  and b is a given vector.
836        #  
837        #  "A" may be a dense or sparse matrix (preferably sparse!)
838        #  or the name of a function such that
839        #               y = A(x)
840        #  returns the product y = Ax for any given vector x.
841        #
842        #  "M" defines a positive-definite preconditioner M = C C'.
843        #  "M" may be a dense or sparse matrix (preferably sparse!)
844        #  or the name of a function such that
845        #  solves the system My = x for any given vector x.
846        #
847        #
848        
849        #------------------------------------------------------------------
850        # Set up y and v for the first Lanczos vector v1.
851        # y  =  beta1 P' v1,  where  P = C**(-1).
852        # v is really P' v1.
853        #------------------------------------------------------------------
854        if x==None:
855          x=0*b
856        else:
857          b += (-1)*Aprod(x)
858    
859        r1    = b
860        y = Msolve(b)
861        beta1 = bilinearform(b,y)
862    
863        if beta1< 0: raise NegativeNorm,"negative norm."
864    
865        #  If b = 0 exactly, stop with x = 0.
866        if beta1==0: return x*0.
867    
868        if beta1> 0:
869          beta1  = math.sqrt(beta1)      
870    
871        #------------------------------------------------------------------
872        # Initialize quantities.
873        # ------------------------------------------------------------------
874        iter   = 0
875        Anorm = 0
876        ynorm = 0
877        oldb   = 0
878        beta   = beta1
879        dbar   = 0
880        epsln  = 0
881        phibar = beta1
882        rhs1   = beta1
883        rhs2   = 0
884        rnorm  = phibar
885        tnorm2 = 0
886        ynorm2 = 0
887        cs     = -1
888        sn     = 0
889        w      = b*0.
890        w2     = b*0.
891        r2     = r1
892        eps    = 0.0001
893    
894        #---------------------------------------------------------------------
895        # Main iteration loop.
896        # --------------------------------------------------------------------
897        while not stoppingcriterium(rnorm,Anorm*ynorm,'MINRES'):    #  checks ||r|| < (||A|| ||x||) * TOL
898    
899        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
900            iter    = iter  +  1
901    
902            #-----------------------------------------------------------------
903            # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
904            # The general iteration is similar to the case k = 1 with v0 = 0:
905            #
906            #   p1      = Operator * v1  -  beta1 * v0,
907            #   alpha1  = v1'p1,
908            #   q2      = p2  -  alpha1 * v1,
909            #   beta2^2 = q2'q2,
910            #   v2      = (1/beta2) q2.
911            #
912            # Again, y = betak P vk,  where  P = C**(-1).
913            #-----------------------------------------------------------------
914            s = 1/beta                 # Normalize previous vector (in y).
915            v = s*y                    # v = vk if P = I
916        
917            y      = Aprod(v)
918        
919            if iter >= 2:
920              y = y - (beta/oldb)*r1
921    
922            alfa   = bilinearform(v,y)              # alphak
923            y      = (- alfa/beta)*r2 + y
924            r1     = r2
925            r2     = y
926            y = Msolve(r2)
927            oldb   = beta                           # oldb = betak
928            beta   = bilinearform(r2,y)             # beta = betak+1^2
929            if beta < 0: raise NegativeNorm,"negative norm."
930    
931            beta   = math.sqrt( beta )
932            tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
933            
934            if iter==1:                 # Initialize a few things.
935              gmax   = abs( alfa )      # alpha1
936              gmin   = gmax             # alpha1
937    
938            # Apply previous rotation Qk-1 to get
939            #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
940            #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
941        
942            oldeps = epsln
943            delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
944            gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
945            epsln  =               sn * beta  # epsln2 = 0         epslnk+1
946            dbar   =            -  cs * beta  # dbar 2 = beta2     dbar k+1
947    
948            # Compute the next plane rotation Qk
949    
950            gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
951            gamma  = max(gamma,eps)
952            cs     = gbar / gamma             # ck
953            sn     = beta / gamma             # sk
954            phi    = cs * phibar              # phik
955            phibar = sn * phibar              # phibark+1
956    
957            # Update  x.
958    
959            denom = 1/gamma
960            w1    = w2
961            w2    = w
962            w     = (v - oldeps*w1 - delta*w2) * denom
963            x     = x  +  phi*w
964    
965            # Go round again.
966    
967            gmax   = max(gmax,gamma)
968            gmin   = min(gmin,gamma)
969            z      = rhs1 / gamma
970            ynorm2 = z*z  +  ynorm2
971            rhs1   = rhs2 -  delta*z
972            rhs2   =      -  epsln*z
973    
974            # Estimate various norms and test for convergence.
975    
976            Anorm  = math.sqrt( tnorm2 )
977            ynorm  = math.sqrt( ynorm2 )
978    
979            rnorm  = phibar
980    
981        return x
982        
983    def NewtonGMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium,x=None, iter_max=100,iter_restart=20):
984    
985        gamma=.9
986        lmaxit=40
987        etamax=.5
988    
989        n = 1 #len(x)
990        iter=0
991        
992        # evaluate f at the initial iterate
993        # compute the stop tolerance
994        #
995        r=b
996        if x==None:
997          x=0*b
998        else:
999          r += (-1)*Aprod(x)
1000    
1001        f0=-Msolve(r)
1002        fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
1003        fnrmo=1
1004        atol=1.e-2
1005        rtol=1.e-4
1006        stop_tol=atol + rtol*fnrm
1007        #
1008        # main iteration loop
1009        #
1010        while not stoppingcriterium(fnrm,stop_tol,'NewtonGMRES',TOL=1.):
1011    
1012                if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1013            #
1014            # keep track of the ratio (rat = fnrm/frnmo)
1015            # of successive residual norms and
1016            # the iteration counter (iter)
1017            #
1018            #rat=fnrm/fnrmo
1019            fnrmo=fnrm
1020            iter+=1
1021            #
1022                # compute the step using a GMRES(m) routine especially designed for this purpose
1023            #
1024                initer=0
1025                while True:
1026                   xc,stopped=FDGMRES(f0, Aprod, Msolve, bilinearform, stoppingcriterium, xc=x, iter_max=lmaxit-initer, iter_restart=iter_restart, TOL=etamax)
1027                   if stopped: break
1028                   initer+=iter_restart
1029            xold=x
1030            x+=xc
1031            f0=-Msolve(Aprod(x))
1032            fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
1033            rat=fnrm/fnrmo
1034    
1035    
1036        #   adjust eta
1037        #
1038            if etamax > 0:
1039                etaold=etamax
1040                etanew=gamma*rat*rat
1041                if gamma*etaold*etaold > .1 :
1042                    etanew=max(etanew,gamma*etaold*etaold)
1043                
1044                etamax=min(etanew,etamax)
1045                etamax=max(etamax,.5*stop_tol/fnrm)
1046    
1047        return x
1048    
1049    def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
1050    
1051    # TFQMR solver for linear systems
1052    #
1053    #
1054    # initialization
1055    #
1056      errtol = math.sqrt(bilinearform(b,b))
1057      norm_b=errtol
1058      kmax  = iter_max
1059      error = []
1060    
1061      if math.sqrt(bilinearform(x,x)) != 0.0:
1062        r = b - Aprod(x)
1063      else:
1064        r = b
1065    
1066      r=Msolve(r)
1067    
1068      u1=0
1069      u2=0
1070      y1=0
1071      y2=0
1072    
1073      w = r
1074      y1 = r
1075      iter = 0
1076      d = 0
1077      
1078      v = Msolve(Aprod(y1))
1079      u1 = v
1080      
1081      theta = 0.0;
1082      eta = 0.0;
1083      tau = math.sqrt(bilinearform(r,r))
1084      error = [ error, tau ]
1085      rho = tau * tau
1086      m=1
1087    #
1088    #  TFQMR iteration
1089    #
1090    #  while ( iter < kmax-1 ):
1091      
1092      while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):
1093        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
1094    
1095        sigma = bilinearform(r,v)
1096    
1097        if ( sigma == 0.0 ):
1098          raise 'TFQMR breakdown, sigma=0'
1099        
1100    
1101        alpha = rho / sigma
1102    
1103        for j in range(2):
1104    #
1105    #   Compute y2 and u2 only if you have to
1106    #
1107          if ( j == 1 ):
1108            y2 = y1 - alpha * v
1109            u2 = Msolve(Aprod(y2))
1110    
1111          m = 2 * (iter+1) - 2 + (j+1)
1112          if j==0:
1113             w = w - alpha * u1
1114             d = y1 + ( theta * theta * eta / alpha ) * d
1115          if j==1:
1116             w = w - alpha * u2
1117             d = y2 + ( theta * theta * eta / alpha ) * d
1118    
1119          theta = math.sqrt(bilinearform(w,w))/ tau
1120          c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1121          tau = tau * theta * c
1122          eta = c * c * alpha
1123          x = x + eta * d
1124    #
1125    #  Try to terminate the iteration at each pass through the loop
1126    #
1127         # if ( tau * math.sqrt ( m + 1 ) <= errtol ):
1128         #   error = [ error, tau ]
1129         #   total_iters = iter
1130         #   break
1131          
1132    
1133        if ( rho == 0.0 ):
1134          raise 'TFQMR breakdown, rho=0'
1135        
1136    
1137        rhon = bilinearform(r,w)
1138        beta = rhon / rho;
1139        rho = rhon;
1140        y1 = w + beta * y2;
1141        u1 = Msolve(Aprod(y1))
1142        v = u1 + beta * ( u2 + beta * v )
1143        error = [ error, tau ]
1144        total_iters = iter
1145        
1146        iter = iter + 1
1147    
1148      return x
1149    
1150    
1151    #############################################
1152    
1153    class ArithmeticTuple(object):
1154       """
1155       tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
1156    
1157       example of usage:
1158    
1159       from esys.escript import Data
1160       from numarray import array
1161       a=Data(...)
1162       b=array([1.,4.])
1163       x=ArithmeticTuple(a,b)
1164       y=5.*x
1165    
1166       """
1167       def __init__(self,*args):
1168           """
1169           initialize object with elements args.
1170    
1171           @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
1172           """
1173           self.__items=list(args)
1174    
1175       def __len__(self):
1176           """
1177           number of items
1178    
1179           @return: number of items
1180           @rtype: C{int}
1181           """
1182           return len(self.__items)
1183    
1184       def __getitem__(self,index):
1185           """
1186           get an item
1187    
1188           @param index: item to be returned
1189           @type index: C{int}
1190           @return: item with index C{index}
1191           """
1192           return self.__items.__getitem__(index)
1193    
1194       def __mul__(self,other):
1195           """
1196           scaling from the right
1197    
1198           @param other: scaling factor
1199           @type other: C{float}
1200           @return: itemwise self*other
1201           @rtype: L{ArithmeticTuple}
1202           """
1203           out=[]
1204           for i in range(len(self)):
1205               out.append(self[i]*other)
1206           return ArithmeticTuple(*tuple(out))
1207    
1208       def __rmul__(self,other):
1209           """
1210           scaling from the left
1211    
1212           @param other: scaling factor
1213           @type other: C{float}
1214           @return: itemwise other*self
1215           @rtype: L{ArithmeticTuple}
1216           """
1217           out=[]
1218           for i in range(len(self)):
1219               out.append(other*self[i])
1220           return ArithmeticTuple(*tuple(out))
1221    
1222    #########################
1223    # Added by Artak
1224    #########################
1225       def __div__(self,other):
1226           """
1227           dividing from the right
1228    
1229           @param other: scaling factor
1230           @type other: C{float}
1231           @return: itemwise self/other
1232           @rtype: L{ArithmeticTuple}
1233           """
1234           out=[]
1235           for i in range(len(self)):
1236               out.append(self[i]/other)
1237           return ArithmeticTuple(*tuple(out))
1238    
1239       def __rdiv__(self,other):
1240           """
1241           dividing from the left
1242    
1243           @param other: scaling factor
1244           @type other: C{float}
1245           @return: itemwise other/self
1246           @rtype: L{ArithmeticTuple}
1247           """
1248           out=[]
1249           for i in range(len(self)):
1250               out.append(other/self[i])
1251           return ArithmeticTuple(*tuple(out))
1252      
1253    ##########################################33
1254    
1255       def __iadd__(self,other):
1256           """
1257           in-place add of other to self
1258    
1259           @param other: increment
1260           @type other: C{ArithmeticTuple}
1261           """
1262           if len(self) != len(other):
1263               raise ValueError,"tuple length must match."
1264           for i in range(len(self)):
1265               self.__items[i]+=other[i]
1266           return self
1267    
1268    class HomogeneousSaddlePointProblem(object):
1269          """
1270          This provides a framwork for solving homogeneous saddle point problem of the form
1271    
1272                 Av+B^*p=f
1273                 Bv    =0
1274    
1275          for the unknowns v and p and given operators A and B and given right hand side f.
1276          B^* is the adjoint operator of B is the given inner product.
1277    
1278          """
1279          def __init__(self,**kwargs):
1280            self.setTolerance()
1281            self.setToleranceReductionFactor()
1282    
1283          def initialize(self):
1284            """
1285            initialize the problem (overwrite)
1286            """
1287            pass
1288          def B(self,v):
1289             """
1290             returns Bv (overwrite)
1291             @rtype: equal to the type of p
1292    
1293             @note: boundary conditions on p should be zero!
1294             """
1295             pass
1296    
1297          def inner(self,p0,p1):
1298             """
1299             returns inner product of two element p0 and p1  (overwrite)
1300            
1301             @type p0: equal to the type of p
1302             @type p1: equal to the type of p
1303             @rtype: C{float}
1304    
1305             @rtype: equal to the type of p
1306             """
1307             pass
1308    
1309          def solve_A(self,u,p):
1310             """
1311             solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
1312    
1313             @rtype: equal to the type of v
1314             @note: boundary conditions on v should be zero!
1315             """
1316             pass
1317    
1318          def solve_prec(self,p):
1319             """
1320             provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
1321    
1322             @rtype: equal to the type of p
1323             """
1324             pass
1325    
1326          def stoppingcriterium(self,Bv,v,p):
1327             """
1328             returns a True if iteration is terminated. (overwrite)
1329    
1330             @rtype: C{bool}
1331             """
1332             pass
1333                
1334          def __inner(self,p,r):
1335             return self.inner(p,r[1])
1336    
1337          def __inner_p(self,p1,p2):
1338             return self.inner(p1,p2)
1339    
1340          def __stoppingcriterium(self,norm_r,r,p):
1341              return self.stoppingcriterium(r[1],r[0],p)
1342    
1343          def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):
1344              return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)
1345    
1346          def setTolerance(self,tolerance=1.e-8):
1347                  self.__tol=tolerance
1348          def getTolerance(self):
1349                  return self.__tol
1350          def setToleranceReductionFactor(self,reduction=0.01):
1351                  self.__reduction=reduction
1352          def getSubProblemTolerance(self):
1353                  return self.__reduction*self.getTolerance()
1354    
1355          def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):
1356                  """
1357                  solves the saddle point problem using initial guesses v and p.
1358    
1359                  @param max_iter: maximum number of iteration steps.
1360                  """
1361                  self.verbose=verbose
1362                  self.show_details=show_details and self.verbose
1363    
1364                  # assume p is known: then v=A^-1(f-B^*p)
1365                  # which leads to BA^-1B^*p = BA^-1f  
1366    
1367              # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)      
1368    
1369              
1370              self.__z=v+self.solve_A(v,p*0)
1371    
1372                  Bz=self.B(self.__z)
1373                  #
1374              #   solve BA^-1B^*p = Bz
1375                  #
1376                  #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1377                  #
1378                  #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1379                  #                           A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)
1380                  #
1381                  self.iter=0
1382              if solver=='GMRES':      
1383                    if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1384                    p=GMRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)
1385                    # solve Au=f-B^*p
1386                    #       A(u-v)=f-B^*p-Av
1387                    #       u=v+(u-v)
1388            u=v+self.solve_A(v,p)
1389    
1390              if solver=='NewtonGMRES':    
1391                    if self.verbose: print "enter NewtonGMRES method (iter_max=%s)"%max_iter
1392                    p=NewtonGMRES(Bz,self.__Aprod_Newton,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1393                    # solve Au=f-B^*p
1394                    #       A(u-v)=f-B^*p-Av
1395                    #       u=v+(u-v)
1396            u=v+self.solve_A(v,p)
1397                    
1398    
1399              if solver=='TFQMR':      
1400                    if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter
1401                    p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1402                    # solve Au=f-B^*p
1403                    #       A(u-v)=f-B^*p-Av
1404                    #       u=v+(u-v)
1405            u=v+self.solve_A(v,p)
1406    
1407              if solver=='MINRES':      
1408                    if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1409                    p=MINRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
1410                    # solve Au=f-B^*p
1411                    #       A(u-v)=f-B^*p-Av
1412                    #       u=v+(u-v)
1413            u=v+self.solve_A(v,p)
1414                  
1415                  if solver=='PCG':
1416                    if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1417                    p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1418                u=r[0]  
1419    
1420                  print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1421    
1422              return u,p
1423    
1424          def __Msolve(self,r):
1425              return self.solve_prec(r[1])
1426    
1427          def __Msolve2(self,r):
1428              return self.solve_prec(r)
1429    
1430    
1431          def __Aprod(self,p):
1432              # return BA^-1B*p
1433              #solve Av =-B^*p as Av =f-Az-B^*p
1434              v=self.solve_A(self.__z,-p)
1435              return ArithmeticTuple(v, self.B(v))
1436    
1437          def __Aprod2(self,p):
1438              # return BA^-1B*p
1439              #solve Av =-B^*p as Av =f-Az-B^*p
1440          v=self.solve_A(self.__z,-p)
1441              return self.B(v)
1442    
1443          def __Aprod_Newton(self,p):
1444              # return BA^-1B*p
1445              #solve Av =-B^*p as Av =f-Az-B^*p
1446          v=self.solve_A(self.__z,-p)
1447              return self.B(v-self.__z)
1448    
1449  class SaddlePointProblem(object):  class SaddlePointProblem(object):
1450     """     """
1451     This implements a solver for a saddlepoint problem     This implements a solver for a saddlepoint problem
# Line 420  class SaddlePointProblem(object): Line 1471  class SaddlePointProblem(object):
1471         @type verbose: C{bool}         @type verbose: C{bool}
1472         @note: this method may be overwritten by a particular saddle point problem         @note: this method may be overwritten by a particular saddle point problem
1473         """         """
1474         if verbose:         if not isinstance(verbose,bool):
1475            self.__verbose=True              raise TypeError("verbose needs to be of type bool.")
1476         else:         self.__verbose=verbose
           self.__verbose=False  
1477         self.relaxation=1.         self.relaxation=1.
1478    
1479     def trace(self,text):     def trace(self,text):
# Line 433  class SaddlePointProblem(object): Line 1483  class SaddlePointProblem(object):
1483         @param text: a text message         @param text: a text message
1484         @type text: C{str}         @type text: C{str}
1485         """         """
1486         if self.__verbose: print "%s: %s"%(str(self),text)         if self.__verbose: print "#s: #s"%(str(self),text)
1487    
1488     def solve_f(self,u,p,tol=1.e-8):     def solve_f(self,u,p,tol=1.e-8):
1489         """         """
# Line 583  class SaddlePointProblem(object): Line 1633  class SaddlePointProblem(object):
1633              norm_u_new = util.Lsup(u_new)              norm_u_new = util.Lsup(u_new)
1634              p_new=p+self.relaxation*g_new              p_new=p+self.relaxation*g_new
1635              norm_p_new = util.sqrt(self.inner(p_new,p_new))              norm_p_new = util.sqrt(self.inner(p_new,p_new))
1636              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))
1637    
1638              if self.iter>1:              if self.iter>1:
1639                 dg2=g_new-g                 dg2=g_new-g
# Line 625  class SaddlePointProblem(object): Line 1675  class SaddlePointProblem(object):
1675  #  #
1676  #      return u,p  #      return u,p
1677                        
1678  # vim: expandtab shiftwidth=4:  def MaskFromBoundaryTag(function_space,*tags):
1679       """
1680       create a mask on the given function space which one for samples
1681       that touch the boundary tagged by tags.
1682    
1683       usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
1684    
1685       @param function_space: a given function space
1686       @type function_space: L{escript.FunctionSpace}
1687       @param tags: boundray tags
1688       @type tags: C{str}
1689       @return: a mask which marks samples used by C{function_space} that are touching the
1690                boundary tagged by any of the given tags.
1691       @rtype: L{escript.Data} of rank 0
1692       """
1693       pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
1694       d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
1695       for t in tags: d.setTaggedValue(t,1.)
1696       pde.setValue(y=d)
1697       out=util.whereNonZero(pde.getRightHandSide())
1698       if out.getFunctionSpace() == function_space:
1699          return out
1700       else:
1701          return util.whereNonZero(util.interpolate(out,function_space))
1702    
1703    
1704    

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

  ViewVC Help
Powered by ViewVC 1.1.26