/[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 1542 by artak, Wed Apr 30 03:40:15 2008 UTC revision 1956 by gross, Mon Nov 3 05:08:42 2008 UTC
# Line 1  Line 1 
1    
2    ########################################################
3  #  #
4  # $Id$  # Copyright (c) 2003-2008 by University of Queensland
5  #  # Earth Systems Science Computational Center (ESSCC)
6  #######################################################  # http://www.uq.edu.au/esscc
7  #  #
8  #           Copyright 2003-2007 by ACceSS MNRF  # Primary Business: Queensland, Australia
9  #       Copyright 2007 by University of Queensland  # Licensed under the Open Software License version 3.0
10  #  # http://www.opensource.org/licenses/osl-3.0.php
11  #                http://esscc.uq.edu.au  #
12  #        Primary Business: Queensland, Australia  ########################################################
13  #  Licensed under the Open Software License version 3.0  
14  #     http://www.opensource.org/licenses/osl-3.0.php  __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15  #  Earth Systems Science Computational Center (ESSCC)
16  #######################################################  http://www.uq.edu.au/esscc
17  #  Primary Business: Queensland, Australia"""
18    __license__="""Licensed under the Open Software License version 3.0
19    http://www.opensource.org/licenses/osl-3.0.php"""
20    __url__="http://www.uq.edu.au/esscc/escript-finley"
21    
22  """  """
23  Provides some tools related to PDEs.  Provides some tools related to PDEs.
# Line 32  Currently includes: Line 37  Currently includes:
37  """  """
38    
39  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
 __copyright__="""  Copyright (c) 2006 by ACcESS MNRF  
                     http://www.access.edu.au  
                 Primary Business: Queensland, Australia"""  
 __license__="""Licensed under the Open Software License version 3.0  
              http://www.opensource.org/licenses/osl-3.0.php"""  
 __url__="http://www.iservo.edu.au/esys"  
 __version__="$Revision$"  
 __date__="$Date$"  
40    
41    
42  import escript  import escript
# Line 474  class IterationHistory(object): Line 471  class IterationHistory(object):
471    
472         """         """
473         self.history.append(norm_r)         self.history.append(norm_r)
474         if self.verbose: print "iter: %s:  inner(rhat,r) = %e"%(len(self.history)-1, self.history[-1])         if self.verbose: print "iter: %s:  inner(rhat,r) = %e"#(len(self.history)-1, self.history[-1])
475         return self.history[-1]<=self.tolerance * self.history[0]         return self.history[-1]<=self.tolerance * self.history[0]
476    
477     def stoppingcriterium2(self,norm_r,norm_b,solver="GMRES",TOL=None):     def stoppingcriterium2(self,norm_r,norm_b,solver="GMRES",TOL=None):
# Line 493  class IterationHistory(object): Line 490  class IterationHistory(object):
490         if TOL==None:         if TOL==None:
491            TOL=self.tolerance            TOL=self.tolerance
492         self.history.append(norm_r)         self.history.append(norm_r)
493         if self.verbose: print "iter: %s:  norm(r) = %e"%(len(self.history)-1, self.history[-1])         if self.verbose: print "iter: %s:  norm(r) = %e"#(len(self.history)-1, self.history[-1])
494         return self.history[-1]<=TOL * norm_b         return self.history[-1]<=TOL * norm_b
495    
496  def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  def PCG(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
# Line 563  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            """
575            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  # Added by Artak          @param x0: value at which derivative is calculated.
634  #################################3          @param v: direction
635            @param v_is_normalised: is true to indicate that C{v} is nomalized (self.norm(v)=0)
636  #Apply a sequence of k Givens rotations, used within gmres codes          @return: derivative of this defect at x0 in the direction of C{v}
637  # vrot=givapp(c, s, vin, k)          @note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is used but this method
638  def givapp(c,s,vin):          maybe oepsnew verwritten to use exact evalution.
639      vrot=vin # warning: vin is altered!!!!          """
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       @param defect: object defining the the function M{F}, C{defect.norm} defines the M{norm} used in the stopping criterion.
664       @type defect: L{Defect}
665       @param x: initial guess for the solution, C{x} is altered.
666       @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 581  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 GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):  def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
    m=iter_restart  
    iter=0  
    while True:  
       if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max  
       x,stopped=GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=x, iter_max=iter_max-iter, iter_restart=m)  
       if stopped: break  
       iter+=iter_restart      
    return x  
   
 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  
    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     while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):     iter=0
