/[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 1552 by gross, Thu May 8 08:52:41 2008 UTC revision 1659 by gross, Fri Jul 18 02:28:13 2008 UTC
# Line 474  class IterationHistory(object): Line 474  class IterationHistory(object):
474    
475         """         """
476         self.history.append(norm_r)         self.history.append(norm_r)
477         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])
478         return self.history[-1]<=self.tolerance * self.history[0]         return self.history[-1]<=self.tolerance * self.history[0]
479    
480     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 493  class IterationHistory(object):
493         if TOL==None:         if TOL==None:
494            TOL=self.tolerance            TOL=self.tolerance
495         self.history.append(norm_r)         self.history.append(norm_r)
496         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])
497         return self.history[-1]<=TOL * norm_b         return self.history[-1]<=TOL * norm_b
498    
499  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 618  def GMRESm(b, Aprod, Msolve, bilinearfor Line 618  def GMRESm(b, Aprod, Msolve, bilinearfor
618        
619     v.append(r/rho)     v.append(r/rho)
620     g[0]=rho     g[0]=rho
621    
622     while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):     while not (stoppingcriterium(rho,norm_b) or iter==iter_restart-1):
623    
624      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
625    
       
626      p=Msolve(Aprod(v[iter]))      p=Msolve(Aprod(v[iter]))
627    
628      v.append(p)      v.append(p)
# Line 630  def GMRESm(b, Aprod, Msolve, bilinearfor Line 630  def GMRESm(b, Aprod, Msolve, bilinearfor
630      v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))        v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
631    
632  # Modified Gram-Schmidt  # Modified Gram-Schmidt
633      for j in range(iter+1):      for j in range(iter+1):
634        h[j][iter]=bilinearform(v[j],v[iter+1])          h[j][iter]=bilinearform(v[j],v[iter+1])  
635        v[iter+1]-=h[j][iter]*v[j]        v[iter+1]-=h[j][iter]*v[j]
636                
637      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]))
638      v_norm2=h[iter+1][iter]      v_norm2=h[iter+1][iter]
639    
   
640  # Reorthogonalize if needed  # Reorthogonalize if needed
641      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)
642       for j in range(iter+1):       for j in range(iter+1):  
643          hr=bilinearform(v[j],v[iter+1])          hr=bilinearform(v[j],v[iter+1])
644              h[j][iter]=h[j][iter]+hr #vhat              h[j][iter]=h[j][iter]+hr
645              v[iter+1] -= hr*v[j]              v[iter+1] -= hr*v[j]
646    
647       v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))         v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))  
# Line 660  def GMRESm(b, Aprod, Msolve, bilinearfor Line 659  def GMRESm(b, Aprod, Msolve, bilinearfor
659              for i in range(iter+1) : h[i][iter]=hhat[i]              for i in range(iter+1) : h[i][iter]=hhat[i]
660    
661      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])
662    
663      if mu!=0 :      if mu!=0 :
664          c[iter]=h[iter][iter]/mu          c[iter]=h[iter][iter]/mu
665          s[iter]=-h[iter+1][iter]/mu          s[iter]=-h[iter+1][iter]/mu
# Line 668  def GMRESm(b, Aprod, Msolve, bilinearfor Line 668  def GMRESm(b, Aprod, Msolve, bilinearfor
668          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])
669    
670  # Update the residual norm  # Update the residual norm
671                  
672          rho=abs(g[iter+1])          rho=abs(g[iter+1])
673      iter+=1      iter+=1
674    
# Line 686  def GMRESm(b, Aprod, Msolve, bilinearfor Line 687  def GMRESm(b, Aprod, Msolve, bilinearfor
687       for i in range(iter-1):       for i in range(iter-1):
688      xhat += v[i]*y[i]      xhat += v[i]*y[i]
689     else : xhat=v[0]     else : xhat=v[0]
690        
691     x += xhat     x += xhat
692     if iter<iter_restart-1:     if iter<iter_restart-1:
693        stopped=True        stopped=True
# Line 1043  def NewtonGMRES(b, Aprod, Msolve, biline Line 1044  def NewtonGMRES(b, Aprod, Msolve, biline
1044              etanew=gamma*rat*rat              etanew=gamma*rat*rat
1045              if gamma*etaold*etaold > .1 :              if gamma*etaold*etaold > .1 :
1046                  etanew=max(etanew,gamma*etaold*etaold)                  etanew=max(etanew,gamma*etaold*etaold)
               
