/[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 2250 by gross, Wed Feb 4 06:27:59 2009 UTC revision 2251 by gross, Fri Feb 6 06:50:39 2009 UTC
# Line 886  def GMRES(r, Aprod, x, bilinearform, ato Line 886  def GMRES(r, Aprod, x, bilinearform, ato
886        if stopped: break        if stopped: break
887        if verbose: print "GMRES: restart."        if verbose: print "GMRES: restart."
888        restarted=True        restarted=True
889     if verbose: print "GMRES: tolerance has reached."     if verbose: print "GMRES: tolerance has been reached."
890     return x     return x
891    
892  def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False):  def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False):
# Line 1465  class HomogeneousSaddlePointProblem(obje Line 1465  class HomogeneousSaddlePointProblem(obje
1465          """          """
1466          pass          pass
1467    
1468        def B(self,v):        def inner_pBv(self,p,v):
         """  
         Returns Bv (overwrite).  
   
         @rtype: equal to the type of p  
         @note: boundary conditions on p should be zero!  
         """  
         raise NotImplementedError, "no operator B implemented."  
   
       def inner_pBv(self,p,Bv):  
1469           """           """
1470           Returns inner product of element p and Bv (overwrite).           Returns inner product of element p and Bv (overwrite).
1471    
1472           @type p: equal to the type of p           @param p: a pressure increment
1473           @type Bv: equal to the type of result of operator B           @param v: a residual
1474           @rtype: equal to the type of p           @return: inner product of element p and Bv
1475             @rtype: C{float}
1476             @note: used if PCG is applied.
1477           """           """
1478           raise NotImplementedError,"no inner product for p implemented."           raise NotImplementedError,"no inner product for p implemented."
1479    
1480        def inner_p(self,p0,p1):        def inner_p(self,p0,p1):
1481           """           """
1482           Returns inner product of element p0 and p1 (overwrite).           Returns inner product of p0 and p1 (overwrite).
1483    
1484           @type p0: equal to the type of p           @param p0: a pressure
1485           @type p1: equal to the type of p           @param p1: a pressure
1486           @rtype: equal to the type of p           @return: inner product of p0 and p1
1487             @rtype: C{float}
1488           """           """
1489           raise NotImplementedError,"no inner product for p implemented."           raise NotImplementedError,"no inner product for p implemented."
1490      
1491          def norm_v(self,v):
1492             """
1493             Returns the norm of v (overwrite).
1494    
1495        def inner_v(self,v0,v1):           @param v: a velovity
1496             @return: norm of v
1497             @rtype: non-negative C{float}
1498           """           """
1499           Returns inner product of two elements v0 and v1 (overwrite).           raise NotImplementedError,"no norm of v implemented."
1500    
1501    
1502           @type v0: equal to the type of v        def getV(self, p, v0):
          @type v1: equal to the type of v  
          @rtype: equal to the type of v  
1503           """           """
1504           raise NotImplementedError,"no inner product for v implemented."           return the value for v for a given p (overwrite)
          pass  
1505    
1506        def solve_A(self,u,p):           @param p: a pressure
1507             @param v0: a initial guess for the value v to return.
1508             @return: v given as M{v= A^{-1} (f-B^*p)}
1509           """           """
1510           Solves M{Av=f-Au-B^*p} with accuracy L{self.getSubProblemTolerance()}           raise NotImplementedError,"no v calculation implemented."
1511            
1512          def norm_Bv(self,v):
1513            """
1514            Returns Bv (overwrite).
1515    
1516            @rtype: equal to the type of p
1517            @note: boundary conditions on p should be zero!
1518            """
1519            raise NotImplementedError, "no operator B implemented."
1520    
1521          def solve_AinvBt(self,p):
1522             """
1523             Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
1524           (overwrite).           (overwrite).
1525    
1526           @rtype: equal to the type of v           @param p: a pressure increment
1527             @return: the solution of M{Av=B^*p}
1528           @note: boundary conditions on v should be zero!           @note: boundary conditions on v should be zero!
1529           """           """
1530           raise NotImplementedError,"no operator A implemented."           raise NotImplementedError,"no operator A implemented."
1531    
1532        def solve_prec(self,p):        def solve_precB(self,v):
1533           """           """
1534           Provides a preconditioner for M{BA^{-1}B^*} with accuracy           Provides a preconditioner for M{BA^{-1}B^*} with accuracy
1535           L{self.getSubProblemTolerance()} (overwrite).           L{self.getSubProblemTolerance()} (overwrite).
# Line 1526  class HomogeneousSaddlePointProblem(obje Line 1540  class HomogeneousSaddlePointProblem(obje
1540           raise NotImplementedError,"no preconditioner for Schur complement implemented."           raise NotImplementedError,"no preconditioner for Schur complement implemented."
1541        #=============================================================        #=============================================================
1542        def __Aprod_PCG(self,p):        def __Aprod_PCG(self,p):
1543            # return (v,Bv) with v=A^-1B*p            return self.solve_AinvBt(p)
           #solve Av =B^*p as Av =f-Az-B^*(-p)  
           v=self.solve_A(self.__z,-p)  
           return ArithmeticTuple(v, self.B(v))  
