Add InterfaceKinetics.advance_coverages integrator options to the cyt…

…hon interface, and test

include/cantera/kinetics/ImplicitSurfChem.h

@@ -63,14 +63,17 @@ class ImplicitSurfChem : public FuncEval
      *           freedom representing the concentration of surface adsorbates.
      * @param rtol   The relative tolerance for the integrator
      * @param atol   The absolute tolerance for the integrator
-     * @param maxStepSize   The maximum step-size the integrator is allowed to take
-     * @param maxSteps   the maximum number of time-steps the integrator can take
+     * @param maxStepSize   The maximum step-size the integrator is allowed to take.
+     *                      If zero, this option is disabled.
+     * @param maxSteps   The maximum number of time-steps the integrator can take.
+     *                   If not supplied, uses the default value in CVODES (500).
      * @param maxErrTestFails   the maximum permissible number of error test failures
+     *                           If not supplied, uses the default value in CVODES (7).
      */
     ImplicitSurfChem(std::vector<InterfaceKinetics*> k,
-                     doublereal rtol=1.e-7, doublereal atol=1.e-14,
-                     doublereal maxStepSize=0, size_t maxSteps=0,
-                     size_t maxErrTestFails=0);
+                     double rtol=1.e-7, double atol=1.e-14,
+                     double maxStepSize=0, size_t maxSteps=500,
+                     size_t maxErrTestFails=7);
 
     virtual ~ImplicitSurfChem() {};
 
@@ -79,6 +82,27 @@ class ImplicitSurfChem : public FuncEval
      */
     virtual void initialize(doublereal t0 = 0.0);
 
+    /*!
+     *  Set the maximum integration step-size.  Note, setting this value to zero
+     *  disables this option
+     */
+    virtual void setMaxStepSize(double maxstep = 0.0);
+
+    /*!
+     *  Set the relative and absolute integration tolerances.
+     */
+    virtual void setTolerances(double rtol=1.e-7, double atol=1.e-14);
+
+    /*!
+     *  Set the maximum number of CVODES integration steps.
+     */
+    virtual void setMaxSteps(size_t maxsteps = 500);
+
+    /*!
+     *  Set the maximum number of CVODES error test failures
+     */
+    virtual void setMaxErrTestFails(size_t maxErrTestFails = 7);
+
     //! Integrate from t0 to t1. The integrator is reinitialized first.
     /*!
      *   This routine does a time accurate solve from t = t0 to t = t1.

include/cantera/kinetics/InterfaceKinetics.h

@@ -242,14 +242,16 @@ class InterfaceKinetics : public Kinetics
      * @param tstep  Time value to advance the surface coverages
      * @param rtol   The relative tolerance for the integrator
      * @param atol   The absolute tolerance for the integrator
-     * @param maxStepSize   The maximum step-size the integrator is allowed to take
-     * @param maxSteps   the maximum number of time-steps the integrator can take
-     *                   before reaching tstep
+     * @param maxStepSize   The maximum step-size the integrator is allowed to take.
+     *                      If zero, this option is disabled.
+     * @param maxSteps   The maximum number of time-steps the integrator can take.
+     *                   If not supplied, uses the default value in CVODES (500).
      * @param maxErrTestFails   the maximum permissible number of error test failures
+     *                           If not supplied, uses the default value in CVODES (7).
      */
