/[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 2061 by jfenwick, Wed Nov 19 03:40:21 2008 UTC revision 2100 by gross, Wed Nov 26 08:13:00 2008 UTC
# Line 426  class MaxIterReached(SolverSchemeExcepti Line 426  class MaxIterReached(SolverSchemeExcepti
426     maxium number of iteration steps is reached.     maxium number of iteration steps is reached.
427     """     """
428     pass     pass
429    class CorrectionFailed(SolverSchemeException):
430       """
431       no convergence has been achieved in the solution correction scheme.
432       """
433       pass
434  class IterationBreakDown(SolverSchemeException):  class IterationBreakDown(SolverSchemeException):
435     """     """
436     iteration scheme econouters an incurable breakdown.     iteration scheme econouters an incurable breakdown.
# Line 437  class NegativeNorm(SolverSchemeException Line 442  class NegativeNorm(SolverSchemeException
442     """     """
443     pass     pass
444    
445  class IterationHistory(object):  def PCG(b, Aprod, Msolve, bilinearform, atol=0, rtol=1.e-8, x=None, iter_max=100, initial_guess=True, verbose=False):
    """  
    The IterationHistory class is used to define a stopping criterium. It keeps track of the  
    residual norms. The stoppingcriterium indicates termination if the residual norm has been reduced by  
    a given tolerance.  
    """  
    def __init__(self,tolerance=math.sqrt(util.EPSILON),verbose=False):  
       """  
       Initialization  
   
       @param tolerance: tolerance  
       @type tolerance: positive C{float}  
       @param verbose: switches on the printing out some information  
       @type verbose: C{bool}  
       """  
       if not tolerance>0.:  
           raise ValueError,"tolerance needs to be positive."  
       self.tolerance=tolerance  
       self.verbose=verbose  
       self.history=[]  
    def stoppingcriterium(self,norm_r,r,x):  
        """  
        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.  
   
         
        @param norm_r: current residual norm  
        @type norm_r: non-negative C{float}  
        @param r: current residual (not used)  
        @param x: current solution approximation (not used)  
        @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.  
        @rtype: C{bool}  
   
        """  
        self.history.append(norm_r)  
        if self.verbose: print "iter: %s:  inner(rhat,r) = %e"#(len(self.history)-1, self.history[-1])  
        return self.history[-1]<=self.tolerance * self.history[0]  
   
    def stoppingcriterium2(self,norm_r,norm_b,solver="GMRES",TOL=None):  
        """  
        returns True if the C{norm_r} is C{tolerance}*C{norm_b}  
   
         
        @param norm_r: current residual norm  
        @type norm_r: non-negative C{float}  
        @param norm_b: norm of right hand side  
        @type norm_b: non-negative C{float}  
        @return: C{True} is the stopping criterium is fullfilled. Otherwise C{False} is returned.  
        @rtype: C{bool}  
   
        """  
        if TOL==None:  
           TOL=self.tolerance  
        self.history.append(norm_r)  
        if self.verbose: print "iter: %s:  norm(r) = %e"#(len(self.history)-1, self.history[-1])  
        return self.history[-1]<=TOL * norm_b  
   
 def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  
446     """     """
447     Solver for     Solver for
448    
# Line 537  def PCG(b, Aprod, Msolve, bilinearform, Line 486  def PCG(b, Aprod, Msolve, bilinearform,
486     iter=0     iter=0
487     if x==None:     if x==None:
488        x=0*b        x=0*b
489          r=b*1.
490          initial_guess=False
491       if initial_guess:
492          r=b+(-1)*Aprod(x)
493          x=x*1.
494     else:     else:
495        b += (-1)*Aprod(x)        x=x*0
496     r=b        r=b*1
497     rhat=Msolve(r)     rhat=Msolve(r)
498     d = rhat     d = rhat
499     rhat_dot_r = bilinearform(rhat, r)     rhat_dot_r = bilinearform(rhat, r)
500     if rhat_dot_r<0: raise NegativeNorm,"negative norm."     if rhat_dot_r<0: raise NegativeNorm,"negative norm."
501       norm_b=math.sqrt(rhat_dot_r)
502       if verbose: print "PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_b, atol+rtol*norm_b)
503    
504     while not stoppingcriterium(math.sqrt(rhat_dot_r),r,x):     while not math.sqrt(rhat_dot_r) <= atol+rtol*norm_b:
505         iter+=1         iter+=1
506         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
507    
# Line 562  def PCG(b, Aprod, Msolve, bilinearform, Line 518  def PCG(b, Aprod, Msolve, bilinearform,
518    
519         rhat_dot_r = rhat_dot_r_new         rhat_dot_r = rhat_dot_r_new
520         if rhat_dot_r<0: raise NegativeNorm,"negative norm."         if rhat_dot_r<0: raise NegativeNorm,"negative norm."
521           if verbose: print "PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))
522       if verbose: print "PCG: tolerance reached after %s steps."%iter
523     return x,r     return x,r
524    
525  class Defect(object):  class Defect(object):
# Line 776  def __FDGMRES(F0, defect, x0, atol, iter Line 733  def __FDGMRES(F0, defect, x0, atol, iter
733    
734          # Modified Gram-Schmidt          # Modified Gram-Schmidt
735      for j in range(iter+1):      for j in range(iter+1):
736           h[j][iter]=defect.bilinearform(v[j],v[iter+1])             h[j,iter]=defect.bilinearform(v[j],v[iter+1])  
737           v[iter+1]-=h[j][iter]*v[j]           v[iter+1]-=h[j,iter]*v[j]
738                
739      h[iter+1][iter]=defect.norm(v[iter+1])      h[iter+1,iter]=defect.norm(v[iter+1])
740      v_norm2=h[iter+1][iter]      v_norm2=h[iter+1,iter]
741    
742          # Reorthogonalize if needed          # Reorthogonalize if needed
743      if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)      if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
744          for j in range(iter+1):            for j in range(iter+1):  
745             hr=defect.bilinearform(v[j],v[iter+1])             hr=defect.bilinearform(v[j],v[iter+1])
746                 h[j][iter]=h[j][iter]+hr                 h[j,iter]=h[j,iter]+hr
747                 v[iter+1] -= hr*v[j]                 v[iter+1] -= hr*v[j]
748    
749          v_norm2=defect.norm(v[iter+1])          v_norm2=defect.norm(v[iter+1])
750          h[iter+1][iter]=v_norm2          h[iter+1,iter]=v_norm2
751          #   watch out for happy breakdown          #   watch out for happy breakdown
752          if not v_norm2 == 0:          if not v_norm2 == 0:
753                  v[iter+1]=v[iter+1]/h[iter+1][iter]                  v[iter+1]=v[iter+1]/h[iter+1,iter]
754    
755          #   Form and store the information for the new Givens rotation          #   Form and store the information for the new Givens rotation
756      if iter > 0 :      if iter > 0 :
757          hhat=numarray.zeros(iter+1,numarray.Float64)          hhat=numarray.zeros(iter+1,numarray.Float64)
758          for i in range(iter+1) : hhat[i]=h[i][iter]          for i in range(iter+1) : hhat[i]=h[i,iter]
759          hhat=__givapp(c[0:iter],s[0:iter],hhat);          hhat=__givapp(c[0:iter],s[0:iter],hhat);
760              for i in range(iter+1) : h[i][iter]=hhat[i]              for i in range(iter+1) : h[i,iter]=hhat[i]
761    
762      mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])      mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
763    
764      if mu!=0 :      if mu!=0 :
765          c[iter]=h[iter][iter]/mu          c[iter]=h[iter,iter]/mu
766          s[iter]=-h[iter+1][iter]/mu          s[iter]=-h[iter+1,iter]/mu
767          h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]          h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
768          h[iter+1][iter]=0.0          h[iter+1,iter]=0.0
769          g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])          g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])
770    
771          # Update the residual norm          # Update the residual norm
# Line 819  def __FDGMRES(F0, defect, x0, atol, iter Line 776  def __FDGMRES(F0, defect, x0, atol, iter
776     # It's time to compute x and leave.             # It's time to compute x and leave.        
777     if iter > 0 :     if iter > 0 :
778       y=numarray.zeros(iter,numarray.Float64)           y=numarray.zeros(iter,numarray.Float64)    
779       y[iter-1] = g[iter-1] / h[iter-1][iter-1]       y[iter-1] = g[iter-1] / h[iter-1,iter-1]
780       if iter > 1 :         if iter > 1 :  
781          i=iter-2            i=iter-2  
782          while i>=0 :          while i>=0 :
783            y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]            y[i] = ( g[i] - numarray.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
784            i=i-1            i=i-1
785       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
786       for i in range(iter-1):       for i in range(iter-1):
# Line 839  def __FDGMRES(F0, defect, x0, atol, iter Line 796  def __FDGMRES(F0, defect, x0, atol, iter
796     return xhat,stopped     return xhat,stopped
797    
798  ##############################################  ##############################################
799  def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):  def GMRES(b, Aprod, bilinearform, atol=0, rtol=1.e-8, x=None, iter_max=100, iter_restart=20, verbose=False):
800  ################################################  ################################################
801     m=iter_restart     m=iter_restart
802     iter=0     iter=0
803     xc=x     if x == None:
804          xc=b*0
805       else:
806          xc=x*1
807       if rtol>0:
808          r_dot_r = bilinearform(b, b)
809          if r_dot_r<0: raise NegativeNorm,"negative norm."
810          atol2=atol+rtol*math.sqrt(r_dot_r)
811          if verbose: print "GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)
812       else:
813          atol2=atol
814          if verbose: print "GMRES: absolute tolerance = %e"%atol2
815      
816     while True:     while True:
817        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
818        xc,stopped=__GMRESm(b*1, Aprod, Msolve, bilinearform, stoppingcriterium, x=xc*1, iter_max=iter_max-iter, iter_restart=m)        xc,stopped=__GMRESm(b, Aprod, bilinearform, atol2, x=xc, iter_max=iter_max-iter, iter_restart=m, verbose=verbose)
       if stopped: break  
