/[escript]/trunk/escript/py_src/pdetools.py
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1312 by ksteube, Mon Sep 24 06:18:44 2007 UTC revision 1489 by artak, Mon Apr 14 04:29:30 2008 UTC
# Line 48  import numarray Line 48  import numarray
48  import util  import util
49  import math  import math
50    
51    ##### Added by Artak
52    # from Numeric import zeros,Int,Float64
53    ###################################
54    
55    
56  class TimeIntegrationManager:  class TimeIntegrationManager:
57    """    """
58    a simple mechanism to manage time dependend values.    a simple mechanism to manage time dependend values.
# Line 435  class NegativeNorm(SolverSchemeException Line 440  class NegativeNorm(SolverSchemeException
440     """     """
441     pass     pass
442    
443  def PCG(b,x,Aprod,Msolve,bilinearform, norm, verbose=True, iter_max=100, tolerance=math.sqrt(util.EPSILON)):  class IterationHistory(object):
444     """     """
445     Solver for     The IterationHistory class is used to define a stopping criterium. It keeps track of the
446       residual norms. The stoppingcriterium indicates termination if the residual norm has been reduced by
447       a given tolerance.
448       """
449       def __init__(self,tolerance=math.sqrt(util.EPSILON),verbose=False):
450          """
451          Initialization
452    
453     M{Ax=b}        @param tolerance: tolerance
454          @type tolerance: positive C{float}
455          @param verbose: switches on the printing out some information
456          @type verbose: C{bool}
457          """
458          if not tolerance>0.:
459              raise ValueError,"tolerance needs to be positive."
460          self.tolerance=tolerance
461          self.verbose=verbose
462          self.history=[]
463       def stoppingcriterium(self,norm_r,r,x):
464           """
465           returns True if the C{norm_r} is C{tolerance}*C{norm_r[0]} where C{norm_r[0]}  is the residual norm at the first call.
466    
467     with a symmetric and positive definite operator A (more details required!).        
468     It uses the conjugate gradient method with preconditioner M providing an approximation of A.         @param norm_r: current residual norm
469           @type norm_r: non-negative C{float}
470           @param r: current residual (not used)
471           @param x: current solution approximation (not used)
472           @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
473           @rtype: C{bool}
474    
475     The iteration is terminated if         """
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     M{norm(r) <= tolerance * norm(b)}     def stoppingcriterium2(self,norm_r,norm_b):
481           """
482           returns True if the C{norm_r} is C{tolerance}*C{norm_b}
483    
484     where C{norm()} defines a norm and        
485           @param norm_r: current residual norm
486           @type norm_r: non-negative C{float}
487           @param norm_b: norm of right hand side
488           @type norm_b: non-negative C{float}
489           @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
490           @rtype: C{bool}
491    
492     M{r = b-Ax}         """
493           self.history.append(norm_r)
494           if self.verbose: print "iter: #s:  norm(r) = #e"#(len(self.history)-1, self.history[-1])
495           return self.history[-1]<=self.tolerance * norm_b
496    
497     the residual.  def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
498       """
499       Solver for
500    
501       M{Ax=b}
502    
503       with a symmetric and positive definite operator A (more details required!).
504       It uses the conjugate gradient method with preconditioner M providing an approximation of A.
505    
506       The iteration is terminated if the C{stoppingcriterium} function return C{True}.
507    
508     For details on the preconditioned conjugate gradient method see the book:     For details on the preconditioned conjugate gradient method see the book:
509    
# Line 461  def PCG(b,x,Aprod,Msolve,bilinearform, n Line 512  def PCG(b,x,Aprod,Msolve,bilinearform, n
512     C. Romine, and H. van der Vorst.     C. Romine, and H. van der Vorst.
513    
514     @param b: the right hand side of the liner system. C{b} is altered.     @param b: the right hand side of the liner system. C{b} is altered.
515     @type b: any object type R supporting inplace add (x+=y) and scaling (x=scalar*y)     @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
    @param x: an initial guess for the solution  
    @type x: any object type S supporting inplace add (x+=y), scaling (x=scalar*y)  
