/[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 2683 by gross, Tue Sep 29 02:20:22 2009 UTC revision 2793 by gross, Tue Dec 1 06:10:10 2009 UTC
# Line 146  class Projector: Line 146  class Projector:
146      self.__pde.setValue(D = 1.)      self.__pde.setValue(D = 1.)
147      return      return
148    def getSolverOptions(self):    def getSolverOptions(self):
149      """      """
150      Returns the solver options of the PDE solver.      Returns the solver options of the PDE solver.
151            
152      :rtype: `linearPDEs.SolverOptions`      :rtype: `linearPDEs.SolverOptions`
153      """      """
154        return self.__pde.getSolverOptions()
155    
156    def __call__(self, input_data):    def __call__(self, input_data):
157      """      """
# Line 421  class Locator: Line 422  class Locator:
422          else:          else:
423             return data             return data
424    
425    
426    def getInfLocator(arg):
427        """
428        Return a Locator for a point with the inf value over all arg.
429        """
430        if not isinstance(arg, escript.Data):
431        raise TypeError,"getInfLocator: Unknown argument type."
432        a_inf=util.inf(arg)
433        loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords
434        x=arg.getFunctionSpace().getX()
435        x_min=x.getTupleForGlobalDataPoint(*loc)
436        return Locator(arg.getFunctionSpace(),x_min)
437    
438    def getSupLocator(arg):
439        """
440        Return a Locator for a point with the sup value over all arg.
441        """
442        if not isinstance(arg, escript.Data):
443        raise TypeError,"getInfLocator: Unknown argument type."
444        a_inf=util.sup(arg)
445        loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords
446        x=arg.getFunctionSpace().getX()
447        x_min=x.getTupleForGlobalDataPoint(*loc)
448        return Locator(arg.getFunctionSpace(),x_min)
449        
450    
451  class SolverSchemeException(Exception):  class SolverSchemeException(Exception):
452     """     """
453     This is a generic exception thrown by solvers.     This is a generic exception thrown by solvers.
# Line 640  class Defect(object): Line 667  class Defect(object):
667          return (F1-F0)/epsnew          return (F1-F0)/epsnew
668    
669  ######################################  ######################################
670  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):
671     """     """
672     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
673     criterion:     criterion:
# Line 665  def NewtonGMRES(defect, x, iter_max=100, Line 692  def NewtonGMRES(defect, x, iter_max=100,
692     :type rtol: positive ``float``     :type rtol: positive ``float``
693     :param gamma: tolerance safety factor for inner iteration     :param gamma: tolerance safety factor for inner iteration
694     :type gamma: positive ``float``, less than 1     :type gamma: positive ``float``, less than 1
695     :param sub_tol_max: upper bound for inner tolerance     :param subtol_max: upper bound for inner tolerance
696     :type sub_tol_max: positive ``float``, less than 1     :type subtol_max: positive ``float``, less than 1
697     :return: an approximation of the solution with the desired accuracy     :return: an approximation of the solution with the desired accuracy
698     :rtype: same type as the initial guess     :rtype: same type as the initial guess
699     """     """
# Line 675  def NewtonGMRES(defect, x, iter_max=100, Line 702  def NewtonGMRES(defect, x, iter_max=100,
702     if rtol<0: raise ValueError,"rtol needs to be non-negative."     if rtol<0: raise ValueError,"rtol needs to be non-negative."
703     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."
704     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
705     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
706    
707     F=defect(x)     F=defect(x)
708     fnrm=defect.norm(F)     fnrm=defect.norm(F)
709     stop_tol=atol + rtol*fnrm     stop_tol=atol + rtol*fnrm
710     sub_tol=sub_tol_max     subtol=subtol_max
711     if verbose: print "NewtonGMRES: initial residual = %e."%fnrm     if verbose: print "NewtonGMRES: initial residual = %e."%fnrm
712     if verbose: print "             tolerance = %e."%sub_tol     if verbose: print "             tolerance = %e."%subtol
713     iter=1     iter=1
714     #     #
715     # main iteration loop     # main iteration loop
# Line 690  def NewtonGMRES(defect, x, iter_max=100, Line 717  def NewtonGMRES(defect, x, iter_max=100,
717     while not fnrm<=stop_tol:     while not fnrm<=stop_tol:
718              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
719              #              #
720          #   adjust sub_tol_          #   adjust subtol_
721          #          #
722              if iter > 1:              if iter > 1:
723             rat=fnrm/fnrmo             rat=fnrm/fnrmo
724                 sub_tol_old=sub_tol                 subtol_old=subtol
725             sub_tol=gamma*rat**2             subtol=gamma*rat**2
726             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)
727             sub_tol=max(min(sub_tol,sub_tol_max), .5*stop_tol/fnrm)             subtol=max(min(subtol,subtol_max), .5*stop_tol/fnrm)
728          #          #
729          # calculate newton increment xc          # calculate newton increment xc
730              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown              #     if iter_max in __FDGMRES is reached MaxIterReached is thrown
# Line 705  def NewtonGMRES(defect, x, iter_max=100, Line 732  def NewtonGMRES(defect, x, iter_max=100,
732              #     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
733              #              #
734              #              #
735              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
736              try:              try:
737                 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)
738              except MaxIterReached:              except MaxIterReached:
739                 raise MaxIterReached,"maximum number of %s steps reached."%iter_max                 raise MaxIterReached,"maximum number of %s steps reached."%iter_max
740              if sub_iter<0:              if sub_iter<0:
# Line 1471  class HomogeneousSaddlePointProblem(obje Line 1498  class HomogeneousSaddlePointProblem(obje
1498    
1499        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
1500        given right hand side *f*. *B^** is the adjoint operator of *B*.        given right hand side *f*. *B^** is the adjoint operator of *B*.
1501          *A* may depend weakly on *v* and *p*.
1502        """        """
1503        def __init__(self, adaptSubTolerance=True, **kwargs):        def __init__(self, **kwargs):
1504      """      """
1505      initializes the saddle point problem      initializes the saddle point problem
       
     :param adaptSubTolerance: If True the tolerance for subproblem is set automatically.  
     :type adaptSubTolerance: ``bool``  
