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

revision 2444 by gross, Wed May 27 08:45:55 2009 UTC revision 2445 by gross, Fri May 29 03:23:25 2009 UTC
# Line 1468  class HomogeneousSaddlePointProblem(obje Line 1468  class HomogeneousSaddlePointProblem(obje
1468          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
1469          self.setSubProblemTolerance()          self.setSubProblemTolerance()
1470
1471
1472        #=============================================================        #=============================================================
1473        def initialize(self):        def initialize(self):
1474          """          """
# Line 1475  class HomogeneousSaddlePointProblem(obje Line 1476  class HomogeneousSaddlePointProblem(obje
1476          """          """
1477          pass          pass
1478
1479        def inner_pBv(self,p,v):        def inner_pBv(self,p,Bv):
1480           """           """
1481           Returns inner product of element p and Bv (overwrite).           Returns inner product of element p and Bv (overwrite).
1482
# Line 1485  class HomogeneousSaddlePointProblem(obje Line 1486  class HomogeneousSaddlePointProblem(obje
1486           @rtype: C{float}           @rtype: C{float}
1487           @note: used if PCG is applied.           @note: used if PCG is applied.
1488           """           """
1489           raise NotImplementedError,"no inner product for p implemented."           raise NotImplementedError,"no inner product for p and Bv implemented."
1490
1491        def inner_p(self,p0,p1):        def inner_p(self,p0,p1):
1492           """           """
# Line 1507  class HomogeneousSaddlePointProblem(obje Line 1508  class HomogeneousSaddlePointProblem(obje
1508           @rtype: non-negative C{float}           @rtype: non-negative C{float}
1509           """           """
1510           raise NotImplementedError,"no norm of v implemented."           raise NotImplementedError,"no norm of v implemented."

