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

revision 1550 by artak, Wed May 7 02:22:57 2008 UTC revision 1552 by gross, Thu May 8 08:52:41 2008 UTC
# Line 1236  class ArithmeticTuple(object): Line 1236  class ArithmeticTuple(object):
1236         """         """
1237         out=[]         out=[]
1238         for i in range(len(self)):         for i in range(len(self)):
1239             out.append(self[i]/other)             out.append(self[i]*(1./other))
1240         return ArithmeticTuple(*tuple(out))         return ArithmeticTuple(*tuple(out))
1241
1242     def __rdiv__(self,other):     def __rdiv__(self,other):
1415                # which leads to BA^-1B^*p = BA^-1f                  # which leads to BA^-1B^*p = BA^-1f
1416
1417            # 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)

1418            self.__z=v+self.solve_A(v,p*0)            self.__z=v+self.solve_A(v,p*0)

1419                Bz=self.B(self.__z)                Bz=self.B(self.__z)
1420                #                #
1421            #   solve BA^-1B^*p = Bz            #   solve BA^-1B^*p = Bz
1422                #                #
#   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1423                #                #
#   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)
1424                #                #
1425                self.iter=0                self.iter=0
1426            if solver=='GMRES':                  if solver=='GMRES':
1472          #u=x[0]          #u=x[0]
1473
1474                if solver=='PCG':                if solver=='PCG':
1475                    #   note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv
1476                    #
1477                    #   with                    Av=Az-B^*p = f - B^*p (v=z on fixed_u_mask)
1478                    #                           A(v-z)= f -Az - B^*p (v-z=0 on fixed_u_mask)
1479                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter                  if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter
1480                  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)
1481              u=r[0]                u=r[0]
1482                    print "DIFF=",util.integrate(p)
1483
1484                print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)                print "RESULT div(u)=",util.Lsup(self.B(u)),util.Lsup(u)
1485
1497
1498        def __Aprod(self,p):        def __Aprod(self,p):
1499            # return BA^-1B*p            # return BA^-1B*p
1500            #solve Av =-B^*p as Av =f-Az-B^*p            #solve Av =B^*p as Av =f-Az-B^*(-p)
1501            v=self.solve_A(self.__z,-p)            v=self.solve_A(self.__z,-p)
1502            return ArithmeticTuple(v, self.B(v))            return ArithmeticTuple(v, self.B(v))
1503