819        iter+=iter_restart            iter+=iter_restart    
820          if stopped: break
821          if verbose: print "GMRES: restart."
822       if verbose: print "GMRES: tolerance has reached."
823     return xc     return xc
824    
825  def __GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):  def __GMRESm(b, Aprod, bilinearform, atol, x, iter_max=100, iter_restart=20, verbose=False):
826     iter=0     iter=0
    r=Msolve(b)  
    r_dot_r = bilinearform(r, r)  
    if r_dot_r<0: raise NegativeNorm,"negative norm."  
    norm_b=math.sqrt(r_dot_r)  
   
    if x==None:  
       x=0*b  
    else:  
       r=Msolve(b-Aprod(x))  
       r_dot_r = bilinearform(r, r)  
       if r_dot_r<0: raise NegativeNorm,"negative norm."  
827        
828     h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)     h=numarray.zeros((iter_restart+1,iter_restart),numarray.Float64)
829     c=numarray.zeros(iter_restart,numarray.Float64)     c=numarray.zeros(iter_restart,numarray.Float64)
830     s=numarray.zeros(iter_restart,numarray.Float64)     s=numarray.zeros(iter_restart,numarray.Float64)
831     g=numarray.zeros(iter_restart,numarray.Float64)     g=numarray.zeros(iter_restart+1,numarray.Float64)
832     v=[]     v=[]
833    
834       r=b-Aprod(x)
835       r_dot_r = bilinearform(r, r)
836       if r_dot_r<0: raise NegativeNorm,"negative norm."
837     rho=math.sqrt(r_dot_r)     rho=math.sqrt(r_dot_r)
838        
839     v.append(r/rho)     v.append(r/rho)
840     g[0]=rho     g[0]=rho
841    
842     while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):     if verbose: print "GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)
843       while not (rho<=atol or iter==iter_restart):
844    
845      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
846    
847      p=Msolve(Aprod(v[iter]))      p=Aprod(v[iter])
   
