/[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 2647 by gross, Fri Sep 4 05:25:25 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 349  class Locator: Line 350  class Locator:
350         Returns the coordinates of the Locator as a string.         Returns the coordinates of the Locator as a string.
351         """         """
352         x=self.getX()         x=self.getX()
353         if instance(x,list):         if isinstance(x,list):
354            out="["            out="["
355            first=True            first=True
356            for xx in x:            for xx in x:
# 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 591  class Defect(object): Line 618  class Defect(object):
618      def __call__(self,x):      def __call__(self,x):
619          return self.eval(x)          return self.eval(x)
620    
621      def setDerivativeIncrementLength(self,inc=math.sqrt(util.EPSILON)):      def setDerivativeIncrementLength(self,inc=1000.*math.sqrt(util.EPSILON)):
622          """          """
623          Sets the relative length of the increment used to approximate the          Sets the relative length of the increment used to approximate the
624          derivative of the defect. The increment is inc*norm(x)/norm(v)*v in the          derivative of the defect. The increment is inc*norm(x)/norm(v)*v in the
# 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.2647  
changed lines
  Added in v.2793

  ViewVC Help
Powered by ViewVC 1.1.26