763       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=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
     p=Msolve(Aprod(v[iter]))  
   
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]+=(-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 #vhat                 h[j][iter]=h[j][iter]+hr
785              v[iter+1] +=(-1.)*hr*v[j]                 v[iter+1] -= hr*v[j]
786    
787            v_norm2=defect.norm(v[iter+1])
788            h[iter+1][iter]=v_norm2
789            #   watch out for happy breakdown
790            if not v_norm2 == 0:
791                    v[iter+1]=v[iter+1]/h[iter+1][iter]
792    
793       v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))            #   Form and store the information for the new Givens rotation
      h[iter+1][iter]=v_norm2  
   
 #   watch out for happy breakdown  
         if v_norm2 != 0:  
          v[iter+1]=v[iter+1]/h[iter+1][iter]  
   
 #   Form and store the information for the new Givens rotation  
794      if iter > 0 :      if iter > 0 :
795          hhat=[]          hhat=numarray.zeros(iter+1,numarray.Float64)
796          for i in range(iter+1) : hhat.append(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])
801    
802      if mu!=0 :      if mu!=0 :
803          c[iter]=h[iter][iter]/mu          c[iter]=h[iter][iter]/mu
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 682  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.  
833    
834       return xhat,stopped
835    
836  # Hardwired difference increment.  ##############################################
837  #  def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
838    epsnew = 1.0e-07  ################################################
839  #     m=iter_restart
840  #  Scale the step.     iter=0
841  #     xc=x
842    norm_w=math.sqrt(bilinearform(w,w))     while True:
843    if norm_w== 0.0:        if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached"%iter_max
844      return x/x        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    epsnew = epsnew / norm_w        iter+=iter_restart    
847       return xc
   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.  
 #  
   
   DEL = x + epsnew * w  
   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)  
   
    r=b  
