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

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):
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).
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
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
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