address broken steady state tests, add options for time step regridding controls, add tests for new options

Author: wandadars

include/cantera/numerics/SteadyStateSystem.h

@@ -209,6 +209,54 @@ class SteadyStateSystem
         m_tfactor = tfactor;
     }
 
+    //! Set the default factor by which the time step is increased after a
+    //! successful step when the Jacobian is re-used.
+    //!
+    //! The default value is 1.5, matching historical behavior.
+    //! @param tfactor  factor time step is multiplied by after a successful step
+    void setTimeStepGrowthFactor(double tfactor) {
+        if (tfactor < 1.0) {
+            throw CanteraError("SteadyStateSystem::setTimeStepGrowthFactor",
+                "Time step growth factor must be >= 1.0. Got {}.", tfactor);
+        }
+        m_tstep_growth = tfactor;
+    }
+
+    //! Get the successful-step time step growth factor.
+    double timeStepGrowthFactor() const {
+        return m_tstep_growth;
+    }
+
+    //! Enable or disable adaptive time-step growth heuristics.
+    //!
+    //! If disabled, steady solvers use the fixed growth factor from
+    //! setTimeStepGrowthFactor(). Disabled by default.
+    void setAdaptiveTimeStepGrowth(bool enabled) {
+        m_adaptive_tstep_growth = enabled;
+    }
+
+    //! Get whether adaptive time-step growth heuristics are enabled.
+    bool adaptiveTimeStepGrowth() const {
+        return m_adaptive_tstep_growth;
+    }
+
+    //! Select which adaptive time-step growth heuristic is used.
+    //!
+    //! Supported values are 1-4.
+    void setTimeStepGrowthHeuristic(int heuristic) {
+        if (heuristic < 1 || heuristic > 4) {
+            throw CanteraError("SteadyStateSystem::setTimeStepGrowthHeuristic",
+                "Time step growth heuristic must be in the range [1, 4]. Got {}.",
+                heuristic);
+        }
+        m_tstep_growth_heuristic = heuristic;
+    }
+
+    //! Get the selected adaptive time-step growth heuristic.
+    int timeStepGrowthHeuristic() const {
+        return m_tstep_growth_heuristic;
+    }
+
     //! Set the maximum number of timeteps allowed before successful steady-state solve
     void setMaxTimeStepCount(int nmax) {
         m_nsteps_max = nmax;
@@ -286,6 +334,17 @@ class SteadyStateSystem
     //! Factor time step is multiplied by if time stepping fails ( < 1 )
     double m_tfactor = 0.5;
 
+    //! Factor time step is multiplied by after successful steps when no new Jacobian
+    //! evaluation is needed.
+    double m_tstep_growth = 1.5;
+
+    //! If `true`, use residual- or iteration-based heuristics to gate time-step
+    //! growth; otherwise use the fixed growth factor.
+    bool m_adaptive_tstep_growth = false;
+
+    //! Selected adaptive growth heuristic (1-4). See timeStepIncreaseFactor().
+    int m_tstep_growth_heuristic = 2;
+
     shared_ptr<vector<double>> m_state; //!< Solution vector
 
     //! Work array used to hold the residual or the new solution

interfaces/cython/cantera/_onedim.pxd

@@ -149,6 +149,12 @@ cdef extern from "cantera/oneD/Sim1D.h":
         void getResidual(double, double*) except +translate_exception
         void setJacAge(int, int)
         void setTimeStepFactor(double)
+        void setTimeStepGrowthFactor(double) except +translate_exception
+        double timeStepGrowthFactor()
+        void setAdaptiveTimeStepGrowth(cbool)
+        cbool adaptiveTimeStepGrowth()
+        void setTimeStepGrowthHeuristic(int) except +translate_exception
+        int timeStepGrowthHeuristic()
         void setMinTimeStep(double)
         void setMaxTimeStep(double)
         void setTimeStepRegridMax(int)

interfaces/cython/cantera/_onedim.pyx

@@ -1420,11 +1420,54 @@ cdef class Sim1D:
 
     def set_time_step_factor(self, tfactor):
         """
-        Set the factor by which the time step will be increased after a
-        successful step, or decreased after an unsuccessful one.
+        Set the factor by which the time step will be decreased after an
+        unsuccessful step.
         """
         self.sim.setTimeStepFactor(tfactor)
 
+    def set_time_step_growth_factor(self, tfactor):
+        """
+        Set the factor by which the time step will be increased after a
+        successful step when the Jacobian is reused.
+
+        The default value is 1.5.
+        """
+        self.sim.setTimeStepGrowthFactor(tfactor)
+
+    def time_step_growth_factor(self):
+        """
+        Get the factor used to increase time step size after successful steps.
+        """
+        return self.sim.timeStepGrowthFactor()
+
+    def set_adaptive_time_step_growth(self, enabled=True):
+        """
+        Enable adaptive time-step growth heuristics.
+
+        Disabled by default to preserve legacy behavior.
+        """
+        self.sim.setAdaptiveTimeStepGrowth(enabled)
+
+    def adaptive_time_step_growth(self):
+        """
+        Get whether adaptive time-step growth heuristics are enabled.
+        """
+        return self.sim.adaptiveTimeStepGrowth()
+
+    def set_time_step_growth_heuristic(self, heuristic):
+        """
+        Select the adaptive time-step growth heuristic.
+
+        Supported values are integers from 1 to 4.
+        """
+        self.sim.setTimeStepGrowthHeuristic(heuristic)
+
+    def time_step_growth_heuristic(self):
+        """
+        Get the selected adaptive time-step growth heuristic.
+        """
+        return self.sim.timeStepGrowthHeuristic()
+
     def set_time_step_regrid(self, max_tries=10):
         """
         Set the maximum number of regrid attempts after a time step failure.

src/numerics/SteadyStateSystem.cpp

@@ -110,9 +110,13 @@ void SteadyStateSystem::solve(int loglevel)
 
 double SteadyStateSystem::timeStepIncreaseFactor(double* x_before, double* x_after)
 {
-    int heuristic = 2; // Can be 1-4
+    if (!m_adaptive_tstep_growth) {
+        return m_tstep_growth;
+    }
+
+    int heuristic = m_tstep_growth_heuristic; // Can be 1-4
     m_work1.resize(m_size);
-    const double grow_factor = 1.25;
+    const double grow_factor = m_tstep_growth;
 
     if (heuristic == 1) {
         // Steady-state residual gate. If the steady-state residual decreased, allow
@@ -134,7 +138,7 @@ double SteadyStateSystem::timeStepIncreaseFactor(double* x_before, double* x_aft
         if (!(ts_after > 0.0) || !(ts_before > ts_after)) {
             return 1.0;
         }
-        const double max_growth = 1.5;
+        const double max_growth = m_tstep_growth;
         const double exponent = 0.2;
         double ratio = ts_before / ts_after;
         double factor = std::pow(ratio, exponent);
@@ -145,7 +149,7 @@ double SteadyStateSystem::timeStepIncreaseFactor(double* x_before, double* x_aft
         const int max_iters_for_growth = 3;
         return (newton().lastIterations() <= max_iters_for_growth) ? grow_factor : 1.0;
     }
-    return 1.0;
+    return m_tstep_growth;
 }
 
 double SteadyStateSystem::timeStep(int nsteps, double dt, double* x, double* r, int loglevel)
@@ -314,7 +318,7 @@ void SteadyStateSystem::initTimeInteg(double dt, double* x)
     m_rdt = 1.0/dt;
 
     // if the stepsize has changed, then update the transient part of the Jacobian
-    if (fabs(rdt_old - m_rdt) > Tiny && m_jac_ok) {
+    if (fabs(rdt_old - m_rdt) > Tiny) {
         m_jac->updateTransient(m_rdt, m_mask.data());
     }
 }

test/python/test_onedim.py

@@ -139,6 +139,32 @@ def test_tolerances(self):
         assert rtol_ss == set((5e-3, 3e-4, 7e-7))
         assert rtol_ts == set((6e-3, 4e-4, 2e-7))
 
+    def test_time_step_growth_controls(self):
+        gas = ct.Solution("h2o2.yaml")
+        left = ct.Inlet1D(gas)
+        flame = ct.FreeFlow(gas)
+        right = ct.Outlet1D(gas)
+        sim = ct.Sim1D((left, flame, right))
+
+        assert sim.time_step_growth_factor() == approx(1.5)
+        assert not sim.adaptive_time_step_growth()
+        assert sim.time_step_growth_heuristic() == 2
+
+        sim.set_time_step_growth_factor(1.2)
+        assert sim.time_step_growth_factor() == approx(1.2)
+
+        sim.set_adaptive_time_step_growth(True)
+        assert sim.adaptive_time_step_growth()
+
+        sim.set_time_step_growth_heuristic(4)
+        assert sim.time_step_growth_heuristic() == 4
+
+        with pytest.raises(ct.CanteraError, match=">= 1.0"):
+            sim.set_time_step_growth_factor(0.9)
+
+        with pytest.raises(ct.CanteraError, match=r"\[1, 4\]"):
+            sim.set_time_step_growth_heuristic(5)
+
     def test_switch_transport(self):
         gas = ct.Solution('h2o2.yaml')
         gas.set_equivalence_ratio(0.9, 'H2', 'O2:0.21, N2:0.79')