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

revision 1550 by artak, Wed May 7 02:22:57 2008 UTC revision 1557 by artak, Mon May 19 04:44:27 2008 UTC
# 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 1204  class ArithmeticTuple(object): Line 1205  class ArithmeticTuple(object):
1205         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1206         """         """
1207         out=[]         out=[]
1208         for i in range(len(self)):         other=1.*other
1209           if isinstance(other,float):
1210        for i in range(len(self)):
1211             out.append(self[i]*other)             out.append(self[i]*other)
1212           else:
1213            for i in range(len(self)):
1214               out.append(self[i]*other[i])
1215         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1216
1217     def __rmul__(self,other):     def __rmul__(self,other):
# Line 1218  class ArithmeticTuple(object): Line 1224  class ArithmeticTuple(object):
1224         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1225         """         """
1226         out=[]         out=[]
1227         for i in range(len(self)):         other=1.*other
1228           if isinstance(other,float):
1229        for i in range(len(self)):
1230             out.append(other*self[i])             out.append(other*self[i])
1231           else:
1232            for i in range(len(self)):
1233               out.append(other[i]*self[i])
1234         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1235
1236  #########################  #########################
# Line 1235  class ArithmeticTuple(object): Line 1246  class ArithmeticTuple(object):
1246         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1247         """         """
1248         out=[]         out=[]
1249         for i in range(len(self)):         other=1.*other
1250             out.append(self[i]/other)         if isinstance(other,float):
1251        for i in range(len(self)):
1252               out.append(self[i]*(1./other))
1253           else:
1254            for i in range(len(self)):
1255               out.append(self[i]*(1./other[i]))
1256         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1257
1258     def __rdiv__(self,other):     def __rdiv__(self,other):
# Line 1249  class ArithmeticTuple(object): Line 1265  class ArithmeticTuple(object):
1265         @rtype: L{ArithmeticTuple}         @rtype: L{ArithmeticTuple}
1266         """         """
1267         out=[]         out=[]
1268         for i in range(len(self)):         other=1.*other
1269             out.append(other/self[i])         if isinstance(other,float):
1270            for i in range(len(self)):
1271               out.append(other*(1./self[i]))
1272           else:
1273            for i in range(len(self)):
1274               out.append(other[i]*(1./self[i]))
1275         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1276
1277  ##########################################33  ##########################################33
# Line 1275  class ArithmeticTuple(object): Line 1296  class ArithmeticTuple(object):
1296         @param other: increment         @param other: increment
1297         @type other: C{ArithmeticTuple}         @type other: C{ArithmeticTuple}
1298         """         """
1299         if not isinstance(other,float):  #      if not isinstance(other,float):
1300          if len(self) != len(other):         if len(self) != len(other):
1301             raise ValueError,"tuple length must match."            raise ValueError,"tuple length must match."
1302          for i in range(len(self)):         for i in range(len(self)):
1303             self.__items[i]+=other[i]            self.__items[i]+=other[i]
1304         else:  #       else:
1305          for i in range(len(self)):  #        for i in range(len(self)):
1306             self.__items[i]+=other  #           self.__items[i]+=other
1307
1308         return self         return self
1309
# Line 1298  class ArithmeticTuple(object): Line 1319  class ArithmeticTuple(object):
1319         for i in range(len(self)):         for i in range(len(self)):
1320             self.__items[i]-=other[i]             self.__items[i]-=other[i]
1321         return self         return self
1322
1323       def __isub__(self,other):
1324           """
1325           subtract other from self
1326
1327           @param other: increment
1328           @type other: C{ArithmeticTuple}
1329           """
1330           if len(self) != len(other):
1331               raise ValueError,"tuple length must match."
1332           for i in range(len(self)):
1333               self.__items[i]-=other[i]
1334           return self
1335
1336     def __neg__(self):     def __neg__(self):
1337         """         """
1449                # which leads to BA^-1B^*p = BA^-1f                  # which leads to BA^-1B^*p = BA^-1f
1450
1451            # 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)

1452            self.__z=v+self.solve_A(v,p*0)            self.__z=v+self.solve_A(v,p*0)
1453                  Bz=self.B(self.__z)
Bz=self.B(self.__z)
1454                #                #
1455            #   solve BA^-1B^*p = Bz            #   solve BA^-1B^*p = Bz
1456                #                #
#   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1457                #                #
#   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)
1458                #                #
1459                self.iter=0                self.iter=0
1460            if solver=='GMRES':                  if solver=='GMRES':
1492
1493            if solver=='GMRESC':                  if solver=='GMRESC':
1494                  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
1495                  p0=self.solve_prec(Bz)                  p0=self.solve_prec1(Bz)
1496             #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)))
1497                 #p-=alfa                  #p-=alfa
1498                  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)
1499                  #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())
1500
1501                  # solve Au=f-B^*p                  # solve Au=f-B^*p
1502                  #       A(u-v)=f-B^*p-Av                  #       A(u-v)=f-B^*p-Av
1503                  #       u=v+(u-v)                  #       u=v+(u-v)
1504              p=x[1]              p=x[1]
1505          #u=v+self.solve_A(v,p)                u=v+self.solve_A(v,p)
1506          #p=x[1]          #p=x[1]
1507          #u=x[0]          #u=x[0]
1508
1509                if solver=='PCG':                if solver=='PCG':
1510                    #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1511                    #
1512                    #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1513                    #                           A(v-z)= f -Az - B^*p (v-z=0 on fixed_u_mask)
1514                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1515                  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)
1516              u=r[0]                u=r[0]
1517                    print "DIFF=",util.integrate(p)
1518
1519                print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)                print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1520
1532
1533        def __Aprod(self,p):        def __Aprod(self,p):
1534            # return BA^-1B*p            # return BA^-1B*p
1535            #solve Av =-B^*p as Av =f-Az-B^*p            #solve Av =B^*p as Av =f-Az-B^*(-p)
1536            v=self.solve_A(self.__z,-p)            v=self.solve_A(self.__z,-p)
1537            return ArithmeticTuple(v, self.B(v))            return ArithmeticTuple(v, self.B(v))
1538
1542
1543        pc=x[1]        pc=x[1]
1544            v=self.solve_A(self.__z,-pc)            v=self.solve_A(self.__z,-pc)
1545        p=self.solve_prec(self.B(v))        p=self.solve_prec1(self.B(v))
1546
1547            return ArithmeticTuple(v,p)            return ArithmeticTuple(v,p)
1548
1549
1550        def __Aprod2(self,p):        def __Aprod2(self,p):
1551            # return BA^-1B*p            # return BA^-1B*p
1552            #solve Av =-B^*p as Av =f-Az-B^*p            #solve Av =B^*p as Av =f-Az-B^*(-p)
1553        v=self.solve_A(self.__z,-p)        v=self.solve_A(self.__z,-p)
1554            return self.B(v)            return self.B(v)
1555