848      v.append(p)      v.append(p)
849    
850      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
851    
852  # Modified Gram-Schmidt  # Modified Gram-Schmidt
853      for j in range(iter+1):      for j in range(iter+1):
854        h[j][iter]=bilinearform(v[j],v[iter+1])          h[j,iter]=bilinearform(v[j],v[iter+1])  
855        v[iter+1]-=h[j][iter]*v[j]        v[iter+1]-=h[j,iter]*v[j]
856                
857      h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))      h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
858      v_norm2=h[iter+1][iter]      v_norm2=h[iter+1,iter]
859    
860  # Reorthogonalize if needed  # Reorthogonalize if needed
861      if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)      if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
862       for j in range(iter+1):         for j in range(iter+1):  
863          hr=bilinearform(v[j],v[iter+1])          hr=bilinearform(v[j],v[iter+1])
864              h[j][iter]=h[j][iter]+hr              h[j,iter]=h[j,iter]+hr
865              v[iter+1] -= hr*v[j]              v[iter+1] -= hr*v[j]
866    
867       v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
868       h[iter+1][iter]=v_norm2       h[iter+1,iter]=v_norm2
869    
870  #   watch out for happy breakdown  #   watch out for happy breakdown
871          if not v_norm2 == 0:          if not v_norm2 == 0:
872           v[iter+1]=v[iter+1]/h[iter+1][iter]           v[iter+1]=v[iter+1]/h[iter+1,iter]
873    
874  #   Form and store the information for the new Givens rotation  #   Form and store the information for the new Givens rotation
875      if iter > 0 :      if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])
876          hhat=numarray.zeros(iter+1,numarray.Float64)      mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
         for i in range(iter+1) : hhat[i]=h[i][iter]  
         hhat=__givapp(c[0:iter],s[0:iter],hhat);  
             for i in range(iter+1) : h[i][iter]=hhat[i]  
   
     mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter])  
877    
878      if mu!=0 :      if mu!=0 :
879          c[iter]=h[iter][iter]/mu          c[iter]=h[iter,iter]/mu
880          s[iter]=-h[iter+1][iter]/mu          s[iter]=-h[iter+1,iter]/mu
881          h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter]          h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
882          h[iter+1][iter]=0.0          h[iter+1,iter]=0.0
883          g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])          g[iter:iter+2]=__givapp(c[iter],s[iter],g[iter:iter+2])
   
