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

revision 1877 by ksteube, Thu Sep 25 06:43:44 2008 UTC revision 1878 by gross, Tue Oct 14 03:39:13 2008 UTC
# Line 560  type like argument C{x}. Line 560  type like argument C{x}.
560
561     return x,r     return x,r
562
563    class Defect(object):
564        """
565        defines a non-linear defect F(x) of a variable x
566        """
567        def __init__(self):
568            """
569            initialize defect
570            """
571            self.setDerivativeIncrementLength()
572
573  ############################      def bilinearform(self, x0, x1):
574  # Added by Artak          """
575  #################################3          returns the inner product of x0 and x1
576            @param x0: a value for x
577            @param x1: a value for x
578            @return: the inner product of x0 and x1
579            @rtype: C{float}
580            """
581            return 0
582
583        def norm(self,x):
584            """
585            the norm of argument C{x}
586
587            @param x: a value for x
588            @return: norm of argument x
589            @rtype: C{float}
590            @note: by default C{sqrt(self.bilinearform(x,x)} is retrurned.
591            """
592            s=self.bilinearform(x,x)
593            if s<0: raise NegativeNorm,"negative norm."
594            return math.sqrt(s)
595
596
597        def eval(self,x):
598            """
599            returns the value F of a given x
600
601            @param x: value for which the defect C{F} is evalulated.
602            @return: value of the defect at C{x}
603            """
604            return 0
605
606        def __call__(self,x):
607            return self.eval(x)
608
609        def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):
610            """
611            sets the relative length of the increment used to approximate the derivative of the defect
612            the increment is inc*norm(x)/norm(v)*v in the direction of v with x as a staring point.
613
614            @param inc: relative increment length
615            @type inc: positive C{float}
616            """
617            if inc<=0: raise ValueError,"positive increment required."
618            self.__inc=inc
619
620        def getDerivativeIncrementLength(self):
621            """
622            returns the relative increment length used to approximate the derivative of the defect
623            @return: value of the defect at C{x}
624            @rtype: positive C{float}
625            """
626            return self.__inc
627
628        def derivative(self, F0, x0, v, v_is_normalised=True):
629            """
630            returns the directional derivative at x0 in the direction of v
631
632            @param F0: value of this defect at x0
633            @param x0: value at which derivative is calculated.
634            @param v: direction
635            @param v_is_normalised: is true to indicate that C{v} is nomalized (self.norm(v)=0)
636            @return: derivative of this defect at x0 in the direction of C{v}
637            @note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is used but this method
638            maybe oepsnew verwritten to use exact evalution.
639            """
640            normx=self.norm(x0)
641            if normx>0:
642                 epsnew = self.getDerivativeIncrementLength() * normx
643            else:
644                 epsnew = self.getDerivativeIncrementLength()
645            if not v_is_normalised:
646                normv=self.norm(v)
647                if normv<=0:
648                   return F0*0
649                else:
650                   epsnew /= normv
651            F1=self.eval(x0 + epsnew * v)
652            return (F1-F0)/epsnew
653
654    ######################################
655    def NewtonGMRES(defect, x, iter_max=100, sub_iter_max=20, atol=0,rtol=1.e-4, sub_tol_max=0.5, gamma=0.9, verbose=False):
656       """
657       solves a non-linear problem M{F(x)=0} for unknown M{x} using the stopping criterion:
658
659       M{norm(F(x) <= atol + rtol * norm(F(x0)}
660
661       where M{x0} is the initial guess.
662
663  #Apply a sequence of k Givens rotations, used within gmres codes     @param defect: object defining the the function M{F}, C{defect.norm} defines the M{norm} used in the stopping criterion.
664  # vrot=givapp(c, s, vin, k)     @type defect: L{Defect}
665  def givapp(c,s,vin):     @param x: initial guess for the solution, C{x} is altered.
666      vrot=vin # warning: vin is altered!!!!     @type x: any object type allowing basic operations such as  L{numarray.NumArray}, L{Data}
667       @param iter_max: maximum number of iteration steps
668       @type iter_max: positive C{int}
669       @param sub_iter_max:
670       @type sub_iter_max:
671       @param atol: absolute tolerance for the solution
672       @type atol: positive C{float}
673       @param rtol: relative tolerance for the solution
674       @type rtol: positive C{float}
675       @param gamma: tolerance safety factor for inner iteration
676       @type gamma: positive C{float}, less than 1
677       @param sub_tol_max: upper bound for inner tolerance.
678       @type sub_tol_max: positive C{float}, less than 1
679       @return: an approximation of the solution with the desired accuracy
680       @rtype: same type as the initial guess.
681       """
682       lmaxit=iter_max
683       if atol<0: raise ValueError,"atol needs to be non-negative."
684       if rtol<0: raise ValueError,"rtol needs to be non-negative."
685       if rtol+atol<=0: raise ValueError,"rtol or atol needs to be non-negative."
686       if gamma<=0 or gamma>=1: raise ValueError,"tolerance safety factor for inner iteration (gamma =%s) needs to be positive and less than 1."%gamma
687       if sub_tol_max<=0 or sub_tol_max>=1: raise ValueError,"upper bound for inner tolerance for inner iteration (sub_tol_max =%s) needs to be positive and less than 1."%sub_tol_max
688
689       F=defect(x)
690       fnrm=defect.norm(F)
691       stop_tol=atol + rtol*fnrm
692       sub_tol=sub_tol_max
693       if verbose: print "NewtonGMRES: initial residual = %e."%fnrm
694       if verbose: print "             tolerance = %e."%sub_tol
695       iter=1
696       #
697       # main iteration loop
698       #
699       while not fnrm<=stop_tol:
700                if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
701                #
702            #   adjust sub_tol_
703            #
704                if iter > 1:
705               rat=fnrm/fnrmo
706                   sub_tol_old=sub_tol
707               sub_tol=gamma*rat**2
708               if gamma*sub_tol_old**2 > .1: sub_tol=max(sub_tol,gamma*sub_tol_old**2)
709               sub_tol=max(min(sub_tol,sub_tol_max), .5*stop_tol/fnrm)
710            #
711            # calculate newton increment xc
712                #     if iter_max in __FDGMRES is reached MaxIterReached is thrown
713                #     if iter_restart -1 is returned as sub_iter
714                #     if  atol is reached sub_iter returns the numer of steps performed to get there
715                #
716                #
717                if verbose: print "             subiteration (GMRES) is called with relative tolerance %e."%sub_tol
718                try:
719                   xc, sub_iter=__FDGMRES(F, defect, x, sub_tol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)
720                except MaxIterReached:
721                   raise MaxIterReached,"maximum number of %s steps reached."%iter_max
722                if sub_iter<0:
723                   iter+=sub_iter_max
724                else:
725                   iter+=sub_iter
726                # ====
727            x+=xc
728                F=defect(x)
729            iter+=1
730                fnrmo, fnrm=fnrm, defect.norm(F)
731                if verbose: print "             step %s: residual %e."%(iter,fnrm)
732       if verbose: print "NewtonGMRES: completed after %s steps."%iter
733       return x
734
735    def __givapp(c,s,vin):
736        """
737        apply a sequence of Givens rotations (c,s) to the recuirsively to the vector vin
738        @warning: C{vin} is altered.
739        """
740        vrot=vin
741      if isinstance(c,float):      if isinstance(c,float):
742          vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]          vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
743      else:      else:
# Line 578  def givapp(c,s,vin): Line 747  def givapp(c,s,vin):
747              vrot[i:i+2]=w1,w2              vrot[i:i+2]=w1,w2
748      return vrot      return vrot
749
750  ##############################################  def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
################################################
m=iter_restart
iter=0
xc=x
while True:
if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
xc,stopped=GMRESm(b*1, Aprod, Msolve, bilinearform, stoppingcriterium, x=xc*1, iter_max=iter_max-iter, iter_restart=m)
if stopped: break
iter+=iter_restart
return xc

def GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
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."

751     h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)     h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
752     c=numarray.zeros(iter_restart,numarray.Float64)     c=numarray.zeros(iter_restart,numarray.Float64)
753     s=numarray.zeros(iter_restart,numarray.Float64)     s=numarray.zeros(iter_restart,numarray.Float64)
754     g=numarray.zeros(iter_restart,numarray.Float64)     g=numarray.zeros(iter_restart,numarray.Float64)
755     v=[]     v=[]
756
757     rho=math.sqrt(r_dot_r)     rho=defect.norm(F0)
758       if rho<=0.: return x0*0
759
760     v.append(r/rho)     v.append(-F0/rho)
761     g[0]=rho     g[0]=rho
762       iter=0
763     while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):     while rho > atol and iter<iter_restart-1:
764
765      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
766
767      p=Msolve(Aprod(v[iter]))          p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)

768      v.append(p)      v.append(p)
769
770      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))        v_norm1=defect.norm(v[iter+1])
771
772  # Modified Gram-Schmidt          # Modified Gram-Schmidt
773      for j in range(iter+1):      for j in range(iter+1):
774        h[j][iter]=bilinearform(v[j],v[iter+1])             h[j][iter]=defect.bilinearform(v[j],v[iter+1])
775        v[iter+1]-=h[j][iter]*v[j]           v[iter+1]-=h[j][iter]*v[j]
776
777      h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))      h[iter+1][iter]=defect.norm(v[iter+1])
778      v_norm2=h[iter+1][iter]      v_norm2=h[iter+1][iter]
779
780  # Reorthogonalize if needed          # Reorthogonalize if needed
781      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)
782       for j in range(iter+1):            for j in range(iter+1):
783          hr=bilinearform(v[j],v[iter+1])             hr=defect.bilinearform(v[j],v[iter+1])
784              h[j][iter]=h[j][iter]+hr                 h[j][iter]=h[j][iter]+hr
785              v[iter+1] -= hr*v[j]                 v[iter+1] -= hr*v[j]
786
787       v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))            v_norm2=defect.norm(v[iter+1])
788       h[iter+1][iter]=v_norm2          h[iter+1][iter]=v_norm2
789            #   watch out for happy breakdown
#   watch out for happy breakdown
790          if not v_norm2 == 0:          if not v_norm2 == 0:
791           v[iter+1]=v[iter+1]/h[iter+1][iter]                  v[iter+1]=v[iter+1]/h[iter+1][iter]
792
793  #   Form and store the information for the new Givens rotation          #   Form and store the information for the new Givens rotation
794      if iter > 0 :      if iter > 0 :
795          hhat=numarray.zeros(iter+1,numarray.Float64)          hhat=numarray.zeros(iter+1,numarray.Float64)
796          for i in range(iter+1) : hhat[i]=h[i][iter]          for i in range(iter+1) : hhat[i]=h[i][iter]
797          hhat=givapp(c[0:iter],s[0:iter],hhat);          hhat=__givapp(c[0:iter],s[0:iter],hhat);
798              for i in range(iter+1) : h[i][iter]=hhat[i]              for i in range(iter+1) : h[i][iter]=hhat[i]
799
800      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])
# Line 662  def GMRESm(b, Aprod, Msolve, bilinearfor Line 804  def GMRESm(b, Aprod, Msolve, bilinearfor
804          s[iter]=-h[iter+1][iter]/mu          s[iter]=-h[iter+1][iter]/mu
805          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]
806          h[iter+1][iter]=0.0          h[iter+1][iter]=0.0
807          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])
808
809  # Update the residual norm          # Update the residual norm

810          rho=abs(g[iter+1])          rho=abs(g[iter+1])
811      iter+=1      iter+=1
812
813  # At this point either iter > iter_max or rho < tol.     # At this point either iter > iter_max or rho < tol.
814  # It's time to compute x and leave.             # It's time to compute x and leave.

815     if iter > 0 :     if iter > 0 :
816       y=numarray.zeros(iter,numarray.Float64)           y=numarray.zeros(iter,numarray.Float64)
817       y[iter-1] = g[iter-1] / h[iter-1][iter-1]       y[iter-1] = g[iter-1] / h[iter-1][iter-1]
# Line 683  def GMRESm(b, Aprod, Msolve, bilinearfor Line 823  def GMRESm(b, Aprod, Msolve, bilinearfor
823       xhat=v[iter-1]*y[iter-1]       xhat=v[iter-1]*y[iter-1]
824       for i in range(iter-1):       for i in range(iter-1):
825      xhat += v[i]*y[i]      xhat += v[i]*y[i]
826     else : xhat=v[0]     else :
827          xhat=v[0] * 0
828
x += xhat
829     if iter<iter_restart-1:     if iter<iter_restart-1:
830        stopped=True        stopped=iter
831     else:     else:
832        stopped=False        stopped=-1

return x,stopped

######################################################
def dirder(x, w, bilinearform, Aprod, Msolve, f0, b ):
######################################################

# DIRDER estimates the directional derivative of a function F.

# Hardwired difference increment.
#
epsnew = 1.0e-07
#
#  Scale the step.
#
norm_w=math.sqrt(bilinearform(w,w))
if norm_w== 0.0:
return x/x

epsnew = epsnew / norm_w

if norm_w > 0.0:
epsnew = epsnew * math.sqrt(bilinearform(x,x))
#
#  DEL and F1 could share the same space if storage
#  is more important than clarity.
#
833
834    DEL = x + epsnew * w     return xhat,stopped
f1 = -Msolve(Aprod(DEL))
z = ( f1 - f0 ) / epsnew
return z

######################################################
def FDGMRES(f0, Aprod, Msolve, bilinearform, stoppingcriterium, xc=None, x=None, iter_max=100, iter_restart=20,TOL=None):
######################################################
b=-f0
b_dot_b = bilinearform(b, b)
if b_dot_b<0: raise NegativeNorm,"negative norm."
norm_b=math.sqrt(b_dot_b)
835
836     r=b  ##############################################
837    def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
838    ################################################
839       m=iter_restart
840       iter=0
841       xc=x
842       while True:
843          if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
844          xc,stopped=__GMRESm(b*1, Aprod, Msolve, bilinearform, stoppingcriterium, x=xc*1, iter_max=iter_max-iter, iter_restart=m)
845          if stopped: break
846          iter+=iter_restart
847       return xc
848
849     if x==None:  def __GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
850        x=0*f0     iter=0
851     else:     r=Msolve(b)
r=-dirder(xc,x,bilinearform,Aprod,Msolve,f0,b)-f0

852     r_dot_r = bilinearform(r, r)     r_dot_r = bilinearform(r, r)
853     if r_dot_r<0: raise NegativeNorm,"negative norm."     if r_dot_r<0: raise NegativeNorm,"negative norm."
854       norm_b=math.sqrt(r_dot_r)
855
856       if x==None:
857          x=0*b
858       else:
859          r=Msolve(b-Aprod(x))
860          r_dot_r = bilinearform(r, r)
861          if r_dot_r<0: raise NegativeNorm,"negative norm."
862
863     h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)     h=numarray.zeros((iter_restart,iter_restart),numarray.Float64)
864     c=numarray.zeros(iter_restart,numarray.Float64)     c=numarray.zeros(iter_restart,numarray.Float64)
# Line 753  def FDGMRES(f0, Aprod, Msolve, bilinearf Line 870  def FDGMRES(f0, Aprod, Msolve, bilinearf
870
871     v.append(r/rho)     v.append(r/rho)
872     g[0]=rho     g[0]=rho
iter=0
873
874     while not (stoppingcriterium(rho,norm_b,solver="FDGMRES",TOL=TOL) or iter==iter_restart-1):     while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
875
876      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
877
878            p=Msolve(Aprod(v[iter]))
p=dirder(xc, v[iter], bilinearform,Aprod,Msolve,f0,b)
879
880      v.append(p)      v.append(p)
881
882      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
883
884  # Modified Gram-Schmidt  # Modified Gram-Schmidt
885      for j in range(iter+1):      for j in range(iter+1):
886        h[j][iter]=bilinearform(v[j],v[iter+1])          h[j][iter]=bilinearform(v[j],v[iter+1])
887        v[iter+1]+=(-1.)*h[j][iter]*v[j]        v[iter+1]-=h[j][iter]*v[j]
888
889      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]))
890      v_norm2=h[iter+1][iter]      v_norm2=h[iter+1][iter]
891

892  # Reorthogonalize if needed  # Reorthogonalize if needed
893      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)
894       for j in range(iter+1):       for j in range(iter+1):
895          hr=bilinearform(v[j],v[iter+1])          hr=bilinearform(v[j],v[iter+1])
896              h[j][iter]=h[j][iter]+hr #vhat              h[j][iter]=h[j][iter]+hr
897              v[iter+1] +=(-1.)*hr*v[j]              v[iter+1] -= hr*v[j]
898
899       v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
900       h[iter+1][iter]=v_norm2       h[iter+1][iter]=v_norm2
901
902  #   watch out for happy breakdown  #   watch out for happy breakdown
903          if v_norm2 != 0:          if not v_norm2 == 0:
904           v[iter+1]=v[iter+1]/h[iter+1][iter]           v[iter+1]=v[iter+1]/h[iter+1][iter]
905
906  #   Form and store the information for the new Givens rotation  #   Form and store the information for the new Givens rotation
907      if iter > 0 :      if iter > 0 :
908          hhat=[]          hhat=numarray.zeros(iter+1,numarray.Float64)
909          for i in range(iter+1) : hhat.append(h[i][iter])          for i in range(iter+1) : hhat[i]=h[i][iter]
910          hhat=givapp(c[0:iter],s[0:iter],hhat);          hhat=__givapp(c[0:iter],s[0:iter],hhat);
911              for i in range(iter+1) : h[i][iter]=hhat[i]              for i in range(iter+1) : h[i][iter]=hhat[i]
912
913      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])
914
915      if mu!=0 :      if mu!=0 :
916          c[iter]=h[iter][iter]/mu          c[iter]=h[iter][iter]/mu
917          s[iter]=-h[iter+1][iter]/mu          s[iter]=-h[iter+1][iter]/mu
918          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]
919          h[iter+1][iter]=0.0          h[iter+1][iter]=0.0
920          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])
921
922  # Update the residual norm  # Update the residual norm
923
924          rho=abs(g[iter+1])          rho=abs(g[iter+1])
925      iter+=1      iter+=1
926
927    # At this point either iter > iter_max or rho < tol.
928    # It's time to compute x and leave.
929
930     if iter > 0 :     if iter > 0 :
931       y=numarray.zeros(iter,numarray.Float64)           y=numarray.zeros(iter,numarray.Float64)
932       y[iter-1] = g[iter-1] / h[iter-1][iter-1]       y[iter-1] = g[iter-1] / h[iter-1][iter-1]
# Line 820  def FDGMRES(f0, Aprod, Msolve, bilinearf Line 939  def FDGMRES(f0, Aprod, Msolve, bilinearf
939       for i in range(iter-1):       for i in range(iter-1):
940      xhat += v[i]*y[i]      xhat += v[i]*y[i]
941     else : xhat=v[0]     else : xhat=v[0]
942
943     x += xhat     x += xhat
944     if iter<iter_restart-1:     if iter<iter_restart-1:
945        stopped=True        stopped=True
# Line 983  def MINRES(b, Aprod, Msolve, bilinearfor Line 1102  def MINRES(b, Aprod, Msolve, bilinearfor
1102
1103      return x      return x
1104
######################################
def NewtonGMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium,x=None, iter_max=100,iter_restart=20,atol=1.e-2,rtol=1.e-4):
#####################################
gamma=.9
lmaxit=100
etamax=.5

n = 1 #len(x)
iter=0

# evaluate f at the initial iterate
# compute the stop tolerance
#
r=b
if x==None:
x=0*b
else:
r += (-1)*Aprod(x)

f0=-Msolve(r)
fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
fnrmo=1
stop_tol=atol + rtol*fnrm
#
# main iteration loop
#
while not stoppingcriterium(fnrm*1,stop_tol,'NewtonGMRES',TOL=1.):

if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
#
# keep track of the ratio (rat = fnrm/frnmo)
# of successive residual norms and
# the iteration counter (iter)
#
#rat=fnrm/fnrmo
fnrmo=fnrm
iter+=1
#
# compute the step using a GMRES(m) routine especially designed for this purpose
#
initer=0
while True:
xc,stopped=FDGMRES(f0*1, Aprod, Msolve, bilinearform, stoppingcriterium, xc=x*1, iter_max=lmaxit-initer, iter_restart=iter_restart, TOL=etamax)
if stopped: break
initer+=iter_restart
x+=xc
f0=-Msolve(Aprod(x))
fnrm=math.sqrt(bilinearform(f0,f0))/math.sqrt(n)
rat=fnrm/fnrmo

#   adjust eta
#
if etamax > 0:
etaold=etamax
etanew=gamma*rat*rat
if gamma*etaold*etaold > .1 :
etanew=max(etanew,gamma*etaold*etaold)
etamax=min(etanew,etamax)
etamax=max(etamax,.5*stop_tol/fnrm)
return x

1105  def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
1106
1107  # TFQMR solver for linear systems  # TFQMR solver for linear systems
# Line 1340  class ArithmeticTuple(object): Line 1397  class ArithmeticTuple(object):
1397
1398  class HomogeneousSaddlePointProblem(object):  class HomogeneousSaddlePointProblem(object):
1399        """        """
1400        This provides a framwork for solving homogeneous saddle point problem of the form        This provides a framwork for solving linear homogeneous saddle point problem of the form
1401
1402               Av+B^*p=f               Av+B^*p=f
1403               Bv    =0               Bv    =0
# Line 1460  class HomogeneousSaddlePointProblem(obje Line 1517  class HomogeneousSaddlePointProblem(obje
1517                  #       u=v+(u-v)                  #       u=v+(u-v)
1518          u=v+self.solve_A(v,p)          u=v+self.solve_A(v,p)
1519
if solver=='NewtonGMRES':
if self.verbose: print "enter NewtonGMRES method (iter_max=%s)"%max_iter
p=NewtonGMRES(Bz,self.__Aprod_Newton,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.,atol=0,rtol=self.getTolerance())
# solve Au=f-B^*p
#       A(u-v)=f-B^*p-Av
#       u=v+(u-v)
u=v+self.solve_A(v,p)

1520            if solver=='TFQMR':                  if solver=='TFQMR':
1521                  if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter                  if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter
1522                  p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)                  p=TFQMR(Bz,self.__Aprod2,self.__Msolve2,self.__inner_p,self.__stoppingcriterium2,iter_max=max_iter, x=p*1.)
# Line 1562  class HomogeneousSaddlePointProblem(obje Line 1610  class HomogeneousSaddlePointProblem(obje
1610            p=self.solve_prec1(self.B(v-self.__z))            p=self.solve_prec1(self.B(v-self.__z))
1611            return ArithmeticTuple(v,p)            return ArithmeticTuple(v,p)
1612
1613
1614    def MaskFromBoundaryTag(domain,*tags):
1615       """
1616       create a mask on the Solution(domain) function space which one for samples
1617       that touch the boundary tagged by tags.
1618
1619       usage: m=MaskFromBoundaryTag(domain,"left", "right")
1620
1621       @param domain: a given domain
1622       @type domain: L{escript.Domain}
1623       @param tags: boundray tags
1624       @type tags: C{str}
1625       @return: a mask which marks samples that are touching the boundary tagged by any of the given tags.
1626       @rtype: L{escript.Data} of rank 0
1627       """
1628       pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1629       d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
1630       for t in tags: d.setTaggedValue(t,1.)
1631       pde.setValue(y=d)
1632       return util.whereNonZero(pde.getRightHandSide())
1633    #============================================================================================================================================
1634  class SaddlePointProblem(object):  class SaddlePointProblem(object):
1635     """     """
1636     This implements a solver for a saddlepoint problem     This implements a solver for a saddlepoint problem
# Line 1592  class SaddlePointProblem(object): Line 1661  class SaddlePointProblem(object):
1661              raise TypeError("verbose needs to be of type bool.")              raise TypeError("verbose needs to be of type bool.")
1662         self.__verbose=verbose         self.__verbose=verbose
1663         self.relaxation=1.         self.relaxation=1.
1664           DeprecationWarning("SaddlePointProblem should not be used anymore.")
1665
1666     def trace(self,text):     def trace(self,text):
1667         """         """
# Line 1765  class SaddlePointProblem(object): Line 1835  class SaddlePointProblem(object):
1835              f, norm_f, u, norm_u, g, norm_g, p, norm_p = f_new, norm_f_new, u_new, norm_u_new, g_new, norm_g_new, p_new, norm_p_new              f, norm_f, u, norm_u, g, norm_g, p, norm_p = f_new, norm_f_new, u_new, norm_u_new, g_new, norm_g_new, p_new, norm_p_new
1836          self.trace("convergence after %s steps."%self.iter)          self.trace("convergence after %s steps."%self.iter)
1837          return u,p          return u,p
#   def solve(self,u0,p0,tolerance=1.e-6,iter_max=10,self.relaxation=1.):
#      tol=1.e-2
#      iter=0
#      converged=False
#      u=u0*1.
#      p=p0*1.
#      while not converged and iter<iter_max:
#          du=self.solve_f(u,p,tol)
#          u-=du
#          norm_du=util.Lsup(du)
#          norm_u=util.Lsup(u)
#
#          dp=self.relaxation*self.solve_g(u,tol)
#          p+=dp
#          norm_dp=util.sqrt(self.inner(dp,dp))
#          norm_p=util.sqrt(self.inner(p,p))
#          print iter,"-th step rel. errror u,p= (%s,%s) (%s,%s)(%s,%s)"%(norm_du,norm_dp,norm_du/norm_u,norm_dp/norm_p,norm_u,norm_p)
#          iter+=1
#
#          converged = (norm_du <= tolerance*norm_u) and  (norm_dp <= tolerance*norm_p)
#      if converged:
#          print "convergence after %s steps."%iter
#      else:
#          raise ArithmeticError("no convergence after %s steps."%iter)
#
#      return u,p

def MaskFromBoundaryTag(domain,*tags):
"""
create a mask on the Solution(domain) function space which one for samples
that touch the boundary tagged by tags.

usage: m=MaskFromBoundaryTag(domain,"left", "right")

@param domain: a given domain
@type domain: L{escript.Domain}
@param tags: boundray tags
@type tags: C{str}
@return: a mask which marks samples that are touching the boundary tagged by any of the given tags.
@rtype: L{escript.Data} of rank 0
"""
pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
d=escript.Scalar(0.,escript.FunctionOnBoundary(function_space.getDomain()))
for t in tags: d.setTaggedValue(t,1.)
pde.setValue(y=d)
return util.whereNonZero(pde.getRightHandSide())

Legend:
 Removed from v.1877 changed lines Added in v.1878

 ViewVC Help Powered by ViewVC 1.1.26