-    void advanceCoverages(doublereal tstep, doublereal rtol=1.e-7,
-                          doublereal atol=1.e-14, doublereal maxStepSize=0,
-                          size_t maxSteps=0, size_t maxErrTestFails=0);
+    void advanceCoverages(doublereal tstep, double rtol=1.e-7,
+                          double atol=1.e-14, double maxStepSize=0,
+                          size_t maxSteps=500, size_t maxErrTestFails=7);
 
     //! Solve for the pseudo steady-state of the surface problem
     /*!

interfaces/cython/cantera/_cantera.pxd

@@ -379,7 +379,7 @@ cdef extern from "cantera/kinetics/Kinetics.h" namespace "Cantera":
 
 cdef extern from "cantera/kinetics/InterfaceKinetics.h":
     cdef cppclass CxxInterfaceKinetics "Cantera::InterfaceKinetics":
-        void advanceCoverages(double) except +translate_exception
+        void advanceCoverages(double, double, double, double, size_t, size_t) except +translate_exception
         void solvePseudoSteadyStateProblem() except +translate_exception
 
 

interfaces/cython/cantera/kinetics.pyx

@@ -376,13 +376,16 @@ cdef class InterfaceKinetics(Kinetics):
             self._phase_indices[phase.name] = i
             self._phase_indices[i] = i
 
-    def advance_coverages(self, double dt):
+    def advance_coverages(self, double dt, double rtol=1e-7, double atol=1e-14,
+                          double max_step_size=0, int max_steps=500,
+                          int max_error_test_failures=7):
         """
         This method carries out a time-accurate advancement of the surface
         coverages for a specified amount of time.
         """
-        (<CxxInterfaceKinetics*>self.kinetics).advanceCoverages(dt)
-
+        (<CxxInterfaceKinetics*>self.kinetics).advanceCoverages(
+            dt, rtol, atol, max_step_size, max_steps, max_error_test_failures)
+
     def advance_coverages_to_steady_state(self):
         """
         This method advances the surface coverages to steady state.

interfaces/cython/cantera/test/test_kinetics.py

@@ -648,7 +648,7 @@ def cathode_curr(E):
             p.TP = T, P
 
         for s in [anode_surf, oxide_surf_a, cathode_surf, oxide_surf_c]:
-            s.advance_coverages(50.0)
+            s.advance_coverages(50.0, max_steps=20000)
 
         # These values are just a regression test with no theoretical basis
         self.assertArrayNear(anode_surf.coverages,

interfaces/cython/cantera/test/test_onedim.py

@@ -896,7 +896,7 @@ def run_reacting_surface(self, xch4, tsurf, mdot, width):
 
         # integrate the coverage equations holding the gas composition fixed
         # to generate a good starting estimate for the coverages.
-        surf_phase.advance_coverages(1.0)
+        surf_phase.advance_coverages(1.0, max_steps=20000)
 
         sim = ct.ImpingingJet(gas=gas, width=width, surface=surf_phase)
         sim.set_refine_criteria(10.0, 0.3, 0.4, 0.0)

interfaces/cython/cantera/test/test_reactor.py

@@ -1542,3 +1542,42 @@ def test_ConstPressureReactor(self):
         self.assertEqual(states.X[-2], 1)
         for i in range(3,7):
             self.assertNear(states.T[i], states.T[2])
+
+
+class AdvanceCoveragesTest(utilities.CanteraTest):
+    def setup(self, model='ptcombust.xml', gas_phase='gas',
+              interface_phase='Pt_surf'):
+        # create gas and interface
+        self.gas = ct.Solution('ptcombust.xml', 'gas')
+        self.surf = ct.Interface('ptcombust.xml', 'Pt_surf', [self.gas])
+
+    def test_advance_coverages_parameters(self):
+        # create gas and interface
+        self.setup()
+
+        # first, test max step size & max steps
+        dt = 1.0
+        max_steps = 10
+        max_step_size = dt / (max_steps + 1)
+        # this should throw an error, as we can't reach dt
+        with self.assertRaises(ct.CanteraError):
+            self.surf.advance_coverages(
+                dt=dt, max_step_size=max_step_size, max_steps=max_steps)
+
+        # next, run with different tolerances
+        self.setup()
+        self.surf.coverages = 'O(S):0.1, PT(S):0.5, H(S):0.4'
+        self.gas.TP = self.surf.TP
+
+        self.surf.advance_coverages(dt=dt, rtol=1e-5, atol=1e-12,
+                                    max_steps=20000)
+        cov = self.surf.coverages[:]
+
+        self.surf.coverages = 'O(S):0.1, PT(S):0.5, H(S):0.4'
+        self.gas.TP = self.surf.TP
+        self.surf.advance_coverages(dt=dt, rtol=1e-7, atol=1e-14,
+                                    max_steps=20000)
+
+        # check that the solutions are similar, but not identical
+        self.assertArrayNear(cov, self.surf.coverages)
+        self.assertTrue(any(cov != self.surf.coverages))

src/kinetics/ImplicitSurfChem.cpp

@@ -18,8 +18,8 @@ namespace Cantera
 {
 
 ImplicitSurfChem::ImplicitSurfChem(
-        vector<InterfaceKinetics*> k, doublereal rtol, doublereal atol,
-        doublereal maxStepSize, size_t maxSteps,
+        vector<InterfaceKinetics*> k, double rtol, double atol,
+        double maxStepSize, size_t maxSteps,
         size_t maxErrTestFails) :
     m_nv(0),
     m_numTotalBulkSpecies(0),
@@ -109,28 +109,49 @@ void ImplicitSurfChem::getState(doublereal* c)
     }
 }
 
-void ImplicitSurfChem::initialize(doublereal t0)
+void ImplicitSurfChem::setMaxStepSize(double maxstep)
 {
-    m_integ->setTolerances(m_rtol, m_atol);
-    if (m_maxstep > 0)
-    {
+    m_maxstep = maxstep;
+    if (m_maxstep > 0) {
         m_integ->setMaxStepSize(m_maxstep);
     }
-    if (m_nmax > 0)
-    {
-        m_integ->setMaxSteps(m_nmax);
-    }
-    if (m_maxErrTestFails > 0)
-    {
-        m_integ->setMaxErrTestFails(m_maxErrTestFails);
-    }
+}
+
+void ImplicitSurfChem::setTolerances(double rtol, double atol)
+{
+    m_rtol = rtol;
+    m_atol = atol;
+    m_integ->setTolerances(m_rtol, m_atol);
+}
+
+void ImplicitSurfChem::setMaxSteps(size_t maxsteps)
+{
+    m_nmax = maxsteps;
+    m_integ->setMaxSteps(m_nmax);
+}
+
+void ImplicitSurfChem::setMaxErrTestFails(size_t maxErrTestFails)
+{
+    m_maxErrTestFails = maxErrTestFails;
+    m_integ->setMaxErrTestFails(m_maxErrTestFails);
+}
+
+void ImplicitSurfChem::initialize(doublereal t0)
+{
+    this->setTolerances(m_rtol, m_atol);
+    this->setMaxStepSize(m_maxstep);
+    this->setMaxSteps(m_nmax);
+    this->setMaxErrTestFails(m_maxErrTestFails);
     m_integ->initialize(t0, *this);
 }
 
 void ImplicitSurfChem::integrate(doublereal t0, doublereal t1)
 {
-    m_integ->initialize(t0, *this);
-    m_integ->setMaxStepSize(t1 - t0);
+    this->initialize(t0);
+    if (fabs(t1 - t0) < m_maxstep || m_maxstep == 0) {
+        // limit max step size on this run to t1 - t0
+        m_integ->setMaxStepSize(t1 - t0);
+    }
     m_integ->integrate(t1);
     updateState(m_integ->solution());
 }

src/kinetics/InterfaceKinetics.cpp

@@ -773,10 +773,12 @@ void InterfaceKinetics::advanceCoverages(doublereal tstep, doublereal rtol,
 {
     if (m_integrator == 0) {
         vector<InterfaceKinetics*> k{this};
-        m_integrator = new ImplicitSurfChem(k, rtol, atol, maxStepSize, maxSteps,
-                                            maxErrTestFails);
-        m_integrator->initialize();
+        m_integrator = new ImplicitSurfChem(k);
     }
+    m_integrator->setTolerances(rtol, atol);
+    m_integrator->setMaxStepSize(maxStepSize);
+    m_integrator->setMaxSteps(maxSteps);
+    m_integrator->setMaxErrTestFails(maxErrTestFails);
     m_integrator->integrate(0.0, tstep);
     delete m_integrator;
     m_integrator = 0;