884  # Update the residual norm  # Update the residual norm
885                                
886          rho=abs(g[iter+1])          rho=abs(g[iter+1])
887            if verbose: print "GMRES: iteration step %s: residual %e"%(iter,rho)
888      iter+=1      iter+=1
889    
890  # At this point either iter > iter_max or rho < tol.  # At this point either iter > iter_max or rho < tol.
891  # It's time to compute x and leave.          # It's time to compute x and leave.        
892    
893       if verbose: print "GMRES: iteration stopped after %s step."%iter
894     if iter > 0 :     if iter > 0 :
895       y=numarray.zeros(iter,numarray.Float64)           y=numarray.zeros(iter,numarray.Float64)    
896       y[iter-1] = g[iter-1] / h[iter-1][iter-1]       y[iter-1] = g[iter-1] / h[iter-1,iter-1]
897       if iter > 1 :         if iter > 1 :  
898          i=iter-2            i=iter-2  
899          while i>=0 :          while i>=0 :
900            y[i] = ( g[i] - numarray.dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i]            y[i] = ( g[i] - numarray.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
901            i=i-1            i=i-1
902       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
903       for i in range(iter-1):       for i in range(iter-1):
904      xhat += v[i]*y[i]      xhat += v[i]*y[i]
905     else : xhat=v[0]     else:
906         xhat=v[0] * 0
907    
908     x += xhat     x += xhat
909     if iter<iter_restart-1:     if iter<iter_restart-1:
# Line 954  def __GMRESm(b, Aprod, Msolve, bilinearf Line 914  def __GMRESm(b, Aprod, Msolve, bilinearf
914     return x,stopped     return x,stopped
915    
916  #################################################  #################################################
917  def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  def MINRES(b, Aprod, Msolve, bilinearform, atol=0, rtol=1.e-8, x=None, iter_max=100):
918  #################################################  #################################################
919      #      #
920      #  minres solves the system of linear equations Ax = b      #  minres solves the system of linear equations Ax = b
# Line 1021  def MINRES(b, Aprod, Msolve, bilinearfor Line 981  def MINRES(b, Aprod, Msolve, bilinearfor
981      #---------------------------------------------------------------------      #---------------------------------------------------------------------
982      # Main iteration loop.      # Main iteration loop.
983      # --------------------------------------------------------------------      # --------------------------------------------------------------------
984      while not stoppingcriterium(rnorm,Anorm*ynorm,'MINRES'):    #  checks ||r|| < (||A|| ||x||) * TOL      while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL
985    
986      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
987          iter    = iter  +  1          iter    = iter  +  1
# Line 1107  def MINRES(b, Aprod, Msolve, bilinearfor Line 1067  def MINRES(b, Aprod, Msolve, bilinearfor
1067    
1068      return x      return x
1069    
1070  def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  def TFQMR(b, Aprod, Msolve, bilinearform, atol=0, rtol=1.e-8, x=None, iter_max=100):
1071    
1072  # TFQMR solver for linear systems  # TFQMR solver for linear systems
1073  #  #
# Line 1144  def TFQMR(b, Aprod, Msolve, bilinearform Line 1104  def TFQMR(b, Aprod, Msolve, bilinearform
1104    tau = math.sqrt(bilinearform(r,r))    tau = math.sqrt(bilinearform(r,r))
1105    error = [ error, tau ]    error = [ error, tau ]
1106    rho = tau * tau    rho = tau * tau
   m=1  
1107  #  #
1108  #  TFQMR iteration  #  TFQMR iteration
1109  #  #
1110  #  while ( iter < kmax-1 ):  #  while ( iter < kmax-1 ):
1111        
1112    while not stoppingcriterium(tau*math.sqrt ( m + 1 ),norm_b,'TFQMR'):    while not tau<=atol+rtol*norm_b:
1113      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
1114    
1115      sigma = bilinearform(r,v)      sigma = bilinearform(r,v)
# Line 1399  class HomogeneousSaddlePointProblem(obje Line 1358  class HomogeneousSaddlePointProblem(obje
1358               Bv    =0               Bv    =0
1359    
1360        for the unknowns v and p and given operators A and B and given right hand side f.        for the unknowns v and p and given operators A and B and given right hand side f.
1361        B^* is the adjoint operator of B is the given inner product.        B^* is the adjoint operator of B.
   
1362        """        """
1363        def __init__(self,**kwargs):        def __init__(self,**kwargs):
1364          self.setTolerance()          self.setTolerance()
1365          self.setToleranceReductionFactor()          self.setAbsoluteTolerance()
1366            self.setSubToleranceReductionFactor()
1367    
1368          #=============================================================
1369        def initialize(self):        def initialize(self):
1370          """          """
1371          initialize the problem (overwrite)          initialize the problem (overwrite)
1372          """          """
1373          pass          pass
1374    
1375        def B(self,v):        def B(self,v):
1376            """
1377            returns Bv (overwrite)
1378            @rtype: equal to the type of p
1379    
1380            @note: boundary conditions on p should be zero!
1381            """
1382            raise NotImplementedError,"no operator B implemented."
1383    
1384          def inner_pBv(self,p,Bv):
1385           """           """
1386           returns Bv (overwrite)           returns inner product of element p and Bv  (overwrite)
1387           @rtype: equal to the type of p          
1388             @type p: equal to the type of p
1389             @type Bv: equal to the type of result of operator B
1390             @rtype: C{float}
1391    
1392           @note: boundary conditions on p should be zero!           @rtype: equal to the type of p
1393           """           """
1394           pass           raise NotImplementedError,"no inner product for p implemented."
1395    
1396        def inner(self,p0,p1):        def inner_p(self,p0,p1):
1397           """           """
1398           returns inner product of two element p0 and p1  (overwrite)           returns inner product of element p0 and p1  (overwrite)
1399                    
1400           @type p0: equal to the type of p           @type p0: equal to the type of p
1401           @type p1: equal to the type of p           @type p1: equal to the type of p
1402           @rtype: equal to the type of p           @rtype: equal to the type of p
1403           """           """
1404           pass           raise NotImplementedError,"no inner product for p implemented."
1405    
1406        def solve_A(self,u,p):        def inner_v(self,v0,v1):
1407           """           """
1408           solves Av=f-Au-B^*p with accuracy self.getReducedTolerance() (overwrite)           returns inner product of two element v0 and v1  (overwrite)
1409            
1410             @type v0: equal to the type of v
1411             @type v1: equal to the type of v
1412             @rtype: C{float}
1413    
1414           @rtype: equal to the type of v           @rtype: equal to the type of v
          @note: boundary conditions on v should be zero!  
1415           """           """
1416             raise NotImplementedError,"no inner product for v implemented."
1417           pass           pass
1418    
1419        def solve_prec(self,p):        def solve_A(self,u,p):
1420           """           """
1421           provides a preconditioner for BA^{-1}B^* with accuracy self.getReducedTolerance() (overwrite)           solves Av=f-Au-B^*p with accuracy self.getSubProblemTolerance() (overwrite)
1422    
1423           @rtype: equal to the type of p           @rtype: equal to the type of v
1424             @note: boundary conditions on v should be zero!
1425           """           """
1426           pass           raise NotImplementedError,"no operator A implemented."
1427    
1428        def stoppingcriterium(self,Bv,v,p):        def solve_prec(self,p):
1429           """           """
1430           returns a True if iteration is terminated. (overwrite)           provides a preconditioner for BA^{-1}B^* with accuracy self.getSubProblemTolerance() (overwrite)
1431    
1432           @rtype: C{bool}           @rtype: equal to the type of p
1433             @note: boundary conditions on p should be zero!
1434           """           """
1435           pass           raise NotImplementedError,"no preconditioner for Schur complement implemented."
1436                      #=============================================================
1437        def __inner(self,p,r):        def __Aprod_PCG(self,p):
1438           return self.inner(p,r[1])            # return (v,Bv) with v=A^-1B*p
1439              #solve Av =B^*p as Av =f-Az-B^*(-p)
1440        def __inner_p(self,p1,p2):            v=self.solve_A(self.__z,-p)
1441           return self.inner(p1,p2)            return ArithmeticTuple(v, self.B(v))
         
       def __inner_a(self,a1,a2):  
          return self.inner_a(a1,a2)  
   
       def __inner_a1(self,a1,a2):  
          return self.inner(a1[1],a2[1])  
   
       def __stoppingcriterium(self,norm_r,r,p):  
           return self.stoppingcriterium(r[1],r[0],p)  
   
       def __stoppingcriterium2(self,norm_r,norm_b,solver='GMRES',TOL=None):  
           return self.stoppingcriterium2(norm_r,norm_b,solver,TOL)  
   
       def setTolerance(self,tolerance=1.e-8):  
               self.__tol=tolerance  
       def getTolerance(self):  
               return self.__tol  
       def setToleranceReductionFactor(self,reduction=0.01):  
               self.__reduction=reduction  
       def getSubProblemTolerance(self):  
               return self.__reduction*self.getTolerance()  
   
       def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='PCG',iter_restart=20):  
               """  
               solves the saddle point problem using initial guesses v and p.  
   
               @param max_iter: maximum number of iteration steps.  
               """  
               self.verbose=verbose  
               self.show_details=show_details and self.verbose  
   
               # assume p is known: then v=A^-1(f-B^*p)  
               # which leads to BA^-1B^*p = BA^-1f    
   
           # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)        
           self.__z=v+self.solve_A(v,p*0)  
               Bz=self.B(self.__z)  
               #  
           #   solve BA^-1B^*p = Bz  
               #  
               #  
               #  
               self.iter=0  
           if solver=='GMRES':        
                 if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter  
                 p=GMRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,iter_restart=iter_restart)  
                 # solve Au=f-B^*p  
                 #       A(u-v)=f-B^*p-Av  
                 #       u=v+(u-v)  
         u=v+self.solve_A(v,p)  
   
           if solver=='TFQMR':        
                 if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter  
                 p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)  
                 # solve Au=f-B^*p  
                 #       A(u-v)=f-B^*p-Av  
                 #       u=v+(u-v)  
         u=v+self.solve_A(v,p)  
   
           if solver=='MINRES':        
                 if self.verbose: print "enter MINRES method (iter_max=%s)"%max_iter  
                 p=MINRES(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)  
                 # solve Au=f-B^*p  
                 #       A(u-v)=f-B^*p-Av  
                 #       u=v+(u-v)  
         u=v+self.solve_A(v,p)  
                 
           if solver=='GMRESC':        
                 if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter  
                 p0=self.solve_prec1(Bz)  
             #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))  
                 #p-=alfa  
                 x=GMRES(ArithmeticTuple(self.__z*1.,p0*1),self.__Anext,self.__Mempty,self.__inner_a,self.__stoppingcriterium2,iter_max=max_iter, x=ArithmeticTuple(v*1,p*1),iter_restart=20)  
                 #x=NewtonGMRES(ArithmeticTuple(self.__z*1.,p0*1),self.__Aprod_Newton2,self.__Mempty,self.__inner_a,self.__stoppingcriterium2,iter_max=max_iter, x=ArithmeticTuple(v*1,p*1),atol=0,rtol=self.getTolerance())  
   
                 # solve Au=f-B^*p  
                 #       A(u-v)=f-B^*p-Av  
                 #       u=v+(u-v)  
             p=x[1]  
         u=v+self.solve_A(v,p)        
         #p=x[1]  
         #u=x[0]  
   
               if solver=='PCG':  
                 #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv  
                 #  
                 #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)  
                 #                           A(v-z)= f -Az - B^*p (v-z=0 on fixed_u_mask)  
                 if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter  
                 p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)  
             u=r[0]    
                 # print "DIFF=",util.integrate(p)  
   
               # print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)  