1506      """      """
1507            self.resetControlParameters()
1508          self.setTolerance()          self.setTolerance()
1509          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
1510      self.__adaptSubTolerance=adaptSubTolerance        def resetControlParameters(self,gamma=0.85, gamma_min=1.e-2,chi_max=0.1, omega_div=0.2, omega_conv=1.1, rtol_min=1.e-7, rtol_max=0.9, chi=1., C_p=1., C_v=1., safety_factor=0.3):
1511        #=============================================================           """
1512        def initialize(self):           sets a control parameter
         """  
         Initializes the problem (overwrite).  
         """  
         pass  
1513    
1514             :param gamma: ``1/(1-gamma)`` controls the perturbation of the converegence rate due to termination errors in the subproblems.
1515             :type gamma: ``float``
1516             :param gamma_min: minimum value for ``gamma``.
1517             :type gamma_min: ``float``
1518             :param chi_max: maximum tolerable converegence rate.
1519             :type chi_max: ``float``
1520             :param omega_div: reduction fact for ``gamma`` if no convergence is detected.
1521             :type omega_div: ``float``
1522             :param omega_conv: raise fact for ``gamma`` if convergence is detected.
1523             :type omega_conv: ``float``
1524             :param rtol_min: minimum relative tolerance used to calculate presssure and velocity increment.
1525             :type rtol_min: ``float``
1526             :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1527             :type rtol_max: ``float``
1528             :param chi: initial convergence measure.
1529             :type chi: ``float``
1530             :param C_p: initial value for constant to adjust pressure tolerance
1531             :type C_p: ``float``
1532             :param C_v: initial value for constant to adjust velocity tolerance
1533             :type C_v: ``float``
1534             :param safety_factor: safety factor for addjustment of pressure and velocity tolerance from stopping criteria
1535             :type safety_factor: ``float``
1536             """
1537             self.setControlParameter(gamma, gamma_min ,chi_max , omega_div , omega_conv, rtol_min , rtol_max, chi,C_p, C_v,safety_factor)
1538    
1539          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):
1540             """
1541             sets a control parameter
1542    
1543             :param gamma: ``1/(1-gamma)`` controls the perturbation of the converegence rate due to termination errors in the subproblems.
1544             :type gamma: ``float``
1545             :param gamma_min: minimum value for ``gamma``.
1546             :type gamma_min: ``float``
1547             :param chi_max: maximum tolerable converegence rate.
1548             :type chi_max: ``float``
1549             :param omega_div: reduction fact for ``gamma`` if no convergence is detected.
1550             :type omega_div: ``float``
1551             :param omega_conv: raise fact for ``gamma`` if convergence is detected.
1552             :type omega_conv: ``float``
1553             :param rtol_min: minimum relative tolerance used to calculate presssure and velocity increment.
1554             :type rtol_min: ``float``
1555             :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1556             :type rtol_max: ``float``
1557             :param chi: initial convergence measure.
1558             :type chi: ``float``
1559             :param C_p: initial value for constant to adjust pressure tolerance
1560             :type C_p: ``float``
1561             :param C_v: initial value for constant to adjust velocity tolerance
1562             :type C_v: ``float``
1563             :param safety_factor: safety factor for addjustment of pressure and velocity tolerance from stopping criteria
1564             :type safety_factor: ``float``
1565             """
1566             if not gamma == None:
1567                if gamma<=0 or gamma>=1:
1568                   raise ValueError,"Initial gamma needs to be positive and less than 1."
1569             else:
1570                gamma=self.__gamma
1571    
1572             if not gamma_min == None:
1573                if gamma_min<=0 or gamma_min>=1:
1574                   raise ValueError,"gamma_min needs to be positive and less than 1."
1575             else:
1576                gamma_min = self.__gamma_min
1577    
1578             if not chi_max == None:
1579                if chi_max<=0 or chi_max>=1:
1580                   raise ValueError,"chi_max needs to be positive and less than 1."
1581             else:
1582                chi_max = self.__chi_max
1583    
1584             if not omega_div == None:
1585                if omega_div<=0 or omega_div >=1:
1586                   raise ValueError,"omega_div needs to be positive and less than 1."
1587             else:
1588                omega_div=self.__omega_div
1589    
1590             if not omega_conv == None:
1591                if omega_conv<1:
1592                   raise ValueError,"omega_conv needs to be greater or equal to 1."
1593             else:
1594                omega_conv=self.__omega_conv
1595    
1596             if not rtol_min == None:
1597                if rtol_min<=0 or rtol_min>=1:
1598                   raise ValueError,"rtol_min needs to be positive and less than 1."
1599             else:
1600                rtol_min=self.__rtol_min
1601    
1602             if not rtol_max == None:
1603                if rtol_max<=0 or rtol_max>=1:
1604                   raise ValueError,"rtol_max needs to be positive and less than 1."
1605             else:
1606                rtol_max=self.__rtol_max
1607    
1608             if not chi == None:
1609                if chi<=0 or chi>1:
1610                   raise ValueError,"chi needs to be positive and less than 1."
1611             else:
1612                chi=self.__chi
1613    
1614             if not C_p == None:
1615                if C_p<1:
1616                   raise ValueError,"C_p need to be greater or equal to 1."
1617             else:
1618                C_p=self.__C_p
1619    
1620             if not C_v == None:
1621                if C_v<1:
1622                   raise ValueError,"C_v need to be greater or equal to 1."
1623             else:
1624                C_v=self.__C_v
1625    
1626             if not safety_factor == None:
1627                if safety_factor<=0 or safety_factor>1:
1628                   raise ValueError,"safety_factor need to be between zero and one."
1629             else:
1630                safety_factor=self.__safety_factor
1631    
1632             if gamma<gamma_min:
1633                   raise ValueError,"gamma = %e needs to be greater or equal gamma_min = %e."%(gamma,gamma_min)
1634             if rtol_max<=rtol_min:
1635                   raise ValueError,"rtol_max = %e needs to be greater rtol_min = %e."%(rtol_max,rtol_min)
1636                
1637             self.__gamma = gamma
1638             self.__gamma_min = gamma_min
1639             self.__chi_max = chi_max
1640             self.__omega_div = omega_div
1641             self.__omega_conv = omega_conv
1642             self.__rtol_min = rtol_min
1643             self.__rtol_max = rtol_max
1644             self.__chi = chi
1645             self.__C_p = C_p
1646             self.__C_v = C_v
1647             self.__safety_factor = safety_factor
1648    
1649          #=============================================================
1650        def inner_pBv(self,p,Bv):        def inner_pBv(self,p,Bv):
1651           """           """
1652           Returns inner product of element p and Bv (overwrite).           Returns inner product of element p and Bv (overwrite).
# Line 1521  class HomogeneousSaddlePointProblem(obje Line 1679  class HomogeneousSaddlePointProblem(obje
1679           :rtype: non-negative ``float``           :rtype: non-negative ``float``
1680           """           """
1681           raise NotImplementedError,"no norm of v implemented."           raise NotImplementedError,"no norm of v implemented."
1682        def getV(self, p, v0):        def getDV(self, p, v, tol):
1683           """           """
1684           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)
1685    
1686           :param p: a pressure           :param p: pressure
1687           :param v0: a initial guess for the value v to return.           :param v: pressure
1688           :return: v given as *v= A^{-1} (f-B^*p)*           :return: dv given as *dv= A^{-1} (f-A v-B^*p)*
1689             :note: Only *A* may depend on *v* and *p*
1690           """           """
1691           raise NotImplementedError,"no v calculation implemented."           raise NotImplementedError,"no dv calculation implemented."
1692    
1693                    
1694        def Bv(self,v):        def Bv(self,v, tol):
1695          """          """
1696          Returns Bv (overwrite).          Returns Bv with accuracy `tol` (overwrite)
1697    
1698          :rtype: equal to the type of p          :rtype: equal to the type of p
1699          :note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
# Line 1550  class HomogeneousSaddlePointProblem(obje Line 1709  class HomogeneousSaddlePointProblem(obje
1709          """          """
1710          raise NotImplementedError, "no norm of Bv implemented."          raise NotImplementedError, "no norm of Bv implemented."
1711    
1712        def solve_AinvBt(self,p):        def solve_AinvBt(self,dp, tol):
1713           """           """
1714           Solves *Av=B^*p* with accuracy `self.getSubProblemTolerance()`           Solves *A dv=B^*dp* with accuracy `tol`
          (overwrite).  
1715    
1716           :param p: a pressure increment           :param dp: a pressure increment
1717           :return: the solution of *Av=B^*p*           :return: the solution of *A dv=B^*dp*
1718           :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.
1719           """           """
1720           raise NotImplementedError,"no operator A implemented."           raise NotImplementedError,"no operator A implemented."
1721    
1722        def solve_prec(self,Bv):        def solve_prec(self,Bv, tol):
1723           """           """
1724           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).  
1725    
1726           :rtype: equal to the type of p           :rtype: equal to the type of p
1727           :note: boundary conditions on p should be zero!           :note: boundary conditions on p should be zero!
1728           """           """
1729           raise NotImplementedError,"no preconditioner for Schur complement implemented."           raise NotImplementedError,"no preconditioner for Schur complement implemented."
       def setSubProblemTolerance(self):  
          """  
      Updates the tolerance for subproblems  
      :note: method is typically the method is overwritten.  
          """  
          pass  
1730        #=============================================================        #=============================================================
1731        def __Aprod_PCG(self,p):        def __Aprod_PCG(self,dp):
1732            dv=self.solve_AinvBt(p)            dv=self.solve_AinvBt(dp, self.__subtol)
1733            return ArithmeticTuple(dv,self.Bv(dv))            return ArithmeticTuple(dv,self.Bv(dv, self.__subtol))
1734    
1735        def __inner_PCG(self,p,r):        def __inner_PCG(self,p,r):
1736           return self.inner_pBv(p,r[1])           return self.inner_pBv(p,r[1])
1737    
1738        def __Msolve_PCG(self,r):        def __Msolve_PCG(self,r):
1739            return self.solve_prec(r[1])            return self.solve_prec(r[1], self.__subtol)
1740        #=============================================================        #=============================================================
1741        def __Aprod_GMRES(self,p):        def __Aprod_GMRES(self,p):
1742            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)
1743        def __inner_GMRES(self,p0,p1):        def __inner_GMRES(self,p0,p1):
1744           return self.inner_p(p0,p1)           return self.inner_p(p0,p1)
1745    
# Line 1604  class HomogeneousSaddlePointProblem(obje Line 1755  class HomogeneousSaddlePointProblem(obje
1755            f=self.inner_p(p,p)            f=self.inner_p(p,p)
1756            if f<0: raise ValueError,"negative pressure norm."            if f<0: raise ValueError,"negative pressure norm."
1757            return math.sqrt(f)            return math.sqrt(f)
       def adaptSubTolerance(self):  
       """  
       Returns True if tolerance adaption for subproblem is choosen.  
       """  
           return self.__adaptSubTolerance  
1758                
1759        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):
1760           """           """
# Line 1630  class HomogeneousSaddlePointProblem(obje Line 1776  class HomogeneousSaddlePointProblem(obje
1776           :type verbose: ``bool``           :type verbose: ``bool``
1777           :type iter_restart: ``int``           :type iter_restart: ``int``
1778           :rtype: ``tuple`` of `Data` objects           :rtype: ``tuple`` of `Data` objects
1779             :note: typically this method is overwritten by a subclass. It provides a wrapper for the ``_solve`` method.
1780             """
1781             return self._solve(v=v,p=p,max_iter=max_iter,verbose=verbose, usePCG=usePCG, iter_restart=iter_restart, max_correction_steps=max_correction_steps)
1782    
1783          def _solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1784             """
1785             see `_solve` method.
1786           """           """
1787           self.verbose=verbose           self.verbose=verbose
1788           rtol=self.getTolerance()           rtol=self.getTolerance()
1789           atol=self.getAbsoluteTolerance()           atol=self.getAbsoluteTolerance()
      if self.adaptSubTolerance(): self.setSubProblemTolerance()  
1790           correction_step=0           correction_step=0
1791           converged=False           converged=False
1792             error=None
1793             chi=None
1794             gamma=self.__gamma
1795             C_p=self.__C_p
1796             C_v=self.__C_v
1797           while not converged:           while not converged:
1798                  if error== None or chi == None:
1799                      gamma_new=gamma/self.__omega_conv
1800                  else:
1801                     if chi < self.__chi_max:
1802                        gamma_new=min(max(gamma*self.__omega_conv,1-chi*error/(self.__safety_factor*ATOL)), 1-chi/self.__chi_max)
1803                     else:
1804                        gamma_new=gamma*self.__omega_div
1805                  gamma=max(gamma_new, self.__gamma_min)
1806                # calculate velocity for current pressure:                # calculate velocity for current pressure:
1807                v=self.getV(p,v)                rtol_v=min(max(gamma/(1.+gamma)/C_v,self.__rtol_min), self.__rtol_max)
1808                Bv=self.Bv(v)                rtol_p=min(max(gamma/C_p, self.__rtol_min), self.__rtol_max)
1809                norm_v=self.norm_v(v)                self.__subtol=rtol_p**2
1810                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)
1811                ATOL=norm_v*rtol+atol                if self.verbose: print "HomogeneousSaddlePointProblem: subtolerance: %e"%self.__subtol
1812                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
1813                  dv1=self.getDV(p,v,rtol_v)
1814                  v1=v+dv1
1815                  Bv1=self.Bv(v1, self.__subtol)
1816                  norm_Bv1=self.norm_Bv(Bv1)
1817                  norm_dv1=self.norm_v(dv1)
1818                  norm_v1=self.norm_v(v1)
1819                  ATOL=norm_v1*rtol+atol
1820                  if self.verbose: print "HomogeneousSaddlePointProblem: step %s: Bv = %e, dv = %e, v=%e"%(correction_step,norm_Bv1, norm_dv1, norm_v1)
1821                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."
1822                if norm_Bv <= ATOL:                if max(norm_Bv1, norm_dv1) <= ATOL:
1823                   converged=True                    converged=True
1824                      v=v1
1825                else:                else:
1826                   correction_step+=1                    # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1
1827                   if correction_step>max_correction_steps:                    if usePCG:
1828                        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)
1829                   dp=self.solve_prec(Bv)                      v2=r[0]
1830                   if usePCG:                      Bv2=r[1]
1831                     norm2=self.inner_pBv(dp,Bv)                    else:
1832                     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)
1833                     norm2=math.sqrt(norm2)                      dv2=self.solve_AinvBt(dp, self.__subtol)
1834                   else:                      v2=v1-dv2
1835                     norm2=self.norm_p(dp)                      Bv2=self.Bv(v2, self.__subtol)
1836                   ATOL_ITER=ATOL/norm_Bv*norm2*0.5                    #
1837                   if self.verbose: print "HomogeneousSaddlePointProblem: tolerance for solver: %e"%ATOL_ITER                    # convergence indicators:
1838                   if usePCG:                    #
1839                         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)
1840                   else:                    norm_dv2=self.norm_v(v2-v)
1841                         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)
1842           if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached."                    norm_Bv2=self.norm_Bv(Bv2)
1843                      if self.verbose: print "HomogeneousSaddlePointProblem: step %s (part 2): Bv = %e, dv = %e, v=%e"%(correction_step,norm_Bv2, norm_dv2, norm_v2)
1844                      if error !=None:
1845                          chi_new=error_new/error
1846                          if self.verbose: print "HomogeneousSaddlePointProblem: step %s: convergence rate = %e, est. error = %e"%(correction_step,chi_new, error_new)
1847                          if chi != None:
1848                              gamma0=max(gamma, 1-chi/chi_new)
1849                              C_p*=gamma0/gamma
1850                              C_v*=gamma0/gamma*(1+gamma)/(1+gamma0)
1851                          chi=chi_new
1852                      else:
1853                          if self.verbose: print "HomogeneousSaddlePointProblem: step %s: est. error = %e"%(correction_step, error_new)
1854    
1855                      error = error_new
1856                      correction_step+=1
1857                      if correction_step>max_correction_steps:
1858                            raise CorrectionFailed,"Given up after %d correction steps."%correction_step
1859                      v,p=v2,p+dp
1860             if self.verbose: print "HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step
1861       return v,p       return v,p
1862    
1863        #========================================================================        #========================================================================
# Line 1709  class HomogeneousSaddlePointProblem(obje Line 1901  class HomogeneousSaddlePointProblem(obje
1901           """           """
1902           return self.__atol           return self.__atol
1903    
       def getSubProblemTolerance(self):  
          """  
          Sets the relative tolerance to solve the subproblem(s).  
          """  
          return max(200.*util.EPSILON,self.getTolerance()**2)  
1904    
1905  def MaskFromBoundaryTag(domain,*tags):  def MaskFromBoundaryTag(domain,*tags):
1906     """     """

Legend:
Removed from v.2683  
changed lines
  Added in v.2793

  ViewVC Help
Powered by ViewVC 1.1.26