/[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 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):
# Line 1415  class HomogeneousSaddlePointProblem(obje Line 1415  class HomogeneousSaddlePointProblem(obje
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':      
# Line 1478  class HomogeneousSaddlePointProblem(obje Line 1472  class HomogeneousSaddlePointProblem(obje
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    
# Line 1498  class HomogeneousSaddlePointProblem(obje Line 1497  class HomogeneousSaddlePointProblem(obje
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    
# Line 1515  class HomogeneousSaddlePointProblem(obje Line 1514  class HomogeneousSaddlePointProblem(obje
1514    
1515        def __Aprod2(self,p):        def __Aprod2(self,p):
1516            # return BA^-1B*p            # return BA^-1B*p
1517            #solve Av =-B^*p as Av =f-Az-B^*p            #solve Av =B^*p as Av =f-Az-B^*(-p)
1518        v=self.solve_A(self.__z,-p)        v=self.solve_A(self.__z,-p)
1519            return self.B(v)            return self.B(v)
1520    

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

  ViewVC Help
Powered by ViewVC 1.1.26