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

revision 2718 by artak, Tue Sep 29 05:47:45 2009 UTC revision 2719 by gross, Wed Oct 14 06:38:03 2009 UTC
# Line 641  class Defect(object): Line 641  class Defect(object):
641          return (F1-F0)/epsnew          return (F1-F0)/epsnew
642
643  ######################################  ######################################
644  def NewtonGMRES(defect, x, iter_max=100, sub_iter_max=20, atol=0,rtol=1.e-4, sub_tol_max=0.5, gamma=0.9, verbose=False):  def NewtonGMRES(defect, x, iter_max=100, sub_iter_max=20, atol=0,rtol=1.e-4, subtol_max=0.5, gamma=0.9, verbose=False):
645     """     """
646     Solves a non-linear problem *F(x)=0* for unknown *x* using the stopping     Solves a non-linear problem *F(x)=0* for unknown *x* using the stopping
647     criterion:     criterion:
# Line 666  def NewtonGMRES(defect, x, iter_max=100, Line 666  def NewtonGMRES(defect, x, iter_max=100,
666     :type rtol: positive ``float``     :type rtol: positive ``float``
667     :param gamma: tolerance safety factor for inner iteration     :param gamma: tolerance safety factor for inner iteration
668     :type gamma: positive ``float``, less than 1     :type gamma: positive ``float``, less than 1
669     :param sub_tol_max: upper bound for inner tolerance     :param subtol_max: upper bound for inner tolerance
670     :type sub_tol_max: positive ``float``, less than 1     :type subtol_max: positive ``float``, less than 1
671     :return: an approximation of the solution with the desired accuracy     :return: an approximation of the solution with the desired accuracy
672     :rtype: same type as the initial guess     :rtype: same type as the initial guess
673     """     """
# Line 676  def NewtonGMRES(defect, x, iter_max=100, Line 676  def NewtonGMRES(defect, x, iter_max=100,
676     if rtol<0: raise ValueError,"rtol needs to be non-negative."     if rtol<0: raise ValueError,"rtol needs to be non-negative."
677     if rtol+atol<=0: raise ValueError,"rtol or atol needs to be non-negative."     if rtol+atol<=0: raise ValueError,"rtol or atol needs to be non-negative."
678     if gamma<=0 or gamma>=1: raise ValueError,"tolerance safety factor for inner iteration (gamma =%s) needs to be positive and less than 1."%gamma     if gamma<=0 or gamma>=1: raise ValueError,"tolerance safety factor for inner iteration (gamma =%s) needs to be positive and less than 1."%gamma
679     if sub_tol_max<=0 or sub_tol_max>=1: raise ValueError,"upper bound for inner tolerance for inner iteration (sub_tol_max =%s) needs to be positive and less than 1."%sub_tol_max     if subtol_max<=0 or subtol_max>=1: raise ValueError,"upper bound for inner tolerance for inner iteration (subtol_max =%s) needs to be positive and less than 1."%subtol_max
680
681     F=defect(x)     F=defect(x)
682     fnrm=defect.norm(F)     fnrm=defect.norm(F)
683     stop_tol=atol + rtol*fnrm     stop_tol=atol + rtol*fnrm
684     sub_tol=sub_tol_max     subtol=subtol_max
685     if verbose: print "NewtonGMRES: initial residual = %e."%fnrm     if verbose: print "NewtonGMRES: initial residual = %e."%fnrm
686     if verbose: print "             tolerance = %e."%sub_tol     if verbose: print "             tolerance = %e."%subtol
687     iter=1     iter=1
688     #     #
689     # main iteration loop     # main iteration loop
# Line 691  def NewtonGMRES(defect, x, iter_max=100, Line 691  def NewtonGMRES(defect, x, iter_max=100,
691     while not fnrm<=stop_tol:     while not fnrm<=stop_tol:
692              if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max              if iter  >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max
693              #              #
695          #          #
696              if iter > 1:              if iter > 1:
697             rat=fnrm/fnrmo             rat=fnrm/fnrmo
698                 sub_tol_old=sub_tol                 subtol_old=subtol
699             sub_tol=gamma*rat**2             subtol=gamma*rat**2
700             if gamma*sub_tol_old**2 > .1: sub_tol=max(sub_tol,gamma*sub_tol_old**2)             if gamma*subtol_old**2 > .1: subtol=max(subtol,gamma*subtol_old**2)
701             sub_tol=max(min(sub_tol,sub_tol_max), .5*stop_tol/fnrm)             subtol=max(min(subtol,subtol_max), .5*stop_tol/fnrm)
702          #          #
703          # calculate newton increment xc          # calculate newton increment xc
704              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown
# Line 706  def NewtonGMRES(defect, x, iter_max=100, Line 706  def NewtonGMRES(defect, x, iter_max=100,
706              #     if  atol is reached sub_iter returns the numer of steps performed to get there              #     if  atol is reached sub_iter returns the numer of steps performed to get there
707              #              #
708              #              #
709              if verbose: print "             subiteration (GMRES) is called with relative tolerance %e."%sub_tol              if verbose: print "             subiteration (GMRES) is called with relative tolerance %e."%subtol
710              try:              try:
711                 xc, sub_iter=__FDGMRES(F, defect, x, sub_tol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)                 xc, sub_iter=__FDGMRES(F, defect, x, subtol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)
712              except MaxIterReached:              except MaxIterReached:
713                 raise MaxIterReached,"maximum number of %s steps reached."%iter_max                 raise MaxIterReached,"maximum number of %s steps reached."%iter_max
714              if sub_iter<0:              if sub_iter<0:
1472
1473        for the unknowns *v* and *p* and given operators *A* and *B* and        for the unknowns *v* and *p* and given operators *A* and *B* and
1474        given right hand side *f*. *B^** is the adjoint operator of *B*.        given right hand side *f*. *B^** is the adjoint operator of *B*.
1475          *A* may depend weakly on *v* and *p*.
1476        """        """
1477        def __init__(self, adaptSubTolerance=True, **kwargs):        def __init__(self, **kwargs):
1478      """      """