1544    
1545        def __inner_PCG(self,p,r):        def __inner_PCG(self,p,v):
1546           return self.inner_pBv(p,r[1])           return self.inner_pBv(p,v)
1547    
1548        def __Msolve_PCG(self,r):        def __Msolve_PCG(self,v):
1549            return self.solve_prec(r[1])            return self.solve_precB(v)
1550        #=============================================================        #=============================================================
1551        def __Aprod_GMRES(self,x):        def __Aprod_GMRES(self,p):
1552            # return w,q  from v, p            return self.solve_precB(self.solve_AinvBt(p))
1553            # solve Aw =Av+B^*p as Aw =f-A(z-v)-B^*(-p)        def __inner_GMRES(self,p0,p1):
1554            #  and  Sq=B(v-w)           return self.inner_p(p0,p1)
           v=x[0]  
           p=x[1]  
           w=self.solve_A(self.__z-v,-p)  
           Bw=self.B(w-v)  
           q=self.solve_prec(Bw)  
           return ArithmeticTuple(w,q)  
   
       def __inner_GMRES(self,a1,a2):  
          return self.inner_p(a1[1],a2[1])+self.inner_v(a1[0],a2[0])  
   
1555        #=============================================================        #=============================================================
1556        def norm(self,v,p):        def norm_p(self,p):
1557          f=self.inner_p(p,p)+self.inner_v(v,v)            """
1558          if f<0:            calculates the norm of C{p}
1559              raise ValueError,"negative norm."            
1560          return math.sqrt(f)            @param p: a pressure
1561              @return: the norm of C{p} using the inner product for pressure
1562              @rtype: C{float}
1563              """
1564              f=self.inner_p(p,p)
1565              if f<0: raise ValueError,"negative pressure norm."
1566              return math.sqrt(f)
1567              
1568    
1569        def solve(self,v,p,max_iter=20, verbose=False, show_details=False, useUzawa=True, iter_restart=20,max_correction_steps=3):        def solve(self,v,p,max_iter=20, verbose=False, show_details=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1570           """           """
1571           Solves the saddle point problem using initial guesses v and p.           Solves the saddle point problem using initial guesses v and p.
1572    
# Line 1566  class HomogeneousSaddlePointProblem(obje Line 1574  class HomogeneousSaddlePointProblem(obje
1574           @param p: initial guess for pressure           @param p: initial guess for pressure
1575           @type v: L{Data}           @type v: L{Data}
1576           @type p: L{Data}           @type p: L{Data}
1577           @param useUzawa: indicates the usage of the Uzawa scheme. Otherwise           @param usePCG: indicates the usage of the PCG rather than GMRES scheme.
                           the problem is solved in coupled form.  
1578           @param max_iter: maximum number of iteration steps per correction           @param max_iter: maximum number of iteration steps per correction
1579                            attempt                            attempt
1580           @param verbose: if True, shows information on the progress of the           @param verbose: if True, shows information on the progress of the
# Line 1575  class HomogeneousSaddlePointProblem(obje Line 1582  class HomogeneousSaddlePointProblem(obje
1582           @param show_details: if True, shows details of the sub problem solver           @param show_details: if True, shows details of the sub problem solver
1583           @param iter_restart: restart the iteration after C{iter_restart} steps           @param iter_restart: restart the iteration after C{iter_restart} steps
1584                                (only used if useUzaw=False)                                (only used if useUzaw=False)
1585           @param max_correction_steps: maximum number of iteration steps in the           @type usePCG: C{bool}
                                       attempt to get |Bv| to zero.  
          @return: new approximations for velocity and pressure  
          @type useUzawa: C{bool}  
1586           @type max_iter: C{int}           @type max_iter: C{int}
1587           @type verbose: C{bool}           @type verbose: C{bool}
1588           @type show_details: C{bool}           @type show_details: C{bool}
1589           @type iter_restart: C{int}           @type iter_restart: C{int}
          @type max_correction_steps: C{int}  
1590           @rtype: C{tuple} of L{Data} objects           @rtype: C{tuple} of L{Data} objects
1591           """           """
1592           self.verbose=verbose           self.verbose=verbose
1593           self.show_details=show_details and self.verbose           self.show_details=show_details and self.verbose
          #=====================================================================  
          # Az=f is solved as A(z-v)=f-Av-B^*0 (typically with z-v=0 on boundary)  
          self.__z=v+self.solve_A(v,p*0)  
          # tolerances:  
1594           rtol=self.getTolerance()           rtol=self.getTolerance()
1595           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
          if useUzawa:  
             p0=self.solve_prec(self.B(self.__z))  
             f0=self.norm(self.__z,p0)  
          else:  
             f0=util.sqrt(self.inner_v(self.__z,self.__z))  
          if not f0>0: return self.__z, p*0  
          ATOL=rtol*f0+atol  
          if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."  
          if self.verbose: print "saddle point solver: initial residual %e, tolerance = %e."%(f0,ATOL)  
          # initialization  
          self.iter=0  
1596           correction_step=0           correction_step=0
1597           converged=False           converged=False
1598           # initial guess:           while not converged:
1599           q=p*1                # calculate velocity for current pressure:
1600           u=v*1                v=self.getV(p,v)
1601           while not converged :                #
1602              if useUzawa:                norm_v=self.norm_v(v)
1603                 # assume p is known: then v=z-A^-1B^*p                norm_Bv=self.norm_Bv(v)
1604                 #                ATOL=norm_v*rtol+atol
1605                 # which leads to BA^-1B^*p = Bz                if self.verbose: print "saddle point solver: norm v= %e, norm_Bv= %e, tolerance = %e."%(norm_v, norm_Bv,ATOL)
1606                 #                if not ATOL>0: raise ValueError,"overall absolute tolerance needs to be positive."
1607                 # note that the residual r=Bz-BA^-1B^*p = B(z-A^-1B^*p) = Bv                if norm_Bv <= ATOL:
1608                 # we use the (v,Bv) to represent the residual                   converged=True
1609                 #                else:
1610                 # the norm of the right hand side Bv = f0                   correction_step+=1
1611                 #                   if correction_step>max_correction_steps:
1612                 #                  and the initial residual                        raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1613                 #                   dp=self.solve_precB(v)
1614                 #     r=Bz-BA^-1B^*q = B(z-A^{-1}B^*q)=Bw with A(w-z)=Az-Az-B^*q = f -A*0 - B^{*}q                   if usePCG:
1615                 #                     norm2=self.inner_pBv(dp,v)
1616                 w=self.solve_A(self.__z,q)+self.__z                     if norm2<0: raise ValueError,"negative PCG norm."
1617                 if self.verbose: print "enter PCG method (iter_max=%s, atol=%e, subtolerance=%e)"%(max_iter,ATOL, self.getSubProblemTolerance())                     norm2=math.sqrt(norm2)
1618                 q,r=PCG(ArithmeticTuple(w,self.B(w)),self.__Aprod_PCG,q,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)                   else:
1619                 u=r[0]                     norm2=self.norm_p(dp)
1620                 Bu=r[1]                   ATOL_ITER=ATOL/norm_Bv*norm2
1621              else:                   if self.verbose: print "saddle point solver: tolerance for solver: %e"%ATOL_ITER
1622                 #                   if usePCG:
1623                 #  with v=dv+z                         p,v0=PCG(v,self.__Aprod_PCG,p,self.__Msolve_PCG,self.__inner_PCG,atol=ATOL_ITER, rtol=0.,iter_max=max_iter, verbose=self.verbose)
1624                 #                   else:
1625                 #   A v + B p = f                         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)
                #   B v       = 0  
                #  
                # apply the preconditioner [[A^{-1} 0][(S^{-1} B A^{-1})  -S^{-1}]]  
                #  
                w=self.solve_A(u,q)  
                if self.verbose: print "enter GMRES (iter_max=%s, atol=%e, subtolerance=%e)"%(max_iter,ATOL,self.getSubProblemTolerance())  
                x=GMRES(ArithmeticTuple(w,self.solve_prec(self.B(u+w))),self.__Aprod_GMRES, ArithmeticTuple(u,q), \  
                          self.__inner_GMRES,atol=ATOL, rtol=0.,iter_max=max_iter, iter_restart=iter_restart, verbose=self.verbose)  
                u=x[0]  
                q=x[1]  
                Bu=self.B(u)  
             # now we check if |Bu| ~ 0 or more precise |Bu|_p  <= rtol * |v|_v  
             nu=self.inner_v(u,u)  
             p2=self.solve_prec(Bu)  
             nBu=self.inner_p(p2,p2)  
             if not nu>=0 and not nBu>=0: raise NegativeNorm,"negative norm."  
             nu=math.sqrt(nu)  
             nBu=math.sqrt(nBu)  
             if self.verbose: print "saddle point solver: norm v= %e (Bv = %e)"%(nu,nBu)  
             QTOL=atol+nu*rtol  
             if nBu <= QTOL:  
                 converged=True  
             else:  
                 ATOL=QTOL/nBu*ATOL*0.3  
                 if self.verbose: print "correction step %s: tolerance reduced to %e."%(correction_step,ATOL)  
                 converged=False  
             correction_step+=1  
             if correction_step>max_correction_steps:  
                raise CorrectionFailed,"Given up after %d correction steps."%correction_step  
1626           if self.verbose: print "saddle point solver: tolerance reached."           if self.verbose: print "saddle point solver: tolerance reached."
1627       return u,q       return v,p
1628    
1629        #========================================================================        #========================================================================
1630        def setTolerance(self,tolerance=1.e-4):        def setTolerance(self,tolerance=1.e-4):

Legend:
Removed from v.2250  
changed lines
  Added in v.2251

  ViewVC Help
Powered by ViewVC 1.1.26