848    
849     if x==None:  def __GMRESm(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100, iter_restart=20):
850        x=0     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 751  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 818  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 827  def FDGMRES(f0, Aprod, Msolve, bilinearf Line 948  def FDGMRES(f0, Aprod, Msolve, bilinearf
948    
949     return x,stopped     return x,stopped
950    
951    #################################################
952  def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  def MINRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
953    #################################################
954      #      #
955      #  minres solves the system of linear equations Ax = b      #  minres solves the system of linear equations Ax = b
956      #  where A is a symmetric matrix (possibly indefinite or singular)      #  where A is a symmetric matrix (possibly indefinite or singular)
# Line 858  def MINRES(b, Aprod, Msolve, bilinearfor Line 980  def MINRES(b, Aprod, Msolve, bilinearfor
980    
981      r1    = b      r1    = b
982      y = Msolve(b)      y = Msolve(b)
983      beta1 = bilinearform(b,y)      beta1 = bilinearform(y,b)
984    
985      if beta1< 0: raise NegativeNorm,"negative norm."      if beta1< 0: raise NegativeNorm,"negative norm."
986    
# Line 920  def MINRES(b, Aprod, Msolve, bilinearfor Line 1042  def MINRES(b, Aprod, Msolve, bilinearfor
1042            y = y - (beta/oldb)*r1            y = y - (beta/oldb)*r1
1043    
1044          alfa   = bilinearform(v,y)              # alphak          alfa   = bilinearform(v,y)              # alphak
1045          y      = (- alfa/beta)*r2 + y          y      += (- alfa/beta)*r2
1046          r1     = r2          r1     = r2
1047          r2     = y          r2     = y
1048          y = Msolve(r2)          y = Msolve(r2)
1049          oldb   = beta                           # oldb = betak          oldb   = beta                           # oldb = betak
1050          beta   = bilinearform(r2,y)             # beta = betak+1^2          beta   = bilinearform(y,r2)             # beta = betak+1^2
1051          if beta < 0: raise NegativeNorm,"negative norm."          if beta < 0: raise NegativeNorm,"negative norm."
1052    
1053          beta   = math.sqrt( beta )          beta   = math.sqrt( beta )
# Line 960  def MINRES(b, Aprod, Msolve, bilinearfor Line 1082  def MINRES(b, Aprod, Msolve, bilinearfor
1082          w1    = w2          w1    = w2
1083          w2    = w          w2    = w
1084          w     = (v - oldeps*w1 - delta*w2) * denom          w     = (v - oldeps*w1 - delta*w2) * denom
1085          x     = x  +  phi*w          x     +=  phi*w
1086    
1087          # Go round again.          # Go round again.
1088    
# Line 979  def MINRES(b, Aprod, Msolve, bilinearfor Line 1101  def MINRES(b, Aprod, Msolve, bilinearfor
1101          rnorm  = phibar          rnorm  = phibar
1102    
1103      return x      return x
       
 def NewtonGMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium,x=None, iter_max=100,iter_restart=20):  
   
     gamma=.9  
     lmaxit=40  
     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  
     atol=1.e-2  
     rtol=1.e-4  
     stop_tol=atol + rtol*fnrm  
     #  
     # main iteration loop  
     #  
     while not stoppingcriterium(fnrm,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, Aprod, Msolve, bilinearform, stoppingcriterium, xc=x, iter_max=lmaxit-initer, iter_restart=iter_restart, TOL=etamax)  
                if stopped: break  
                initer+=iter_restart  
         xold=x  
         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  
1104    
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    
# Line 1201  class ArithmeticTuple(object): Line 1257  class ArithmeticTuple(object):
1257         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1258         """         """
1259         out=[]         out=[]
1260         for i in range(len(self)):         try:  
1261             out.append(self[i]*other)             l=len(other)
1262               if l!=len(self):
1263                   raise ValueError,"length of of arguments don't match."
1264               for i in range(l): out.append(self[i]*other[i])
1265           except TypeError:
1266           for i in range(len(self)): out.append(self[i]*other)
1267         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1268    
1269     def __rmul__(self,other):     def __rmul__(self,other):
# Line 1215  class ArithmeticTuple(object): Line 1276  class ArithmeticTuple(object):
1276         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1277         """         """
1278         out=[]         out=[]
1279         for i in range(len(self)):         try:  
1280             out.append(other*self[i])             l=len(other)
1281               if l!=len(self):
1282                   raise ValueError,"length of of arguments don't match."
1283               for i in range(l): out.append(other[i]*self[i])
1284           except TypeError:
1285           for i in range(len(self)): out.append(other*self[i])
1286         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1287    
 #########################  
 # Added by Artak  
 #########################  
1288     def __div__(self,other):     def __div__(self,other):
1289         """         """
1290         dividing from the right         dividing from the right
# Line 1231  class ArithmeticTuple(object): Line 1294  class ArithmeticTuple(object):
1294         @return: itemwise self/other         @return: itemwise self/other
1295         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1296         """         """
1297         out=[]         return self*(1/other)
        for i in range(len(self)):  
            out.append(self[i]/other)  
        return ArithmeticTuple(*tuple(out))  
1298    
1299     def __rdiv__(self,other):     def __rdiv__(self,other):
1300         """         """
# Line 1246  class ArithmeticTuple(object): Line 1306  class ArithmeticTuple(object):
1306         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1307         """         """
1308         out=[]         out=[]
1309         for i in range(len(self)):         try:  
1310             out.append(other/self[i])             l=len(other)
1311               if l!=len(self):
1312                   raise ValueError,"length of of arguments don't match."
1313               for i in range(l): out.append(other[i]/self[i])
1314           except TypeError:
1315           for i in range(len(self)): out.append(other/self[i])
1316         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1317        
 ##########################################33  
   
1318     def __iadd__(self,other):     def __iadd__(self,other):
1319         """         """
1320         in-place add of other to self         in-place add of other to self
# Line 1265  class ArithmeticTuple(object): Line 1328  class ArithmeticTuple(object):
1328             self.__items[i]+=other[i]             self.__items[i]+=other[i]
1329         return self         return self
1330    
1331       def __add__(self,other):
1332           """
1333           add other to self
1334    
1335           @param other: increment
1336           @type other: C{ArithmeticTuple}
1337           """
1338           out=[]
1339           try:  
1340               l=len(other)
1341               if l!=len(self):
1342                   raise ValueError,"length of of arguments don't match."
1343               for i in range(l): out.append(self[i]+other[i])
1344           except TypeError:
1345           for i in range(len(self)): out.append(self[i]+other)
1346           return ArithmeticTuple(*tuple(out))
1347    
1348       def __sub__(self,other):
1349           """
1350           subtract other from self
1351    
1352           @param other: increment
1353           @type other: C{ArithmeticTuple}
1354           """
1355           out=[]
1356           try:  
1357               l=len(other)
1358               if l!=len(self):
1359                   raise ValueError,"length of of arguments don't match."
1360               for i in range(l): out.append(self[i]-other[i])
1361           except TypeError:
1362           for i in range(len(self)): out.append(self[i]-other)
1363           return ArithmeticTuple(*tuple(out))
1364      
1365       def __isub__(self,other):
1366           """
1367           subtract other from self
1368    
1369           @param other: increment
1370           @type other: C{ArithmeticTuple}
1371           """
1372           if len(self) != len(other):
1373               raise ValueError,"tuple length must match."
1374           for i in range(len(self)):
1375               self.__items[i]-=other[i]
1376           return self
1377    
1378       def __neg__(self):
1379           """
1380           negate
1381    
1382           """
1383           out=[]
1384           for i in range(len(self)):
1385               out.append(-self[i])
1386           return ArithmeticTuple(*tuple(out))
1387    
1388    
1389  class HomogeneousSaddlePointProblem(object):  class HomogeneousSaddlePointProblem(object):
1390        """        """
1391        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
1392    
1393               Av+B^*p=f               Av+B^*p=f
1394               Bv    =0               Bv    =0
# Line 1336  class HomogeneousSaddlePointProblem(obje Line 1457  class HomogeneousSaddlePointProblem(obje
1457    
1458        def __inner_p(self,p1,p2):        def __inner_p(self,p1,p2):
1459           return self.inner(p1,p2)           return self.inner(p1,p2)
1460          
1461          def __inner_a(self,a1,a2):
1462             return self.inner_a(a1,a2)
1463    
1464          def __inner_a1(self,a1,a2):
1465             return self.inner(a1[1],a2[1])
1466    
1467        def __stoppingcriterium(self,norm_r,r,p):        def __stoppingcriterium(self,norm_r,r,p):
1468            return self.stoppingcriterium(r[1],r[0],p)            return self.stoppingcriterium(r[1],r[0],p)
# Line 1365  class HomogeneousSaddlePointProblem(obje Line 1492  class HomogeneousSaddlePointProblem(obje
1492                # which leads to BA^-1B^*p = BA^-1f                  # which leads to BA^-1B^*p = BA^-1f  
1493    
1494            # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)                  # Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask)      
   
             
1495            self.__z=v+self.solve_A(v,p*0)            self.__z=v+self.solve_A(v,p*0)
1496                  Bz=self.B(self.__z)
               Bz=self.B(self.__z)  
1497                #                #
1498            #   solve BA^-1B^*p = Bz            #   solve BA^-1B^*p = Bz
1499                #                #
               #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv  
1500                #                #
               #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)  
               #                           A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask)  
1501                #                #
1502                self.iter=0                self.iter=0
1503            if solver=='GMRES':                  if solver=='GMRES':      
# Line 1387  class HomogeneousSaddlePointProblem(obje Line 1508  class HomogeneousSaddlePointProblem(obje
1508                  #       u=v+(u-v)                  #       u=v+(u-v)
1509          u=v+self.solve_A(v,p)          u=v+self.solve_A(v,p)
1510    
           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.)  
                 # solve Au=f-B^*p  
                 #       A(u-v)=f-B^*p-Av  
                 #       u=v+(u-v)  
         u=v+self.solve_A(v,p)  
                   
   