1442    
1443            return u,p        def __inner_PCG(self,p,r):
1444             return self.inner_pBv(p,r[1])
1445    
1446        def __Msolve(self,r):        def __Msolve_PCG(self,r):
1447            return self.solve_prec(r[1])            return self.solve_prec(r[1])
1448          #=============================================================
1449          def __Aprod_GMRES(self,x):
1450              # return w,q  from v, p
1451              # solve Aw =Av+B^*p as Aw =f-A(z-v)-B^*(-p)
1452              #  and  Sq=B(v-w)
1453              v=x[0]
1454              p=x[1]
1455              w=self.solve_A(self.__z-v,-p)
1456              Bw=self.B(w-v)
1457          q=self.solve_prec(Bw)
1458              return ArithmeticTuple(w,q)
1459    
1460          def __inner_GMRES(self,a1,a2):
1461             return self.inner_p(a1[1],a2[1])+self.inner_v(a1[0],a2[0])
1462    
1463          #=============================================================
1464          def norm(self,v,p):
1465            f=self.inner_p(p,p)+self.inner_v(v,v)
1466            if f<0:
1467                raise ValueError,"negative norm."
1468            return math.sqrt(f)
1469    
1470        def __Msolve2(self,r):        def solve(self,v,p,max_iter=20, verbose=False, show_details=False, useUzawa=True, iter_restart=20,max_correction_steps=3):
1471            return self.solve_prec(r*1)           """
1472             solves the saddle point problem using initial guesses v and p.
1473    
1474        def __Mempty(self,r):           @param v: initial guess for velocity
1475            return r           @param p: initial guess for pressure
1476             @type v: L{Data}
1477             @type p: L{Data}
1478             @param useUzawa: indicate the usage of the Uzawa scheme. Otherwise the problem is solved in coupled form.
1479             @param max_iter: maximum number of iteration steps per correction attempt.
1480             @param verbose: show information on the progress of the saddlepoint problem solver.
1481             @param show_details: show details of the sub problems solves
1482             @param iter_restart: restart the iteration after C{iter_restart} steps (only used if useUzaw=False)
1483             @param max_correction_steps: maximum number of iteration steps in the attempt get |Bv| to zero.
1484             @return: new approximations for velocity and pressure
1485             @type useUzawa: C{bool}
1486             @type max_iter: C{int}
1487             @type verbose: C{bool}
1488             @type show_details: C{bool}
1489             @type iter_restart: C{int}
1490             @type max_correction_steps: C{int}
1491             @rtype: C{tuple} of L{Data} objects
1492             """
1493             self.verbose=verbose
1494             self.show_details=show_details and self.verbose
1495             #=================================================================================
1496             # Az=f is solved as A(z-v)=f-Av-B^*0 (typically with z-v=0 on boundary)
1497             self.setSubProblemTolerance(self.getTolerance(), apply_reduction=True)
1498             self.__z=v+self.solve_A(v,p*0)
1499             Bz=self.B(self.__z)
1500             p0=self.solve_prec(Bz)
1501             # tolerances:
1502             rtol=self.getTolerance()
1503             atol=self.getAbsoluteTolerance()
1504             f0=self.norm(self.__z,p0)
1505             if not f0>0:
1506                return self.__z, p*0
1507             ATOL=rtol*f0+atol
1508             if not ATOL>0:
1509                raise ValueError,"overall absolute tolerance needs to be positive."
1510             if self.verbose: print "saddle point solver: initial residual %e, tolerance = %e."%(f0,ATOL)
1511             # initialization
1512             self.iter=0
1513             correction_step=0
1514             converged=False
1515             # initial guess:
1516             q=p*1
1517             if not useUzawa:
1518               u=v
1519             while not converged :
1520                if useUzawa:
1521                   self.setSubProblemTolerance(ATOL/f0, apply_reduction=True)
1522                   # assume p is known: then v=z-A^-1B^*p
1523                   #      
1524                   # which leads to BA^-1B^*p = Bz
1525                   #
1526                   # note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bu
1527                   #
1528                   # note that q=BA^-1B^*p = B (A^-1B^*p)=Bw  with Aw=B^*p = f - Az - B^*(-p)
1529                   #
1530                   if self.verbose: print "enter PCG method (iter_max=%s, atol=%e, subtolerance=%e)"%(max_iter,ATOL, self.getSubProblemTolerance())
1531                   q,r=PCG(ArithmeticTuple(self.__z,Bz),self.__Aprod_PCG,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, x=q, verbose=self.verbose)
1532               u=r[0]  
1533               Bu=r[1]  
1534                else:
1535                   self.setSubProblemTolerance(ATOL/f0, apply_reduction=True)
1536                   #
1537                   #  with v=dv+z
1538                   #
1539                   #   A dv + B p = A (v-z) + Bp = 0
1540                   #   B dv       =  - Bz
1541                   #
1542                   # apply the preconditioner [[A^{-1} 0][(S^{-1} B A^{-1})  -S^{-1}]] to the problem
1543                   # the right hand side is then [0, S^{-1} Bz ]
1544                   #
1545                   #          
1546                   if self.verbose: print "enter GMRES (iter_max=%s, atol=%e, subtolerance=%e)"%(max_iter,ATOL,self.getSubProblemTolerance())
1547                   x=GMRES(ArithmeticTuple(0*self.__z,p0),self.__Aprod_GMRES,self.__inner_GMRES,atol=ATOL, rtol=0.,iter_max=max_iter, x=ArithmeticTuple(u-self.__z,q),iter_restart=iter_restart, verbose=self.verbose)
1548                   self.setSubProblemTolerance(ATOL/f0, apply_reduction=False)
1549               u=self.__z+x[0]
1550                   q=x[1]
1551               Bu=self.B(u)
1552                # now we check if |Bu| ~ 0 or more precise |Bu|  <= rtol * |v|
1553                nu=self.inner_v(u,u)
1554                nz=self.inner_v(self.__z,self.__z)
1555                p2=self.solve_prec(Bu)
1556                nBu=self.inner_p(p2,p2)
1557                if not nu>=0 and not nBu>=0: raise NegativeNorm,"negative norm."
1558                nu=math.sqrt(nu)
1559                nBu=math.sqrt(nBu)
1560                if self.verbose: print "saddle point solver: norm v= %e (Bv = %e)"%(nu,nBu)
1561                QTOL=atol+nu*rtol
1562                if nBu <= QTOL:
1563                    converged=True
1564                else:
1565                    ATOL=QTOL/nBu*ATOL*0.3
1566                    if self.verbose: print "correction step %s: tolerance reduced to %e."%(correction_step,ATOL)
1567                    converged=False
1568                correction_step+=1
1569                if correction_step>max_correction_steps:
1570                   raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1571             if self.verbose: print "saddle point solver: tolerance reached."
1572         return u,q
1573    
1574          #==========================================================================================================================
1575          def setTolerance(self,tolerance=1.e-4):
1576             """
1577             sets the relative tolerance for (v,p)
1578    
1579        def __Aprod(self,p):           @param tolerance: tolerance to be used
1580            # return BA^-1B*p           @type tolerance: non-negative C{float}
1581            #solve Av =B^*p as Av =f-Az-B^*(-p)           """
1582            v=self.solve_A(self.__z,-p)           if tolerance<0:
1583            return ArithmeticTuple(v, self.B(v))               raise ValueError,"tolerance must be positive."
1584             self.__rtol=tolerance
1585          def getTolerance(self):
1586             """
1587             returns the relative tolerance
1588    
1589        def __Anext(self,x):           @return: relative tolerance
1590            # return next v,p           @rtype: C{float}
1591            #solve Av =-B^*p as Av =f-Az-B^*p           """
1592             return self.__rtol
1593          def setAbsoluteTolerance(self,tolerance=0.):
1594             """
1595             sets the absolute tolerance
1596    
1597        pc=x[1]           @param tolerance: tolerance to be used
1598            v=self.solve_A(self.__z,-pc)           @type tolerance: non-negative C{float}
1599        p=self.solve_prec1(self.B(v))           """
1600             if tolerance<0:
1601                 raise ValueError,"tolerance must be non-negative."
1602             self.__atol=tolerance
1603          def getAbsoluteTolerance(self):
1604             """
1605             returns the absolute tolerance
1606    
1607            return ArithmeticTuple(v,p)           @return: absolute tolerance
1608             @rtype: C{float}
1609             """
1610             return self.__atol
1611          def setSubToleranceReductionFactor(self,reduction=None):
1612             """
1613             sets the reduction factor for the tolerance used to solve the subproblem in each iteration step.
1614             if set to C{None} the square of the problem tolerance is used.
1615        
1616             @param reduction: reduction factor.
1617             @type reduction: positive C{float} less or equal one or C{None}
1618             """
1619             if not reduction == None:
1620                  if reduction<=0 or reduction>1:
1621                      raise ValueError,"reduction factor (=%e) must be positive and less or equal one."%reduction
1622             self.__reduction=reduction
1623          def getSubToleranceReductionFactor(self):
1624             """
1625             returns the subproblem reduction factor
1626    
1627             @return: subproblem reduction factor
1628             @rtype: C{float} or C{None}
1629             """
1630             return self.__reduction
1631    
1632        def __Aprod2(self,p):        def setSubProblemTolerance(self,rtol,apply_reduction=False):
1633            # return BA^-1B*p           """
1634            #solve Av =B^*p as Av =f-Az-B^*(-p)           sets the relative tolerance to solve the subproblem(s).
       v=self.solve_A(self.__z,-p)  
           return self.B(v)  