516     @param Aprod: returns the value Ax     @param Aprod: returns the value Ax
517     @type Aprod: function C{Aprod(x)} where C{x} is of object type S. The returned object needs to be of type R.     @type Aprod: function C{Aprod(x)} where C{x} is of the same object like argument C{x}. The returned object needs to be of the same type like argument C{b}.
518     @param Msolve: solves Mx=r     @param Msolve: solves Mx=r
519     @type Msolve: function C{Msolve(r)} where C{r} is of object type R. The returned object needs to be of tupe S.     @type Msolve: function C{Msolve(r)} where C{r} is of the same type like argument C{b}. The returned object needs to be of the same
520    type like argument C{x}.
521     @param bilinearform: inner product C{<x,r>}     @param bilinearform: inner product C{<x,r>}
522     @type bilinearform: function C{bilinearform(x,r)} where C{x} is of object type S and C{r} is of object type R. The returned value is a C{float}.     @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same type like argument C{x} and C{r} is . The returned value is a C{float}.
523     @param norm: norm calculation for the residual C{r=b-Ax}.     @param stoppingcriterium: function which returns True if a stopping criterium is meet. C{stoppingcriterium} has the arguments C{norm_r}, C{r} and C{x} giving the current norm of the residual (=C{sqrt(bilinearform(Msolve(r),r)}), the current residual and the current solution approximation. C{stoppingcriterium} is called in each iteration step.
524     @type norm: function C{norm(r)} where C{r} is of object type R. The returned value is a C{float}.     @type stoppingcriterium: function that returns C{True} or C{False}
525     @param verbose: switches on the printing out some information     @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
526     @type verbose: C{bool}     @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
527     @param iter_max: maximum number of iteration steps.     @param iter_max: maximum number of iteration steps.
528     @type iter_max: C{int}     @type iter_max: C{int}
529     @param tolerance: tolerance     @return: the solution approximation and the corresponding residual
530     @type tolerance: positive C{float}     @rtype: C{tuple}
531     @return: the solution apprximation and the corresponding residual     @warning: C{b} and C{x} are altered.
    @rtype: C{tuple} of an S type and and an R type object.A  
    @warning: C{b} ans C{x} are altered.  
532     """     """
    if verbose:  
         print "Enter PCG for solving Ax=b\n\t iter_max =%s\t tolerance   =%e"%(iter_max,tolerance)  
533     iter=0     iter=0
534     normb = norm(b)     if x==None:
535     if normb<0: raise NegativeNorm        x=0*b
536       else:
537     b += (-1)*Aprod(x)        b += (-1)*Aprod(x)
538     r=b     r=b
539     rhat=Msolve(r)     rhat=Msolve(r)
540     d = rhat;     d = rhat
541     rhat_dot_r = bilinearform(rhat, r)     rhat_dot_r = bilinearform(rhat, r)
542       if rhat_dot_r<0: raise NegativeNorm,"negative norm."
543    
544     while True:     while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
545         normr=norm(r)         iter+=1
        if normr<0: raise NegativeNorm  
        if verbose: print "iter: %s: norm(r) = %e, tolerance*norm(b) = %e"%(iter, normr,tolerance * normb)  
        if normr <= tolerance * normb: return x,r  
   
        iter+=1 # k = iter = 1 first time through  
546         if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max         if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
547    
548         q=Aprod(d)         q=Aprod(d)
# Line 515  def PCG(b,x,Aprod,Msolve,bilinearform, n Line 557  def PCG(b,x,Aprod,Msolve,bilinearform, n
557         d=rhat         d=rhat
558    
559         rhat_dot_r = rhat_dot_r_new         rhat_dot_r = rhat_dot_r_new
560           if rhat_dot_r<0: raise NegativeNorm,"negative norm."
561    
562       return x,r
563    
564    
565    ############################
566    # Added by Artak
567    #################################3
568    
569    #Apply a sequence of k Givens rotations, used within gmres codes
570    # vrot=givapp(c, s, vin, k)
571    def givapp(c,s,vin):
572        vrot=vin # warning: vin is altered!!!!
573        if isinstance(c,float):
574            vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
575        else:
576            for i in range(len(c)):
577                w1=c[i]*vrot[i]-s[i]*vrot[i+1]
578            w2=s[i]*vrot[i]+c[i]*vrot[i+1]
579                vrot[i:i+2]=w1,w2
580        return vrot
581    
582    def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=10):
583       m=iter_restart
584       iter=0
585       while True:
586          if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
587          x,stopped=GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=x, iter_max=iter_max-iter, iter_restart=m)
588          iter+=iter_restart
589          if stopped: break
590       return x
591    
592    def GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=10):
593       iter=0
594       r=Msolve(b)
595       r_dot_r = bilinearform(r, r)
596       if r_dot_r<0: raise NegativeNorm,"negative norm."
597       norm_b=math.sqrt(r_dot_r)
598    
599       if x==None:
600          x=0
601       else:
602          r=Msolve(b-Aprod(x))
603          r_dot_r = bilinearform(r, r)
604          if r_dot_r<0: raise NegativeNorm,"negative norm."
605      
606       h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
607       c=numarray.zeros(iter_restart,numarray.Float64)
608       s=numarray.zeros(iter_restart,numarray.Float64)
609       g=numarray.zeros(iter_restart,numarray.Float64)
610       v=[]
611    
612       rho=math.sqrt(r_dot_r)
613      
614       v.append(r/rho)
615       g[0]=rho
616    
617       while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
618    
619        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
620    
621        
622        p=Msolve(Aprod(v[iter]))
623    
624        v.append(p)
625    
626        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
627    
628    # Modified Gram-Schmidt
629        for j in range(iter+1):
630          h[j][iter]=bilinearform(v[j],v[iter+1])  
631          v[iter+1]+=(-1.)*h[j][iter]*v[j]
632          
633        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
634        v_norm2=h[iter+1][iter]
635    
636    
637    # Reorthogonalize if needed
638        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
639         for j in range(iter+1):
640            hr=bilinearform(v[j],v[iter+1])
641                h[j][iter]=h[j][iter]+hr #vhat
642                v[iter+1] +=(-1.)*hr*v[j]
643    
644         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
645         h[iter+1][iter]=v_norm2
646    
647    #   watch out for happy breakdown
648            if v_norm2 != 0:
649             v[iter+1]=v[iter+1]/h[iter+1][iter]
650    
651    #   Form and store the information for the new Givens rotation
652        if iter > 0 :
653            hhat=[]
654            for i in range(iter+1) : hhat.append(h[i][iter])
655            hhat=givapp(c[0:iter],s[0:iter],hhat);
656                for i in range(iter+1) : h[i][iter]=hhat[i]
657    
658        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
659        if mu!=0 :
660            c[iter]=h[iter][iter]/mu
661            s[iter]=-h[iter+1][iter]/mu
662            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
663            h[iter+1][iter]=0.0
664            g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
665    
666    # Update the residual norm
667            rho=abs(g[iter+1])
668        iter+=1
669    
670    # At this point either iter > iter_max or rho < tol.
671    # It's time to compute x and leave.        
672    
673       if iter > 0 :
674         y=numarray.zeros(iter,numarray.Float64)    
675         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
676         if iter > 1 :  
677            i=iter-2  
678            while i>=0 :
679              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
680              i=i-1
681         xhat=v[iter-1]*y[iter-1]
682         for i in range(iter-1):
683        xhat += v[i]*y[i]
684       else : xhat=v[0]
685        
686       x += xhat
687       if iter<iter_restart-1:
688          stopped=True
689       else:
690          stopped=False
691    
692       return x,stopped
693    
694    def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
695    
696        #
697        #  minres solves the system of linear equations Ax = b
698        #  where A is a symmetric matrix (possibly indefinite or singular)
699        #  and b is a given vector.
700        #  
701        #  "A" may be a dense or sparse matrix (preferably sparse!)
702        #  or the name of a function such that
703        #               y = A(x)
704        #  returns the product y = Ax for any given vector x.
705        #
706        #  "M" defines a positive-definite preconditioner M = C C'.
707        #  "M" may be a dense or sparse matrix (preferably sparse!)
708        #  or the name of a function such that
709        #  solves the system My = x for any given vector x.
710        #
711        #
712        
713        #------------------------------------------------------------------
714        # Set up y and v for the first Lanczos vector v1.
715        # y  =  beta1 P' v1,  where  P = C**(-1).
716        # v is really P' v1.
717        #------------------------------------------------------------------
718        if x==None:
719          x=0*b
720        else:
721          b += (-1)*Aprod(x)
722    
723        r1    = b
724        y = Msolve(b)
725        beta1 = bilinearform(b,y)
726    
727        if beta1< 0: raise NegativeNorm,"negative norm."
728    
729        #  If b = 0 exactly, stop with x = 0.
730        if beta1==0: return x*0.
731    
732        if beta1> 0:
733          beta1  = math.sqrt(beta1)      
734    
735        #------------------------------------------------------------------
736        # Initialize quantities.
737        # ------------------------------------------------------------------
738        iter   = 0
739        Anorm = 0
740        ynorm = 0
741        oldb   = 0
742        beta   = beta1
743        dbar   = 0
744        epsln  = 0
745        phibar = beta1
746        rhs1   = beta1
747        rhs2   = 0
748        rnorm  = phibar
749        tnorm2 = 0
750        ynorm2 = 0
751        cs     = -1
752        sn     = 0
753        w      = b*0.
754        w2     = b*0.
755        r2     = r1
756        eps    = 0.0001
757    
758        #---------------------------------------------------------------------
759        # Main iteration loop.
760        # --------------------------------------------------------------------
761        while not stoppingcriterium(rnorm,Anorm*ynorm):    #  checks ||r|| < (||A|| ||x||) * TOL
762    
763        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
764            iter    = iter  +  1
765    
766            #-----------------------------------------------------------------
767            # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
768            # The general iteration is similar to the case k = 1 with v0 = 0:
769            #
770            #   p1      = Operator * v1  -  beta1 * v0,
771            #   alpha1  = v1'p1,
772            #   q2      = p2  -  alpha1 * v1,
773            #   beta2^2 = q2'q2,
774            #   v2      = (1/beta2) q2.
775            #
776            # Again, y = betak P vk,  where  P = C**(-1).
777            #-----------------------------------------------------------------
778            s = 1/beta                 # Normalize previous vector (in y).
779            v = s*y                    # v = vk if P = I
780        
781            y      = Aprod(v)
782        
783            if iter >= 2:
784              y = y - (beta/oldb)*r1
785    
786            alfa   = bilinearform(v,y)              # alphak
787            y      = (- alfa/beta)*r2 + y
788            r1     = r2
789            r2     = y
790            y = Msolve(r2)
791            oldb   = beta                           # oldb = betak
792            beta   = bilinearform(r2,y)             # beta = betak+1^2
793            if beta < 0: raise NegativeNorm,"negative norm."
794    
795            beta   = math.sqrt( beta )
796            tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
797            
798            if iter==1:                 # Initialize a few things.
799              gmax   = abs( alfa )      # alpha1
800              gmin   = gmax             # alpha1
801    
802            # Apply previous rotation Qk-1 to get
803            #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
804            #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
805        
806            oldeps = epsln
807            delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
808            gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
809            epsln  =               sn * beta  # epsln2 = 0         epslnk+1
810            dbar   =            -  cs * beta  # dbar 2 = beta2     dbar k+1
811    
812            # Compute the next plane rotation Qk
813    
814            gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
815            gamma  = max(gamma,eps)
816            cs     = gbar / gamma             # ck
817            sn     = beta / gamma             # sk
818            phi    = cs * phibar              # phik
819            phibar = sn * phibar              # phibark+1
820    
821            # Update  x.
822    
823            denom = 1/gamma
824            w1    = w2
825            w2    = w
826            w     = (v - oldeps*w1 - delta*w2) * denom
827            x     = x  +  phi*w
828    
829            # Go round again.
830    
831            gmax   = max(gmax,gamma)
832            gmin   = min(gmin,gamma)
833            z      = rhs1 / gamma
834            ynorm2 = z*z  +  ynorm2
835            rhs1   = rhs2 -  delta*z
836            rhs2   =      -  epsln*z
837    
838            # Estimate various norms and test for convergence.
839    
840            Anorm  = math.sqrt( tnorm2 )
841            ynorm  = math.sqrt( ynorm2 )
842    
843            rnorm  = phibar
844    
845        return x
846        
847    
848    def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
849    
850    # TFQMR solver for linear systems
851    #
852    #
853    # initialization
854    #
855      errtol = math.sqrt(bilinearform(b,b))
856      norm_b=errtol
857      kmax  = iter_max
858      error = []
859    
860      if math.sqrt(bilinearform(x,x)) != 0.0:
861        r = b - Aprod(x)
862      else:
863        r = b
864    
865      r=Msolve(r)
866    
867      u1=0
868      u2=0
869      y1=0
870      y2=0
871    
872      w = r
873      y1 = r
874      iter = 0
875      d = 0
876      
877      v = Msolve(Aprod(y1))
878      u1 = v
879      
880      theta = 0.0;
881      eta = 0.0;
882      tau = math.sqrt(bilinearform(r,r))
883      error = [ error, tau ]
884      rho = tau * tau
885      m=1
886    #
887    #  TFQMR iteration
888    #
889    #  while ( iter < kmax-1 ):
890      
891      while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b):
892        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
893    
894        sigma = bilinearform(r,v)
895    
896        if ( sigma == 0.0 ):
897          raise 'TFQMR breakdown, sigma=0'
898        
899    
900        alpha = rho / sigma
901    
902        for j in range(2):
903    #
904    #   Compute y2 and u2 only if you have to
905    #
906          if ( j == 1 ):
907            y2 = y1 - alpha * v
908            u2 = Msolve(Aprod(y2))
909    
910          m = 2 * (iter+1) - 2 + (j+1)
911          if j==0:
912             w = w - alpha * u1
913             d = y1 + ( theta * theta * eta / alpha ) * d
914          if j==1:
915             w = w - alpha * u2
916             d = y2 + ( theta * theta * eta / alpha ) * d
917    
918          theta = math.sqrt(bilinearform(w,w))/ tau
919          c = 1.0 / math.sqrt ( 1.0 + theta * theta )
920          tau = tau * theta * c
921          eta = c * c * alpha
922          x = x + eta * d
923    #
924    #  Try to terminate the iteration at each pass through the loop
925    #
926         # if ( tau * math.sqrt ( m + 1 ) <= errtol ):
927         #   error = [ error, tau ]
928         #   total_iters = iter
929         #   break
930          
931    
932        if ( rho == 0.0 ):
933          raise 'TFQMR breakdown, rho=0'
934        
935    
936        rhon = bilinearform(r,w)
937        beta = rhon / rho;
938        rho = rhon;
939        y1 = w + beta * y2;
940        u1 = Msolve(Aprod(y1))
941        v = u1 + beta * ( u2 + beta * v )
942        error = [ error, tau ]
943        total_iters = iter
944        
945        iter = iter + 1
946    
947      print iter
948      return x
949    
950    
951    #############################################
952    
953    class ArithmeticTuple(object):
954       """
955       tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
956    
957       example of usage:
958    
959       from esys.escript import Data
960       from numarray import array
961       a=Data(...)
962       b=array([1.,4.])
963       x=ArithmeticTuple(a,b)
964       y=5.*x
965    
966       """
967       def __init__(self,*args):
968           """
969           initialize object with elements args.
970    
971           @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
972           """
973           self.__items=list(args)
974    
975       def __len__(self):
976           """
977           number of items
978    
979           @return: number of items
980           @rtype: C{int}
981           """
982           return len(self.__items)
983    
984       def __getitem__(self,index):
985           """
986           get an item
987    
988           @param index: item to be returned
989           @type index: C{int}
990           @return: item with index C{index}
991           """
992           return self.__items.__getitem__(index)
993    
994       def __mul__(self,other):
995           """
996           scaling from the right
997    
998           @param other: scaling factor
999           @type other: C{float}
1000           @return: itemwise self*other
1001           @rtype: L{ArithmeticTuple}
1002           """
1003           out=[]
1004           for i in range(len(self)):
1005               out.append(self[i]*other)
1006           return ArithmeticTuple(*tuple(out))
1007    
1008       def __rmul__(self,other):
1009           """
1010           scaling from the left
1011    
1012           @param other: scaling factor
1013           @type other: C{float}
1014           @return: itemwise other*self
1015           @rtype: L{ArithmeticTuple}
1016           """
1017           out=[]
1018           for i in range(len(self)):
1019               out.append(other*self[i])
1020           return ArithmeticTuple(*tuple(out))
1021    
1022    #########################
1023    # Added by Artak
1024    #########################
1025       def __div__(self,other):
1026           """
1027           dividing from the right
1028    
1029           @param other: scaling factor
1030           @type other: C{float}
1031           @return: itemwise self/other
1032           @rtype: L{ArithmeticTuple}
1033           """
1034           out=[]
1035           for i in range(len(self)):
1036               out.append(self[i]/other)
1037           return ArithmeticTuple(*tuple(out))
1038    
1039       def __rdiv__(self,other):
1040           """
1041           dividing from the left
1042    
1043           @param other: scaling factor
1044           @type other: C{float}
1045           @return: itemwise other/self
1046           @rtype: L{ArithmeticTuple}
1047           """
1048           out=[]
1049           for i in range(len(self)):
1050               out.append(other/self[i])
1051           return ArithmeticTuple(*tuple(out))
1052      
1053    ##########################################33
1054    
1055       def __iadd__(self,other):
1056           """
1057           in-place add of other to self
1058    
1059           @param other: increment
1060           @type other: C{ArithmeticTuple}
1061           """
1062           if len(self) != len(other):
1063               raise ValueError,"tuple length must match."
1064           for i in range(len(self)):
1065               self.__items[i]+=other[i]
1066           return self
1067    
1068    class HomogeneousSaddlePointProblem(object):
1069          """
1070          This provides a framwork for solving homogeneous saddle point problem of the form
1071    
1072                 Av+B^*p=f
1073                 Bv    =0
1074    
1075          for the unknowns v and p and given operators A and B and given right hand side f.
1076          B^* is the adjoint operator of B is the given inner product.
1077    
1078          """
1079          def __init__(self,**kwargs):
1080            self.setTolerance()
1081            self.setToleranceReductionFactor()
1082    
1083          def initialize(self):
1084            """
1085            initialize the problem (overwrite)
1086            """
1087            pass
1088          def B(self,v):
1089             """
1090             returns Bv (overwrite)
1091             @rtype: equal to the type of p
1092    
1093             @note: boundary conditions on p should be zero!
1094             """
1095             pass
1096    
1097          def inner(self,p0,p1):
1098             """
1099             returns inner product of two element p0 and p1  (overwrite)
1100            
1101             @type p0: equal to the type of p
1102             @type p1: equal to the type of p
1103             @rtype: C{float}
1104    
1105             @rtype: equal to the type of p
1106             """
1107             pass
1108    
1109          def solve_A(self,u,p):
1110             """
1111             solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
1112    
1113             @rtype: equal to the type of v
1114             @note: boundary conditions on v should be zero!
1115             """
1116             pass
1117    
1118          def solve_prec(self,p):
1119             """
1120             provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
1121    
1122             @rtype: equal to the type of p
1123             """
1124             pass
1125    
1126          def stoppingcriterium(self,Bv,v,p):
1127             """
1128             returns a True if iteration is terminated. (overwrite)
1129    
1130             @rtype: C{bool}
1131             """
1132             pass
1133                
1134          def __inner(self,p,r):
1135             return self.inner(p,r[1])
1136    
1137          def __inner_p(self,p1,p2):
1138             return self.inner(p1,p2)
1139    
1140          def __stoppingcriterium(self,norm_r,r,p):
1141              return self.stoppingcriterium(r[1],r[0],p)
1142    
1143          def __stoppingcriterium_GMRES(self,norm_r,norm_b):
1144              return self.stoppingcriterium_GMRES(norm_r,norm_b)
1145    
1146          def __stoppingcriterium_MINRES(self,norm_r,norm_Ax):
1147              return self.stoppingcriterium_MINRES(norm_r,norm_Ax)
1148    
1149    
1150          def setTolerance(self,tolerance=1.e-8):
1151                  self.__tol=tolerance
1152          def getTolerance(self):
1153                  return self.__tol
1154          def setToleranceReductionFactor(self,reduction=0.01):
1155                  self.__reduction=reduction
1156          def getSubProblemTolerance(self):
1157                  return self.__reduction*self.getTolerance()
1158    
1159          def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=10):
1160                  """
1161                  solves the saddle point problem using initial guesses v and p.
1162    
1163                  @param max_iter: maximum number of iteration steps.
1164                  """
1165                  self.verbose=verbose
1166                  self.show_details=show_details and self.verbose
1167    
1168                  # assume p is known: then v=A^-1(f-B^*p)
1169                  # which leads to BA^-1B^*p = BA^-1f  
1170    
1171              # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)      
1172    
1173              
1174              self.__z=v+self.solve_A(v,p*0)
1175    
1176                  Bz=self.B(self.__z)
1177                  #
1178              #   solve BA^-1B^*p = Bz
1179                  #
1180                  #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1181                  #
1182                  #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1183                  #                           A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)
1184                  #
1185                  self.iter=0
1186              if solver=='GMRES':      
1187                    if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1188                    p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)
1189                    # solve Au=f-B^*p
1190                    #       A(u-v)=f-B^*p-Av
1191                    #       u=v+(u-v)
1192            u=v+self.solve_A(v,p)
1193    
1194              if solver=='TFQMR':      
1195                    if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1196                    p=TFQMR(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1.)
1197                    # solve Au=f-B^*p
1198                    #       A(u-v)=f-B^*p-Av
1199                    #       u=v+(u-v)
1200            u=v+self.solve_A(v,p)
1201    
1202              if solver=='MINRES':      
1203                    if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1204                    p=MINRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_MINRES,iter_max=max_iter, x=p*1.)
1205                    # solve Au=f-B^*p
1206                    #       A(u-v)=f-B^*p-Av
1207                    #       u=v+(u-v)
1208            u=v+self.solve_A(v,p)
1209                  
1210                  if solver=='PCG':
1211                    if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1212                    p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1213                u=r[0]  
1214    
1215                  print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1216    
1217              return u,p
1218    
1219          def __Msolve(self,r):
1220              return self.solve_prec(r[1])
1221    
1222          def __Msolve_GMRES(self,r):
1223              return self.solve_prec(r)
1224    
1225    
1226          def __Aprod(self,p):
1227              # return BA^-1B*p
1228              #solve Av =-B^*p as Av =f-Az-B^*p
1229              v=self.solve_A(self.__z,-p)
1230              return ArithmeticTuple(v, self.B(v))
1231    
1232          def __Aprod_GMRES(self,p):
1233              # return BA^-1B*p
1234              #solve Av =-B^*p as Av =f-Az-B^*p
1235          v=self.solve_A(self.__z,-p)
1236              return self.B(v)
1237    
1238  class SaddlePointProblem(object):  class SaddlePointProblem(object):
1239     """     """
# Line 554  class SaddlePointProblem(object): Line 1272  class SaddlePointProblem(object):
1272         @param text: a text message         @param text: a text message
1273         @type text: C{str}         @type text: C{str}
1274         """         """
1275         if self.__verbose: print "%s: %s"%(str(self),text)         if self.__verbose: print "#s: #s"%(str(self),text)
1276    
1277     def solve_f(self,u,p,tol=1.e-8):     def solve_f(self,u,p,tol=1.e-8):
1278         """         """
# Line 771  def MaskFromBoundaryTag(function_space,* Line 1489  def MaskFromBoundaryTag(function_space,*
1489     else:     else:
1490        return util.whereNonZero(util.interpolate(out,function_space))        return util.whereNonZero(util.interpolate(out,function_space))
1491    
1492    
1493    

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

  ViewVC Help
Powered by ViewVC 1.1.26