1511        def getV(self, p, v0):        def getV(self, p, v0):
1512           """           """
1513           return the value for v for a given p (overwrite)           return the value for v for a given p (overwrite)
# Line 1518  class HomogeneousSaddlePointProblem(obje Line 1517  class HomogeneousSaddlePointProblem(obje
1517           @return: v given as M{v= A^{-1} (f-B^*p)}           @return: v given as M{v= A^{-1} (f-B^*p)}
1518           """           """
1519           raise NotImplementedError,"no v calculation implemented."           raise NotImplementedError,"no v calculation implemented."
1520
1521
1522        def norm_Bv(self,v):        def Bv(self,v):
1523          """          """
1524          Returns Bv (overwrite).          Returns Bv (overwrite).
1525
# Line 1528  class HomogeneousSaddlePointProblem(obje Line 1528  class HomogeneousSaddlePointProblem(obje
1528          """          """
1529          raise NotImplementedError, "no operator B implemented."          raise NotImplementedError, "no operator B implemented."
1530
1531          def norm_Bv(self,Bv):
1532            """
1533            Returns the norm of Bv (overwrite).
1534
1535            @rtype: equal to the type of p
1536            @note: boundary conditions on p should be zero!
1537            """
1538            raise NotImplementedError, "no norm of Bv implemented."
1539
1540        def solve_AinvBt(self,p):        def solve_AinvBt(self,p):
1541           """           """
1542           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
# Line 1539  class HomogeneousSaddlePointProblem(obje Line 1548  class HomogeneousSaddlePointProblem(obje
1548           """           """
1549           raise NotImplementedError,"no operator A implemented."           raise NotImplementedError,"no operator A implemented."
1550
1551        def solve_precB(self,v):        def solve_prec(self,Bv):
1552           """           """
1553           Provides a preconditioner for M{BA^{-1}B^*} with accuracy           Provides a preconditioner for M{BA^{-1}B^*} applied to Bv with accuracy
1554           L{self.getSubProblemTolerance()} (overwrite).           L{self.getSubProblemTolerance()} (overwrite).
1555
1556           @rtype: equal to the type of p           @rtype: equal to the type of p
# Line 1550  class HomogeneousSaddlePointProblem(obje Line 1559  class HomogeneousSaddlePointProblem(obje
1559           raise NotImplementedError,"no preconditioner for Schur complement implemented."           raise NotImplementedError,"no preconditioner for Schur complement implemented."
1560        #=============================================================        #=============================================================
1561        def __Aprod_PCG(self,p):        def __Aprod_PCG(self,p):
1562            return self.solve_AinvBt(p)            dv=self.solve_AinvBt(p)
1563              return ArithmeticTuple(dv,self.Bv(dv))
1564
1565        def __inner_PCG(self,p,v):        def __inner_PCG(self,p,r):
1566           return self.inner_pBv(p,v)           return self.inner_pBv(p,r[1])
1567
1568        def __Msolve_PCG(self,v):        def __Msolve_PCG(self,r):
1569            return self.solve_precB(v)            return self.solve_prec(r[1])
1570        #=============================================================        #=============================================================
1571    # rename solve_prec and change argument v to Bv
1572    # chnage the argument of inner_pBv to v->Bv
1573    # add Bv
1574    # inner p still needed?
1575    # change norm_Bv argument to Bv
1576        def __Aprod_GMRES(self,p):        def __Aprod_GMRES(self,p):
1577            return self.solve_precB(self.solve_AinvBt(p))            return self.solve_prec(self.Bv(self.solve_AinvBt(p)))
1578        def __inner_GMRES(self,p0,p1):        def __inner_GMRES(self,p0,p1):
1579           return self.inner_p(p0,p1)           return self.inner_p(p0,p1)
1580        #=============================================================        #=============================================================
# Line 1608  class HomogeneousSaddlePointProblem(obje Line 1623  class HomogeneousSaddlePointProblem(obje
1623           while not converged:           while not converged:
1624                # calculate velocity for current pressure:                # calculate velocity for current pressure:
1625                v=self.getV(p,v)                v=self.getV(p,v)
1626                  Bv=self.Bv(v)
1627                norm_v=self.norm_v(v)                norm_v=self.norm_v(v)
1628                norm_Bv=self.norm_Bv(v)                norm_Bv=self.norm_Bv(Bv)
1629                ATOL=norm_v*rtol+atol                ATOL=norm_v*rtol+atol
1630                if self.verbose: print "HomogeneousSaddlePointProblem: norm v= %e, norm_Bv= %e, tolerance = %e."%(norm_v, norm_Bv,ATOL)                if self.verbose: print "HomogeneousSaddlePointProblem: norm v= %e, norm_Bv= %e, tolerance = %e."%(norm_v, norm_Bv,ATOL)
1631                if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."                if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."
# Line 1619  class HomogeneousSaddlePointProblem(obje Line 1635  class HomogeneousSaddlePointProblem(obje
1635                   correction_step+=1                   correction_step+=1
1636                   if correction_step>max_correction_steps:                   if correction_step>max_correction_steps:
1637                        raise CorrectionFailed,"Given up after %d correction steps."%correction_step                        raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1638                   dp=self.solve_precB(v)                   dp=self.solve_prec(Bv)
1639                   if usePCG:                   if usePCG:
1640                     norm2=self.inner_pBv(dp,v)                     norm2=self.inner_pBv(dp,Bv)
1641                     if norm2<0: raise ValueError,"negative PCG norm."                     if norm2<0: raise ValueError,"negative PCG norm."
1642                     norm2=math.sqrt(norm2)                     norm2=math.sqrt(norm2)
1643                   else:                   else:
# Line 1629  class HomogeneousSaddlePointProblem(obje Line 1645  class HomogeneousSaddlePointProblem(obje
1645                   ATOL_ITER=ATOL/norm_Bv*norm2*0.5                   ATOL_ITER=ATOL/norm_Bv*norm2*0.5
1646                   if self.verbose: print "HomogeneousSaddlePointProblem: tolerance for solver: %e"%ATOL_ITER                   if self.verbose: print "HomogeneousSaddlePointProblem: tolerance for solver: %e"%ATOL_ITER
1647                   if usePCG:                   if usePCG:
1648                         p,v0,a_norm=PCG(v,self.__Aprod_PCG,p,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL_ITER, rtol=0.,iter_max=max_iter, verbose=self.verbose)                         p,v0,a_norm=PCG(ArithmeticTuple(v,Bv),self.__Aprod_PCG,p,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL_ITER, rtol=0.,iter_max=max_iter, verbose=self.verbose)
1649                   else:                   else:
1650                         p=GMRES(dp,self.__Aprod_GMRES, p, self.__inner_GMRES,atol=ATOL_ITER, rtol=0.,iter_max=max_iter, iter_restart=iter_restart, verbose=self.verbose)                         p=GMRES(dp,self.__Aprod_GMRES, p, self.__inner_GMRES,atol=ATOL_ITER, rtol=0.,iter_max=max_iter, iter_restart=iter_restart, verbose=self.verbose)
1651           if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached."           if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached."

Legend:
 Removed from v.2444 changed lines Added in v.2445