1635    
1636        def __Aprod_Newton(self,p):           @param rtol: relative tolerence
1637            # return BA^-1B*p - Bz           @type rtol: positive C{float}
1638            #solve Av =-B^*p as Av =f-Az-B^*p           @param apply_reduction: if set the the reduction factor is applied.
1639        v=self.solve_A(self.__z,-p)           @type apply_reduction: C{bool}
1640            return self.B(v-self.__z)           """
1641             if rtol<=0:
1642        def __Aprod_Newton2(self,x):               raise ValueError,"tolerance must be positive."
1643            # return BA^-1B*p - Bz           if apply_reduction:
1644            #solve Av =-B^*p as Av =f-Az-B^*p               if self.__reduction == None:
1645            pc=x[1]                   self.__sub_tol=max(rtol**2,util.EPSILON*10)
1646        v=self.solve_A(self.__z,-pc)               else:
1647            p=self.solve_prec1(self.B(v-self.__z))                   self.__sub_tol=max(self.__reduction*rtol,util.EPSILON*10)
1648            return ArithmeticTuple(v,p)           else:
1649                 self.__sub_tol=rtol
1650          def getSubProblemTolerance(self):
1651             """
1652             returns the subproblem reduction factor
1653    
1654             @return: subproblem reduction factor
1655             @rtype: C{float}
1656             """
1657             return self.__sub_tol
1658    
1659  def MaskFromBoundaryTag(domain,*tags):  def MaskFromBoundaryTag(domain,*tags):
1660     """     """
1661     creates a mask on the Solution(domain) function space which one for samples     create a mask on the Solution(domain) function space which one for samples
1662     that touch the boundary tagged by tags.     that touch the boundary tagged by tags.
1663    
1664     usage: m=MaskFromBoundaryTag(domain,"left", "right")     usage: m=MaskFromBoundaryTag(domain,"left", "right")

Legend:
Removed from v.2061  
changed lines
  Added in v.2100

  ViewVC Help
Powered by ViewVC 1.1.26