/[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 1481 by artak, Wed Apr 9 00:45:47 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=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*b
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       v.append(r/rho)
614       g[0]=rho
615    
616       while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
617    
618        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
619    
620        
621        p=Msolve(Aprod(v[iter]))
622    
623        v.append(p)
624    
625        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
626    
627    # Modified Gram-Schmidt
628        for j in range(iter+1):
629          h[j][iter]=bilinearform(v[j],v[iter+1])  
630          v[iter+1]+=(-1.)*h[j][iter]*v[j]
631          
632        h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
633        v_norm2=h[iter+1][iter]
634    
635    
636    # Reorthogonalize if needed
637        if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
638         for j in range(iter+1):
639            hr=bilinearform(v[j],v[iter+1])
640                h[j][iter]=h[j][iter]+hr #vhat
641                v[iter+1] +=(-1.)*hr*v[j]
642    
643         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
644         h[iter+1][iter]=v_norm2
645    
646    #   watch out for happy breakdown
647            if v_norm2 != 0:
648             v[iter+1]=v[iter+1]/h[iter+1][iter]
649    
650    #   Form and store the information for the new Givens rotation
651        if iter > 0 :
652            hhat=[]
653            for i in range(iter+1) : hhat.append(h[i][iter])
654            hhat=givapp(c[0:iter],s[0:iter],hhat);
655                for i in range(iter+1) : h[i][iter]=hhat[i]
656    
657        mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])
658        if mu!=0 :
659            c[iter]=h[iter][iter]/mu
660            s[iter]=-h[iter+1][iter]/mu
661            h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]
662            h[iter+1][iter]=0.0
663            g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2])
664    
665    # Update the residual norm
666            rho=abs(g[iter+1])
667        iter+=1
668    
669    # At this point either iter > iter_max or rho < tol.
670    # It's time to compute x and leave.        
671    
672       if iter > 0 :
673         y=numarray.zeros(iter,numarray.Float64)    
674         y[iter-1] = g[iter-1] / h[iter-1][iter-1]
675         if iter > 1 :  
676            i=iter-2  
677            while i>=0 :
678              y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]
679              i=i-1
680         xhat=v[iter-1]*y[iter-1]
681         for i in range(iter-1):
682        xhat += v[i]*y[i]
683       else : xhat=v[0]
684        
685       x += xhat
686       if iter!=iter_restart-1:
687          stopped=True
688       else:
689          stopped=False
690    
691       return x,stopped
692    
693    def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
694    
695        #
696        #  minres solves the system of linear equations Ax = b
697        #  where A is a symmetric matrix (possibly indefinite or singular)
698        #  and b is a given vector.
699        #  
700        #  "A" may be a dense or sparse matrix (preferably sparse!)
701        #  or the name of a function such that
702        #               y = A(x)
703        #  returns the product y = Ax for any given vector x.
704        #
705        #  "M" defines a positive-definite preconditioner M = C C'.
706        #  "M" may be a dense or sparse matrix (preferably sparse!)
707        #  or the name of a function such that
708        #  solves the system My = x for any given vector x.
709        #
710        #
711    
712        #  Initialize                              
713    
714        iter   = 0
715        Anorm = 0
716        ynorm = 0
717        x=x*0
718        #------------------------------------------------------------------
719        # Set up y and v for the first Lanczos vector v1.
720        # y  =  beta1 P' v1,  where  P = C**(-1).
721        # v is really P' v1.
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)       # Normalize y to get v1 later.
734    
735        #------------------------------------------------------------------
736        # Initialize other quantities.
737        # ------------------------------------------------------------------
738        oldb   = 0
739        beta   = beta1
740        dbar   = 0
741        epsln  = 0
742        phibar = beta1
743        rhs1   = beta1
744        rhs2   = 0
745        rnorm  = phibar
746        tnorm2 = 0
747        ynorm2 = 0
748        cs     = -1
749        sn     = 0
750        w      = b*0.
751        w2     = b*0.
752        r2     = r1
753        eps    = 0.0001
754    
755        #---------------------------------------------------------------------
756        # Main iteration loop.
757        # --------------------------------------------------------------------
758        while not stoppingcriterium(rnorm,Anorm*ynorm):    #  ||r|| / (||A|| ||x||)
759    
760        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
761            iter    = iter  +  1
762    
763            #-----------------------------------------------------------------
764            # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
765            # The general iteration is similar to the case k = 1 with v0 = 0:
766            #
767            #   p1      = Operator * v1  -  beta1 * v0,
768            #   alpha1  = v1'p1,
769            #   q2      = p2  -  alpha1 * v1,
770            #   beta2^2 = q2'q2,
771            #   v2      = (1/beta2) q2.
772            #
773            # Again, y = betak P vk,  where  P = C**(-1).
774            #-----------------------------------------------------------------
775            s = 1/beta                 # Normalize previous vector (in y).
776            v = s*y                    # v = vk if P = I
777        
778            y      = Aprod(v)
779        
780            if iter >= 2:
781              y = y - (beta/oldb)*r1
782    
783            alfa   = bilinearform(v,y)              # alphak
784            y      = (- alfa/beta)*r2 + y
785            r1     = r2
786            r2     = y
787            y = Msolve(r2)
788            oldb   = beta                           # oldb = betak
789            beta   = bilinearform(r2,y)             # beta = betak+1^2
790            if beta < 0: raise NegativeNorm,"negative norm."
791    
792            beta   = math.sqrt( beta )
793            tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
794            
795            if iter==1:                 # Initialize a few things.
796              gmax   = abs( alfa )      # alpha1
797              gmin   = gmax             # alpha1
798    
799            # Apply previous rotation Qk-1 to get
800            #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
801            #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
802        
803            oldeps = epsln
804            delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
805            gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
806            epsln  =               sn * beta  # epsln2 = 0         epslnk+1
807            dbar   =            -  cs * beta  # dbar 2 = beta2     dbar k+1
808    
809            # Compute the next plane rotation Qk
810    
811            gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
812            gamma  = max(gamma,eps)
813            cs     = gbar / gamma             # ck
814            sn     = beta / gamma             # sk
815            phi    = cs * phibar              # phik
816            phibar = sn * phibar              # phibark+1
817    
818            # Update  x.
819    
820            denom = 1/gamma
821            w1    = w2
822            w2    = w
823            w     = (v - oldeps*w1 - delta*w2) * denom
824            x     = x  +  phi*w
825    
826            # Go round again.
827    
828            gmax   = max(gmax,gamma)
829            gmin   = min(gmin,gamma)
830            z      = rhs1 / gamma
831            ynorm2 = z*z  +  ynorm2
832            rhs1   = rhs2 -  delta*z
833            rhs2   =      -  epsln*z
834    
835            # Estimate various norms and test for convergence.
836    
837            Anorm  = math.sqrt( tnorm2 )
838            ynorm  = math.sqrt( ynorm2 )
839    
840            rnorm  = phibar
841    
842        # Return final answer.
843        return x
844        
845    #############################################
846    
847    class ArithmeticTuple(object):
848       """
849       tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float.
850    
851       example of usage:
852    
853       from esys.escript import Data
854       from numarray import array
855       a=Data(...)
856       b=array([1.,4.])
857       x=ArithmeticTuple(a,b)
858       y=5.*x
859    
860       """
861       def __init__(self,*args):
862           """
863           initialize object with elements args.
864    
865           @param args: tuple of object that support implace add (x+=y) and scaling (x=a*y)
866           """
867           self.__items=list(args)
868    
869       def __len__(self):
870           """
871           number of items
872    
873           @return: number of items
874           @rtype: C{int}
875           """
876           return len(self.__items)
877    
878       def __getitem__(self,index):
879           """
880           get an item
881    
882           @param index: item to be returned
883           @type index: C{int}
884           @return: item with index C{index}
885           """
886           return self.__items.__getitem__(index)
887    
888       def __mul__(self,other):
889           """
890           scaling from the right
891    
892           @param other: scaling factor
893           @type other: C{float}
894           @return: itemwise self*other
895           @rtype: L{ArithmeticTuple}
896           """
897           out=[]
898           for i in range(len(self)):
899               out.append(self[i]*other)
900           return ArithmeticTuple(*tuple(out))
901    
902       def __rmul__(self,other):
903           """
904           scaling from the left
905    
906           @param other: scaling factor
907           @type other: C{float}
908           @return: itemwise other*self
909           @rtype: L{ArithmeticTuple}
910           """
911           out=[]
912           for i in range(len(self)):
913               out.append(other*self[i])
914           return ArithmeticTuple(*tuple(out))
915    
916    #########################
917    # Added by Artak
918    #########################
919       def __div__(self,other):
920           """
921           dividing from the right
922    
923           @param other: scaling factor
924           @type other: C{float}
925           @return: itemwise self/other
926           @rtype: L{ArithmeticTuple}
927           """
928           out=[]
929           for i in range(len(self)):
930               out.append(self[i]/other)
931           return ArithmeticTuple(*tuple(out))
932    
933       def __rdiv__(self,other):
934           """
935           dividing from the left
936    
937           @param other: scaling factor
938           @type other: C{float}
939           @return: itemwise other/self
940           @rtype: L{ArithmeticTuple}
941           """
942           out=[]
943           for i in range(len(self)):
944               out.append(other/self[i])
945           return ArithmeticTuple(*tuple(out))
946      
947    ##########################################33
948    
949       def __iadd__(self,other):
950           """
951           in-place add of other to self
952    
953           @param other: increment
954           @type other: C{ArithmeticTuple}
955           """
956           if len(self) != len(other):
957               raise ValueError,"tuple length must match."
958           for i in range(len(self)):
959               self.__items[i]+=other[i]
960           return self
961    
962    class HomogeneousSaddlePointProblem(object):
963          """
964          This provides a framwork for solving homogeneous saddle point problem of the form
965    
966                 Av+B^*p=f
967                 Bv    =0
968    
969          for the unknowns v and p and given operators A and B and given right hand side f.
970          B^* is the adjoint operator of B is the given inner product.
971    
972          """
973          def __init__(self,**kwargs):
974            self.setTolerance()
975            self.setToleranceReductionFactor()
976    
977          def initialize(self):
978            """
979            initialize the problem (overwrite)
980            """
981            pass
982          def B(self,v):
983             """
984             returns Bv (overwrite)
985             @rtype: equal to the type of p
986    
987             @note: boundary conditions on p should be zero!
988             """
989             pass
990    
991          def inner(self,p0,p1):
992             """
993             returns inner product of two element p0 and p1  (overwrite)
994            
995             @type p0: equal to the type of p
996             @type p1: equal to the type of p
997             @rtype: C{float}
998    
999             @rtype: equal to the type of p
1000             """
1001             pass
1002    
1003          def solve_A(self,u,p):
1004             """
1005             solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)
1006    
1007             @rtype: equal to the type of v
1008             @note: boundary conditions on v should be zero!
1009             """
1010             pass
1011    
1012          def solve_prec(self,p):
1013             """
1014             provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)
1015    
1016             @rtype: equal to the type of p
1017             """
1018             pass
1019    
1020          def stoppingcriterium(self,Bv,v,p):
1021             """
1022             returns a True if iteration is terminated. (overwrite)
1023    
1024             @rtype: C{bool}
1025             """
1026             pass
1027                
1028          def __inner(self,p,r):
1029             return self.inner(p,r[1])
1030    
1031          def __inner_p(self,p1,p2):
1032             return self.inner(p1,p2)
1033    
1034          def __stoppingcriterium(self,norm_r,r,p):
1035              return self.stoppingcriterium(r[1],r[0],p)
1036    
1037          def __stoppingcriterium_GMRES(self,norm_r,norm_b):
1038              return self.stoppingcriterium_GMRES(norm_r,norm_b)
1039    
1040          def __stoppingcriterium_MINRES(self,norm_r,norm_Ax):
1041              return self.stoppingcriterium_MINRES(norm_r,norm_Ax)
1042    
1043    
1044          def setTolerance(self,tolerance=1.e-8):
1045                  self.__tol=tolerance
1046          def getTolerance(self):
1047                  return self.__tol
1048          def setToleranceReductionFactor(self,reduction=0.01):
1049                  self.__reduction=reduction
1050          def getSubProblemTolerance(self):
1051                  return self.__reduction*self.getTolerance()
1052    
1053          def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG'):
1054                  """
1055                  solves the saddle point problem using initial guesses v and p.
1056    
1057                  @param max_iter: maximum number of iteration steps.
1058                  """
1059                  self.verbose=verbose
1060                  self.show_details=show_details and self.verbose
1061    
1062                  # assume p is known: then v=A^-1(f-B^*p)
1063                  # which leads to BA^-1B^*p = BA^-1f  
1064    
1065              # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)      
1066    
1067              
1068              self.__z=v+self.solve_A(v,p*0)
1069    
1070                  Bz=self.B(self.__z)
1071                  #
1072              #   solve BA^-1B^*p = Bz
1073                  #
1074                  #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1075                  #
1076                  #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1077                  #                           A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)
1078                  #
1079                  self.iter=0
1080              if solver=='GMRES':      
1081                    if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter
1082                    p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1.)
1083                    # solve Au=f-B^*p
1084                    #       A(u-v)=f-B^*p-Av
1085                    #       u=v+(u-v)
1086            u=v+self.solve_A(v,p)
1087    
1088              if solver=='MINRES':      
1089                    if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter
1090                    p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_MINRES,iter_max=max_iter, x=p*1.)
1091                    # solve Au=f-B^*p
1092                    #       A(u-v)=f-B^*p-Av
1093                    #       u=v+(u-v)
1094            u=v+self.solve_A(v,p)
1095        
1096                  if solver=='PCG':
1097                    if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1098                    p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1099                u=r[0]  
1100    
1101                  print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1102    
1103              return u,p
1104    
1105          def __Msolve(self,r):
1106              return self.solve_prec(r[1])
1107    
1108          def __Msolve_GMRES(self,r):
1109              return self.solve_prec(r)
1110    
1111    
1112          def __Aprod(self,p):
1113              # return BA^-1B*p
1114              #solve Av =-B^*p as Av =f-Az-B^*p
1115              v=self.solve_A(self.__z,-p)
1116              return ArithmeticTuple(v, self.B(v))
1117    
1118          def __Aprod_GMRES(self,p):
1119              # return BA^-1B*p
1120              #solve Av =-B^*p as Av =f-Az-B^*p
1121          v=self.solve_A(self.__z,-p)
1122              return self.B(v)
1123    
1124  class SaddlePointProblem(object):  class SaddlePointProblem(object):
1125     """     """
1126     This implements a solver for a saddlepoint problem     This implements a solver for a saddlepoint problem
# Line 580  class SaddlePointProblem(object): Line 1308  class SaddlePointProblem(object):
1308              norm_u_new = util.Lsup(u_new)              norm_u_new = util.Lsup(u_new)
1309              p_new=p+self.relaxation*g_new              p_new=p+self.relaxation*g_new
1310              norm_p_new = util.sqrt(self.inner(p_new,p_new))              norm_p_new = util.sqrt(self.inner(p_new,p_new))
1311              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))
1312    
1313              if self.iter>1:              if self.iter>1:
1314                 dg2=g_new-g                 dg2=g_new-g
# Line 622  class SaddlePointProblem(object): Line 1350  class SaddlePointProblem(object):
1350  #  #
1351  #      return u,p  #      return u,p
1352                        
1353  # vim: expandtab shiftwidth=4:  def MaskFromBoundaryTag(function_space,*tags):
1354       """
1355       create a mask on the given function space which one for samples
1356       that touch the boundary tagged by tags.
1357    
1358       usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")
1359    
1360       @param function_space: a given function space
1361       @type function_space: L{escript.FunctionSpace}
1362       @param tags: boundray tags
1363       @type tags: C{str}
1364       @return: a mask which marks samples used by C{function_space} that are touching the
1365                boundary tagged by any of the given tags.
1366       @rtype: L{escript.Data} of rank 0
1367       """
1368       pde=linearPDEs.LinearPDE(function_space.getDomain(),numEquations=1, numSolutions=1)
1369       d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
1370       for t in tags: d.setTaggedValue(t,1.)
1371       pde.setValue(y=d)
1372       out=util.whereNonZero(pde.getRightHandSide())
1373       if out.getFunctionSpace() == function_space:
1374          return out
1375       else:
1376          return util.whereNonZero(util.interpolate(out,function_space))
1377    
1378    
1379    

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

  ViewVC Help
Powered by ViewVC 1.1.26