:param adaptSubTolerance: If True the tolerance for subproblem is set automatically.
1480      """      """
1481            self.resetControlParameters()
1482          self.setTolerance()          self.setTolerance()
1483          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
1484      self.__adaptSubTolerance=adaptSubTolerance        def resetControlParameters(self,gamma=0.85, gamma_min=1.e-8,chi_max=0.1, omega_div=0.2, omega_conv=1.1, rtol_min=1.e-7, rtol_max=0.1, chi=1., C_p=1., C_v=1., safety_factor=0.3):
1485             """
1486             sets a control parameter
1487
1488             :param gamma: ``1/(1-gamma)`` controls the perturbation of the converegence rate due to termination errors in the subproblems.
1489             :type gamma: ``float``
1490             :param gamma_min: minimum value for ``gamma``.
1491             :type gamma_min: ``float``
1492             :param chi_max: maximum tolerable converegence rate.
1493             :type chi_max: ``float``
1494             :param omega_div: reduction fact for ``gamma`` if no convergence is detected.
1495             :type omega_div: ``float``
1496             :param omega_conv: raise fact for ``gamma`` if convergence is detected.
1497             :type omega_conv: ``float``
1498             :param rtol_min: minimum relative tolerance used to calculate presssure and velocity increment.
1499             :type rtol_min: ``float``
1500             :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1501             :type rtol_max: ``float``
1502             :param chi: initial convergence measure.
1503             :type chi: ``float``
1504             :param C_p: initial value for constant to adjust pressure tolerance
1505             :type C_p: ``float``
1506             :param C_v: initial value for constant to adjust velocity tolerance
1507             :type C_v: ``float``
1508             :param safety_factor: safety factor for addjustment of pressure and velocity tolerance from stopping criteria
1509             :type safety_factor: ``float``
1510             """
1511             self.setControlParameter(gamma, gamma_min ,chi_max , omega_div , omega_conv, rtol_min , rtol_max, chi,C_p, C_v,safety_factor)
1512
1513          def setControlParameter(self,gamma=None, gamma_min=None ,chi_max=None, omega_div=None, omega_conv=None, rtol_min=None, rtol_max=None, chi=None, C_p=None, C_v=None, safety_factor=None):
1514             """
1515             sets a control parameter
1516
1517             :param gamma: ``1/(1-gamma)`` controls the perturbation of the converegence rate due to termination errors in the subproblems.
1518             :type gamma: ``float``
1519             :param gamma_min: minimum value for ``gamma``.
1520             :type gamma_min: ``float``
1521             :param chi_max: maximum tolerable converegence rate.
1522             :type chi_max: ``float``
1523             :param omega_div: reduction fact for ``gamma`` if no convergence is detected.
1524             :type omega_div: ``float``
1525             :param omega_conv: raise fact for ``gamma`` if convergence is detected.
1526             :type omega_conv: ``float``
1527             :param rtol_min: minimum relative tolerance used to calculate presssure and velocity increment.
1528             :type rtol_min: ``float``
1529             :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1530             :type rtol_max: ``float``
1531             :param chi: initial convergence measure.
1532             :type chi: ``float``
1533             :param C_p: initial value for constant to adjust pressure tolerance
1534             :type C_p: ``float``
1535             :param C_v: initial value for constant to adjust velocity tolerance
1536             :type C_v: ``float``
1537             :param safety_factor: safety factor for addjustment of pressure and velocity tolerance from stopping criteria
1538             :type safety_factor: ``float``
1539             """
1540             if not gamma == None:
1541                if gamma<=0 or gamma>=1:
1542                   raise ValueError,"Initial gamma needs to be positive and less than 1."
1543             else:
1544                gamma=self.__gamma
1545
1546             if not gamma_min == None:
1547                if gamma_min<=0 or gamma_min>=1:
1548                   raise ValueError,"gamma_min needs to be positive and less than 1."
1549             else:
1550                gamma_min = self.__gamma_min
1551
1552             if not chi_max == None:
1553                if chi_max<=0 or chi_max>=1:
1554                   raise ValueError,"chi_max needs to be positive and less than 1."
1555             else:
1556                chi_max = self.__chi_max
1557
1558             if not omega_div == None:
1559                if omega_div<=0 or omega_div >=1:
1560                   raise ValueError,"omega_div needs to be positive and less than 1."
1561             else:
1562                omega_div=self.__omega_div
1563
1564             if not omega_conv == None:
1565                if omega_conv<1:
1566                   raise ValueError,"omega_conv needs to be greater or equal to 1."
1567             else:
1568                omega_conv=self.__omega_conv
1569
1570             if not rtol_min == None:
1571                if rtol_min<=0 or rtol_min>=1:
1572                   raise ValueError,"rtol_min needs to be positive and less than 1."
1573             else:
1574                rtol_min=self.__rtol_min
1575
1576             if not rtol_max == None:
1577                if rtol_max<=0 or rtol_max>=1:
1578                   raise ValueError,"rtol_max needs to be positive and less than 1."
1579             else:
1580                rtol_max=self.__rtol_max
1581
1582             if not chi == None:
1583                if chi<=0 or chi>1:
1584                   raise ValueError,"chi needs to be positive and less than 1."
1585             else:
1586                chi=self.__chi
1587
1588             if not C_p == None:
1589                if C_p<1:
1590                   raise ValueError,"C_p need to be greater or equal to 1."
1591             else:
1592                C_p=self.__C_p
1593
1594             if not C_v == None:
1595                if C_v<1:
1596                   raise ValueError,"C_v need to be greater or equal to 1."
1597             else:
1598                C_v=self.__C_v
1599
1600             if not safety_factor == None:
1601                if safety_factor<=0 or safety_factor>1:
1602                   raise ValueError,"safety_factor need to be between zero and one."
1603             else:
1604                safety_factor=self.__safety_factor
1605
1606             if gamma<gamma_min:
1607                   raise ValueError,"gamma = %e needs to be greater or equal gamma_min = %e."%(gamma,gamma_min)
1608             if rtol_max<=rtol_min:
1609                   raise ValueError,"rtol_max = %e needs to be greater rtol_min = %e."%(rtol_max,rtol_min)
1610
1611             self.__gamma = gamma
1612             self.__gamma_min = gamma_min
1613             self.__chi_max = chi_max
1614             self.__omega_div = omega_div
1615             self.__omega_conv = omega_conv
1616             self.__rtol_min = rtol_min
1617             self.__rtol_max = rtol_max
1618             self.__chi = chi
1619             self.__C_p = C_p
1620             self.__C_v = C_v
1621             self.__safety_factor = safety_factor
1622
1623        #=============================================================        #=============================================================
1624        def initialize(self):        def initialize(self):
1625          """          """
1659           :rtype: non-negative ``float``           :rtype: non-negative ``float``
1660           """           """
1661           raise NotImplementedError,"no norm of v implemented."           raise NotImplementedError,"no norm of v implemented."
1662        def getV(self, p, v0):        def getDV(self, p, v, tol):
1663           """           """
1664           return the value for v for a given p (overwrite)           return a correction to the value for a given v and a given p with accuracy `tol` (overwrite)
1665
1666           :param p: a pressure           :param p: pressure
1667           :param v0: a initial guess for the value v to return.           :param v: pressure
1668           :return: v given as *v= A^{-1} (f-B^*p)*           :return: dv given as *dv= A^{-1} (f-A v-B^*p)*
1669             :note: Only *A* may depend on *v* and *p*
1670           """           """
1671           raise NotImplementedError,"no v calculation implemented."           raise NotImplementedError,"no dv calculation implemented."
1672
1673
1674        def Bv(self,v):        def Bv(self,v, tol):
1675          """          """
1676          Returns Bv (overwrite).          Returns Bv with accuracy `tol` (overwrite)
1677
1678          :rtype: equal to the type of p          :rtype: equal to the type of p
1679          :note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1689          """          """
1690          raise NotImplementedError, "no norm of Bv implemented."          raise NotImplementedError, "no norm of Bv implemented."
1691
1692        def solve_AinvBt(self,p):        def solve_AinvBt(self,dp, tol):
1693           """           """
1694           Solves *Av=B^*p* with accuracy `self.getSubProblemTolerance()`           Solves *A dv=B^*dp* with accuracy `tol`
(overwrite).
1695
1696           :param p: a pressure increment           :param dp: a pressure increment
1697           :return: the solution of *Av=B^*p*           :return: the solution of *A dv=B^*dp*
1698           :note: boundary conditions on v should be zero!           :note: boundary conditions on dv should be zero! *A* is the operator used in ``getDV`` and must not be altered.
1699           """           """
1700           raise NotImplementedError,"no operator A implemented."           raise NotImplementedError,"no operator A implemented."
1701
1702        def solve_prec(self,Bv):        def solve_prec(self,Bv, tol):
1703           """           """
1704           Provides a preconditioner for *(BA^{-1}B^ * )* applied to Bv with accuracy           Provides a preconditioner for *(BA^{-1}B^ * )* applied to Bv with accuracy `tol`
`self.getSubProblemTolerance()` (overwrite).
1705
1706           :rtype: equal to the type of p           :rtype: equal to the type of p
1707           :note: boundary conditions on p should be zero!           :note: boundary conditions on p should be zero!
1708           """           """
1709           raise NotImplementedError,"no preconditioner for Schur complement implemented."           raise NotImplementedError,"no preconditioner for Schur complement implemented."
def setSubProblemTolerance(self):
"""
:note: method is typically the method is overwritten.
"""
pass
1710        #=============================================================        #=============================================================
1711        def __Aprod_PCG(self,p):        def __Aprod_PCG(self,dp):
1712            dv=self.solve_AinvBt(p)            dv=self.solve_AinvBt(dp, self.__subtol)
1713            return ArithmeticTuple(dv,self.Bv(dv))            return ArithmeticTuple(dv,self.Bv(dv, self.__subtol))
1714
1715        def __inner_PCG(self,p,r):        def __inner_PCG(self,p,r):
1716           return self.inner_pBv(p,r[1])           return self.inner_pBv(p,r[1])
1717
1718        def __Msolve_PCG(self,r):        def __Msolve_PCG(self,r):
1719            return self.solve_prec(r[1])            return self.solve_prec(r[1], self.__subtol)
1720        #=============================================================        #=============================================================
1721        def __Aprod_GMRES(self,p):        def __Aprod_GMRES(self,p):
1722            return self.solve_prec(self.Bv(self.solve_AinvBt(p)))            return self.solve_prec(self.Bv(self.solve_AinvBt(p, self.__subtol), self.__subtol), self.__subtol)
1723        def __inner_GMRES(self,p0,p1):        def __inner_GMRES(self,p0,p1):
1724           return self.inner_p(p0,p1)           return self.inner_p(p0,p1)
1725
1735            f=self.inner_p(p,p)            f=self.inner_p(p,p)
1736            if f<0: raise ValueError,"negative pressure norm."            if f<0: raise ValueError,"negative pressure norm."
1737            return math.sqrt(f)            return math.sqrt(f)
"""
Returns True if tolerance adaption for subproblem is choosen.
"""
1738
1739        def solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):        def solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1740           """           """
1760           self.verbose=verbose           self.verbose=verbose
1761           rtol=self.getTolerance()           rtol=self.getTolerance()
1762           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
1763           correction_step=0           correction_step=0
1764           converged=False           converged=False
1765             error=None
1766             chi=None
1767             gamma=self.__gamma
1768             C_p=self.__C_p
1769             C_v=self.__C_v
1770           while not converged:           while not converged:
1771                  if error== None or chi == None:
1772                      gamma_new=gamma/self.__omega_conv
1773                  else:
1774                     if chi < self.__chi_max:
1775                        gamma_new=min(max(gamma*self.__omega_conv,1-chi*error/(self.__safety_factor*ATOL)), 1-chi/self.__chi_max)
1776                     else:
1777                        gamma_new=gamma*self.__omega_div
1778                  gamma=max(gamma_new, self.__gamma_min)
1779                # calculate velocity for current pressure:                # calculate velocity for current pressure:
1780                v=self.getV(p,v)                rtol_v=min(max(gamma/(1.+gamma)/C_v,self.__rtol_min), self.__rtol_max)
1781                Bv=self.Bv(v)                rtol_p=min(max(gamma/C_p, self.__rtol_min), self.__rtol_max)
1782                norm_v=self.norm_v(v)                self.__subtol=rtol_p**2
1783                norm_Bv=self.norm_Bv(Bv)                if self.verbose: print "HomogeneousSaddlePointProblem: step %s: gamma = %e, rtol_v= %e, rtol_p=%e"%(correction_step,gamma,rtol_v,rtol_p)
1784                ATOL=norm_v*rtol+atol                if self.verbose: print "HomogeneousSaddlePointProblem: subtolerance: %e"%self.__subtol
1785                if self.verbose: print "HomogeneousSaddlePointProblem: norm v= %e, norm_Bv= %e, tolerance = %e."%(norm_v, norm_Bv,ATOL)                # calculate velocity for current pressure: A*dv= F-A*v-B*p
1786                  dv1=self.getDV(p,v,rtol_v)
1787                  v1=v+dv1
1788                  Bv1=self.Bv(v1, self.__subtol)
1789                  norm_Bv1=self.norm_Bv(Bv1)
1790                  norm_dv1=self.norm_v(dv1)
1791                  norm_v1=self.norm_v(v1)
1792                  ATOL=norm_v1*rtol+atol
1793                  if self.verbose: print "HomogeneousSaddlePointProblem: step %s: Bv = %e, dv = %e, v=%e"%(correction_step,norm_Bv1, norm_dv1, norm_v1)
1794                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."
1795                if norm_Bv <= ATOL:                if max(norm_Bv1, norm_dv1) <= ATOL:
1796                   converged=True                    converged=True
1797                      v=v1
1798                else:                else:
1799                   correction_step+=1                    # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1
1800                   if correction_step>max_correction_steps:                    if usePCG:
1801                        raise CorrectionFailed,"Given up after %d correction steps."%correction_step                      dp,r,a_norm=PCG(ArithmeticTuple(v1,Bv1),self.__Aprod_PCG,0*p,self.__Msolve_PCG,self.__inner_PCG,atol=0, rtol=rtol_p,iter_max=max_iter, verbose=self.verbose)
1802                   dp=self.solve_prec(Bv)                      v2=r[0]
1803                   if usePCG:                      Bv2=r[1]
1804                     norm2=self.inner_pBv(dp,Bv)                    else:
1805                     if norm2<0: raise ValueError,"negative PCG norm."                      dp=GMRES(self.solve_prec(Bv1,self.__subtol),self.__Aprod_GMRES, 0*p, self.__inner_GMRES,atol=0, rtol=rtol_p,iter_max=max_iter, iter_restart=iter_restart, verbose=self.verbose)
1806                     norm2=math.sqrt(norm2)                      dv2=self.solve_AinvBt(dp, self.__subtol)
1807                   else:                      v2=v1-dv2
1808                     norm2=self.norm_p(dp)                      Bv2=self.Bv(v2, self.__subtol)
1809                   ATOL_ITER=ATOL/norm_Bv*norm2*0.5                    #
1810                   if self.verbose: print "HomogeneousSaddlePointProblem: tolerance for solver: %e"%ATOL_ITER                    # convergence indicators:
1811                   if usePCG:                    #
1812                         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)                    norm_v2=self.norm_v(v2)
1813                   else:                    norm_dv2=self.norm_v(v2-v)
1814                         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)                    error_new=max(norm_dv2, norm_Bv1)
1815           if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached."                    norm_Bv2=self.norm_Bv(Bv2)
1816                      if self.verbose: print "HomogeneousSaddlePointProblem: step %s (part 2): Bv = %e, dv = %e, v=%e"%(correction_step,norm_Bv2, norm_dv2, norm_v2)
1817                      if error !=None:
1818                          chi_new=error_new/error
1819                          if self.verbose: print "HomogeneousSaddlePointProblem: step %s: convergence rate = %e"%(correction_step,chi_new)
1820                          if chi != None:
1821                              gamma0=max(gamma, 1-chi/chi_new)
1822                              C_p*=gamma0/gamma
1823                              C_v*=gamma0/gamma*(1+gamma)/(1+gamma0)
1824                          chi=chi_new
1825                      error = error_new
1826                      correction_step+=1
1827                      if correction_step>max_correction_steps:
1828                            raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1829                      v,p=v2,p+dp
1830             if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step
1831       return v,p       return v,p
1832
1833        #========================================================================        #========================================================================
1871           """           """
1872           return self.__atol           return self.__atol
1873
def getSubProblemTolerance(self):
"""
Sets the relative tolerance to solve the subproblem(s).
"""
return max(200.*util.EPSILON,self.getTolerance()**2)
1874