/[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 1122 by gross, Tue May 1 03:21:04 2007 UTC revision 1514 by artak, Wed Apr 16 00:15:44 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 393  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):
481           """
482           returns True if the C{norm_r} is C{tolerance}*C{norm_b}
483    
484          
485           @param norm_r: current residual norm
486           @type norm_r: non-negative C{float}
487           @param norm_b: norm of right hand side
488           @type norm_b: non-negative C{float}
489           @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.
490           @rtype: C{bool}
491    
492           """
493           self.history.append(norm_r)
494           if self.verbose: print "iter: #s:  norm(r) = #e"#(len(self.history)-1, self.history[-1])
495           return self.history[-1]<=self.tolerance * norm_b
496    
497    def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
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:
509    
510       Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
511       T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
512       C. Romine, and H. van der Vorst.
513    
514       @param b: the right hand side of the liner system. C{b} is altered.
515       @type b: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
516       @param Aprod: returns the value Ax
517       @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
519       @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>}
522       @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 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 stoppingcriterium: function that returns C{True} or C{False}
525       @param x: an initial guess for the solution. If no C{x} is given 0*b is used.
526       @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
527       @param iter_max: maximum number of iteration steps.
528       @type iter_max: C{int}
529       @return: the solution approximation and the corresponding residual
530       @rtype: C{tuple}
531       @warning: C{b} and C{x} are altered.
532       """
533       iter=0
534       if x==None:
535          x=0*b
536       else:
537          b += (-1)*Aprod(x)
538       r=b
539       rhat=Msolve(r)
540       d = rhat
541       rhat_dot_r = bilinearform(rhat, r)
542       if rhat_dot_r<0: raise NegativeNorm,"negative norm."
543    
544       while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):
545           iter+=1
546           if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
547    
548           q=Aprod(d)
549           alpha = rhat_dot_r / bilinearform(d, q)
550           x += alpha * d
551           r += (-alpha) * q
552    
553           rhat=Msolve(r)
554           rhat_dot_r_new = bilinearform(rhat, r)
555           beta = rhat_dot_r_new / rhat_dot_r
556           rhat+=beta * d
557           d=rhat
558    
559           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=20):
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=20):
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     """     """
1240     This implements a solver for a saddlepoint problem     This implements a solver for a saddlepoint problem
# Line 430  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 580  class SaddlePointProblem(object): Line 1422  class SaddlePointProblem(object):
1422              norm_u_new = util.Lsup(u_new)              norm_u_new = util.Lsup(u_new)
1423              p_new=p+self.relaxation*g_new              p_new=p+self.relaxation*g_new
1424              norm_p_new = util.sqrt(self.inner(p_new,p_new))              norm_p_new = util.sqrt(self.inner(p_new,p_new))
1425              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))
1426    
1427              if self.iter>1:              if self.iter>1:
1428                 dg2=g_new-g                 dg2=g_new-g
# Line 622  class SaddlePointProblem(object): Line 1464  class SaddlePointProblem(object):
1464  #  #
1465  #      return u,p  #      return u,p
1466                        
1467  # vim: expandtab shiftwidth=4:  def MaskFromBoundaryTag(function_space,*tags):
1468       """
1469       create a mask on the given function space which one for samples
1470       that touch the boundary tagged by tags.
1471    
1472       usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
1473    
1474       @param function_space: a given function space
1475       @type function_space: L{escript.FunctionSpace}
1476       @param tags: boundray tags
1477       @type tags: C{str}
1478       @return: a mask which marks samples used by C{function_space} that are touching the
1479                boundary tagged by any of the given tags.
1480       @rtype: L{escript.Data} of rank 0
1481       """
1482       pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
1483       d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
1484       for t in tags: d.setTaggedValue(t,1.)
1485       pde.setValue(y=d)
1486       out=util.whereNonZero(pde.getRightHandSide())
1487       if out.getFunctionSpace() == function_space:
1488          return out
1489       else:
1490          return util.whereNonZero(util.interpolate(out,function_space))
1491    
1492    
1493    

Legend:
Removed from v.1122  
changed lines
  Added in v.1514

  ViewVC Help
Powered by ViewVC 1.1.26