1047              etamax=min(etanew,etamax)              etamax=min(etanew,etamax)
1048              etamax=max(etamax,.5*stop_tol/fnrm)              etamax=max(etamax,.5*stop_tol/fnrm)
   
1049      return x      return x
1050    
1051  def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):  def TFQMR(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100):
# Line 1204  class ArithmeticTuple(object): Line 1203  class ArithmeticTuple(object):
1203         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1204         """         """
1205         out=[]         out=[]
1206         for i in range(len(self)):         other=1.*other
1207           if isinstance(other,float):
1208        for i in range(len(self)):
1209             out.append(self[i]*other)             out.append(self[i]*other)
1210           else:
1211            for i in range(len(self)):
1212               out.append(self[i]*other[i])
1213         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1214    
1215     def __rmul__(self,other):     def __rmul__(self,other):
# Line 1218  class ArithmeticTuple(object): Line 1222  class ArithmeticTuple(object):
1222         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1223         """         """
1224         out=[]         out=[]
1225         for i in range(len(self)):         other=1.*other
1226           if isinstance(other,float):
1227        for i in range(len(self)):
1228             out.append(other*self[i])             out.append(other*self[i])
1229           else:
1230            for i in range(len(self)):
1231               out.append(other[i]*self[i])
1232         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1233    
1234  #########################  #########################
# Line 1235  class ArithmeticTuple(object): Line 1244  class ArithmeticTuple(object):
1244         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1245         """         """
1246         out=[]         out=[]
1247         for i in range(len(self)):         other=1.*other
1248           if isinstance(other,float):
1249        for i in range(len(self)):
1250             out.append(self[i]*(1./other))             out.append(self[i]*(1./other))
1251           else:
1252            for i in range(len(self)):
1253               out.append(self[i]*(1./other[i]))
1254         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1255    
1256     def __rdiv__(self,other):     def __rdiv__(self,other):
# Line 1249  class ArithmeticTuple(object): Line 1263  class ArithmeticTuple(object):
1263         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1264         """         """
1265         out=[]         out=[]
1266         for i in range(len(self)):         other=1.*other
1267             out.append(other/self[i])         if isinstance(other,float):
1268            for i in range(len(self)):
1269               out.append(other*(1./self[i]))
1270           else:
1271            for i in range(len(self)):
1272               out.append(other[i]*(1./self[i]))
1273         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1274        
1275  ##########################################33  ##########################################33
# Line 1275  class ArithmeticTuple(object): Line 1294  class ArithmeticTuple(object):
1294         @param other: increment         @param other: increment
1295         @type other: C{ArithmeticTuple}         @type other: C{ArithmeticTuple}
1296         """         """
1297         if not isinstance(other,float):  #      if not isinstance(other,float):
1298          if len(self) != len(other):         if len(self) != len(other):
1299             raise ValueError,"tuple length must match."            raise ValueError,"tuple length must match."
1300          for i in range(len(self)):         for i in range(len(self)):
1301             self.__items[i]+=other[i]            self.__items[i]+=other[i]
1302         else:  #       else:
1303          for i in range(len(self)):  #        for i in range(len(self)):
1304             self.__items[i]+=other  #           self.__items[i]+=other
1305    
1306         return self         return self
1307    
# Line 1298  class ArithmeticTuple(object): Line 1317  class ArithmeticTuple(object):
1317         for i in range(len(self)):         for i in range(len(self)):
1318             self.__items[i]-=other[i]             self.__items[i]-=other[i]
1319         return self         return self
1320      
1321       def __isub__(self,other):
1322           """
1323           subtract other from self
1324    
1325           @param other: increment
1326           @type other: C{ArithmeticTuple}
1327           """
1328           if len(self) != len(other):
1329               raise ValueError,"tuple length must match."
1330           for i in range(len(self)):
1331               self.__items[i]-=other[i]
1332           return self
1333    
1334     def __neg__(self):     def __neg__(self):
1335         """         """
# Line 1416  class HomogeneousSaddlePointProblem(obje Line 1448  class HomogeneousSaddlePointProblem(obje
1448    
1449            # 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)      
1450            self.__z=v+self.solve_A(v,p*0)            self.__z=v+self.solve_A(v,p*0)
1451                Bz=self.B(self.__z)                Bz=self.B(self.__z)
1452                #                #
1453            #   solve BA^-1B^*p = Bz            #   solve BA^-1B^*p = Bz
1454                #                #
# Line 1458  class HomogeneousSaddlePointProblem(obje Line 1490  class HomogeneousSaddlePointProblem(obje
1490                                
1491            if solver=='GMRESC':                  if solver=='GMRESC':      
1492                  if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter                  if self.verbose: print "enter GMRES coupled method (iter_max=%s)"%max_iter
1493                  p0=self.solve_prec(Bz)                  p0=self.solve_prec1(Bz)
1494             #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))              #alfa=(1/self.vol)*util.integrate(util.interpolate(p,escript.Function(self.domain)))
1495                 #p-=alfa                  #p-=alfa
1496                  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=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)
1497                  #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())                  #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())
1498    
1499                  # solve Au=f-B^*p                  # solve Au=f-B^*p
1500                  #       A(u-v)=f-B^*p-Av                  #       A(u-v)=f-B^*p-Av
1501                  #       u=v+(u-v)                  #       u=v+(u-v)
1502              p=x[1]              p=x[1]
1503          #u=v+self.solve_A(v,p)                u=v+self.solve_A(v,p)      
1504          #p=x[1]          #p=x[1]
1505          #u=x[0]          #u=x[0]
1506    
# Line 1479  class HomogeneousSaddlePointProblem(obje Line 1512  class HomogeneousSaddlePointProblem(obje
1512                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1513                  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)
1514              u=r[0]                u=r[0]  
1515                  print "DIFF=",util.integrate(p)                  # print "DIFF=",util.integrate(p)
1516    
1517                print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)                # print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1518    
1519            return u,p            return u,p
1520    
# Line 1507  class HomogeneousSaddlePointProblem(obje Line 1540  class HomogeneousSaddlePointProblem(obje
1540    
1541        pc=x[1]        pc=x[1]
1542            v=self.solve_A(self.__z,-pc)            v=self.solve_A(self.__z,-pc)
1543        p=self.solve_prec(self.B(v))        p=self.solve_prec1(self.B(v))
1544    
1545            return ArithmeticTuple(v,p)            return ArithmeticTuple(v,p)
1546    
# Line 1529  class HomogeneousSaddlePointProblem(obje Line 1562  class HomogeneousSaddlePointProblem(obje
1562            #solve Av =-B^*p as Av =f-Az-B^*p            #solve Av =-B^*p as Av =f-Az-B^*p
1563            pc=x[1]            pc=x[1]
1564        v=self.solve_A(self.__z,-pc)        v=self.solve_A(self.__z,-pc)
1565            p=self.solve_prec(self.B(v-self.__z))            p=self.solve_prec1(self.B(v-self.__z))
1566            return ArithmeticTuple(v,p)            return ArithmeticTuple(v,p)
1567    
1568  class SaddlePointProblem(object):  class SaddlePointProblem(object):
# Line 1557  class SaddlePointProblem(object): Line 1590  class SaddlePointProblem(object):
1590         @type verbose: C{bool}         @type verbose: C{bool}
1591         @note: this method may be overwritten by a particular saddle point problem         @note: this method may be overwritten by a particular saddle point problem
1592         """         """
1593           print "SaddlePointProblem should not be used anymore!"
1594         if not isinstance(verbose,bool):         if not isinstance(verbose,bool):
1595              raise TypeError("verbose needs to be of type bool.")              raise TypeError("verbose needs to be of type bool.")
1596         self.__verbose=verbose         self.__verbose=verbose
# Line 1569  class SaddlePointProblem(object): Line 1603  class SaddlePointProblem(object):
1603         @param text: a text message         @param text: a text message
1604         @type text: C{str}         @type text: C{str}
1605         """         """
1606         if self.__verbose: print "#s: #s"%(str(self),text)         if self.__verbose: print "%s: %s"%(str(self),text)
1607    
1608     def solve_f(self,u,p,tol=1.e-8):     def solve_f(self,u,p,tol=1.e-8):
1609         """         """

Legend:
Removed from v.1552  
changed lines
  Added in v.1659

  ViewVC Help
Powered by ViewVC 1.1.26