1511            if solver=='TFQMR':                  if solver=='TFQMR':      
1512                  if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter                  if self.verbose: print "enter TFQMR method (iter_max=%s)"%max_iter
1513                  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 1412  class HomogeneousSaddlePointProblem(obje Line 1524  class HomogeneousSaddlePointProblem(obje
1524                  #       u=v+(u-v)                  #       u=v+(u-v)
1525          u=v+self.solve_A(v,p)          u=v+self.solve_A(v,p)
1526                                
1527              if solver=='GMRESC':      
1528                    if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter
1529                    p0=self.solve_prec1(Bz)
1530                #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))
1531                    #p-=alfa
1532                    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)
1533                    #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())
1534    
1535                    # solve Au=f-B^*p
1536                    #       A(u-v)=f-B^*p-Av
1537                    #       u=v+(u-v)
1538                p=x[1]
1539            u=v+self.solve_A(v,p)      
1540            #p=x[1]
1541            #u=x[0]
1542    
1543                if solver=='PCG':                if solver=='PCG':
1544                    #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1545                    #
1546                    #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1547                    #                           A(v-z)= f -Az - B^*p (v-z=0 on fixed_u_mask)
1548                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1549                  p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)                  p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p)
1550              u=r[0]                u=r[0]  
1551                    # print "DIFF=",util.integrate(p)
1552    
1553                print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)                # print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1554    
1555            return u,p            return u,p
1556    
# Line 1425  class HomogeneousSaddlePointProblem(obje Line 1558  class HomogeneousSaddlePointProblem(obje
1558            return self.solve_prec(r[1])            return self.solve_prec(r[1])
1559    
1560        def __Msolve2(self,r):        def __Msolve2(self,r):
1561            return self.solve_prec(r)            return self.solve_prec(r*1)
1562    
1563          def __Mempty(self,r):
1564              return r
1565    
1566    
1567        def __Aprod(self,p):        def __Aprod(self,p):
1568            # return BA^-1B*p            # return BA^-1B*p
1569            #solve Av =-B^*p as Av =f-Az-B^*p            #solve Av =B^*p as Av =f-Az-B^*(-p)
1570            v=self.solve_A(self.__z,-p)            v=self.solve_A(self.__z,-p)
1571            return ArithmeticTuple(v, self.B(v))            return ArithmeticTuple(v, self.B(v))
1572    
1573          def __Anext(self,x):
1574              # return next v,p
1575              #solve Av =-B^*p as Av =f-Az-B^*p
1576    
1577          pc=x[1]
1578              v=self.solve_A(self.__z,-pc)
1579          p=self.solve_prec1(self.B(v))
1580    
1581              return ArithmeticTuple(v,p)
1582    
1583    
1584        def __Aprod2(self,p):        def __Aprod2(self,p):
1585            # return BA^-1B*p            # return BA^-1B*p
1586            #solve Av =-B^*p as Av =f-Az-B^*p            #solve Av =B^*p as Av =f-Az-B^*(-p)
1587        v=self.solve_A(self.__z,-p)        v=self.solve_A(self.__z,-p)
1588            return self.B(v)            return self.B(v)
1589    
1590        def __Aprod_Newton(self,p):        def __Aprod_Newton(self,p):
1591            # return BA^-1B*p            # return BA^-1B*p - Bz
1592            #solve Av =-B^*p as Av =f-Az-B^*p            #solve Av =-B^*p as Av =f-Az-B^*p
1593        v=self.solve_A(self.__z,-p)        v=self.solve_A(self.__z,-p)
1594            return self.B(v-self.__z)            return self.B(v-self.__z)
1595    
1596          def __Aprod_Newton2(self,x):
1597              # return BA^-1B*p - Bz
1598              #solve Av =-B^*p as Av =f-Az-B^*p
1599              pc=x[1]
1600          v=self.solve_A(self.__z,-pc)
1601              p=self.solve_prec1(self.B(v-self.__z))
1602              return ArithmeticTuple(v,p)
1603    
1604    
1605    def MaskFromBoundaryTag(domain,*tags):
1606       """
1607       creates a mask on the Solution(domain) function space which one for samples
1608       that touch the boundary tagged by tags.
1609    
1610       usage: m=MaskFromBoundaryTag(domain,"left", "right")
1611    
1612       @param domain: a given domain
1613       @type domain: L{escript.Domain}
1614       @param tags: boundray tags
1615       @type tags: C{str}
1616       @return: a mask which marks samples that are touching the boundary tagged by any of the given tags.
1617       @rtype: L{escript.Data} of rank 0
1618       """
1619       pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1620       d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
1621       for t in tags: d.setTaggedValue(t,1.)
1622       pde.setValue(y=d)
1623       return util.whereNonZero(pde.getRightHandSide())
1624    #============================================================================================================================================
1625  class SaddlePointProblem(object):  class SaddlePointProblem(object):
1626     """     """
1627     This implements a solver for a saddlepoint problem     This implements a solver for a saddlepoint problem
# Line 1471  class SaddlePointProblem(object): Line 1647  class SaddlePointProblem(object):
1647         @type verbose: C{bool}         @type verbose: C{bool}
1648         @note: this method may be overwritten by a particular saddle point problem         @note: this method may be overwritten by a particular saddle point problem
1649         """         """
1650           print "SaddlePointProblem should not be used anymore!"
1651         if not isinstance(verbose,bool):         if not isinstance(verbose,bool):
1652              raise TypeError("verbose needs to be of type bool.")              raise TypeError("verbose needs to be of type bool.")
1653         self.__verbose=verbose         self.__verbose=verbose
1654         self.relaxation=1.         self.relaxation=1.
1655           DeprecationWarning("SaddlePointProblem should not be used anymore.")
1656    
1657     def trace(self,text):     def trace(self,text):
1658         """         """
# Line 1633  class SaddlePointProblem(object): Line 1811  class SaddlePointProblem(object):
1811              norm_u_new = util.Lsup(u_new)              norm_u_new = util.Lsup(u_new)
1812              p_new=p+self.relaxation*g_new              p_new=p+self.relaxation*g_new
1813              norm_p_new = util.sqrt(self.inner(p_new,p_new))              norm_p_new = util.sqrt(self.inner(p_new,p_new))
1814                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))
1815    
             self.trace("step %d: f = %e, f/u = %e, g = %e, g/p = %e, relaxation = %e."%(self.iter, norm_f_new ,norm_f_new/norm_u_new, norm_g_new, norm_g_new/norm_p_new, self.relaxation))  
1816              if self.iter>1:              if self.iter>1:
1817                 dg2=g_new-g                 dg2=g_new-g
1818                 df2=f_new-f                 df2=f_new-f
# Line 1648  class SaddlePointProblem(object): Line 1826  class SaddlePointProblem(object):
1826              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
1827          self.trace("convergence after %s steps."%self.iter)          self.trace("convergence after %s steps."%self.iter)
1828          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(function_space,*tags):  
    """  
    create a mask on the given function space which one for samples  
    that touch the boundary tagged by tags.  
   
    usage: m=MaskFromBoundaryTag(Solution(domain),"left", "right")  
   
    @param function_space: a given function space  
    @type function_space: L{escript.FunctionSpace}  
    @param tags: boundray tags  
    @type tags: C{str}  
    @return: a mask which marks samples used by C{function_space} that are touching the  
             boundary tagged by any of the given tags.  
    @rtype: L{escript.Data} of rank 0  
    """  
    pde=linearPDEs.LinearPDE(function_space.getDomain(),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)  
    out=util.whereNonZero(pde.getRightHandSide())  
    if out.getFunctionSpace() == function_space:  
       return out  
    else:  
       return util.whereNonZero(util.interpolate(out,function_space))  
   
   
   

Legend:
Removed from v.1542  
changed lines
  Added in v.1956

  ViewVC Help
Powered by ViewVC 1.1.26