Timestep regridding for 1D solver

include/cantera/numerics/SteadyStateSystem.h

@@ -6,6 +6,8 @@
 #ifndef CT_STEADYSTATESYSTEM_H
 #define CT_STEADYSTATESYSTEM_H
 
+#include <cmath>
+
 #include "cantera/base/ct_defs.h"
 #include "SystemJacobian.h"
 
@@ -73,6 +75,11 @@ class SteadyStateSystem
     //! x. On return, array r contains the steady-state residual values.
     double ssnorm(span<const double> x, span<double> r);
 
+    //! Transient max norm (infinity norm) of the residual evaluated using solution
+    //! x and the current timestep (rdt). On return, array r contains the
+    //! transient residual values.
+    double tsnorm(span<const double> x, span<double> r);
+
     //! Total solution vector length;
     size_t size() const {
         return m_size;
@@ -194,6 +201,75 @@ class SteadyStateSystem
         m_tfactor = tfactor;
     }
 
+    //! Set the growth factor used after successful timesteps when the Jacobian is
+    //! re-used.
+    //!
+    //! When adaptive growth is disabled (default), this fixed factor is applied after
+    //! each successful step that does not require a new Jacobian.
+    //!
+    //! When adaptive growth is enabled, this value is used as:
+    //! - the accepted growth factor for heuristics 1, 2, and 4; and
+    //! - the maximum allowed growth factor for heuristic 3.
+    //!
+    //! The default value is 1.5, matching historical behavior.
+    //! @param tfactor  Finite growth factor applied to successful timesteps;
+    //!     must be >= 1.0.
+    void setTimeStepGrowthFactor(double tfactor) {
+        if (!std::isfinite(tfactor) || tfactor < 1.0) {
+            throw CanteraError("SteadyStateSystem::setTimeStepGrowthFactor",
+                "Time step growth factor must be finite and >= 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.
+    //!
+    //! Adaptive heuristics are only applied after successful timesteps that reuse the
+    //! Jacobian. If disabled, the fixed growth factor from setTimeStepGrowthFactor()
+    //! is used directly. Disabled by default.
+    //!
+    //! @param enabled  Enable (`true`) or disable (`false`) adaptive growth heuristics.
+    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 the adaptive time-step growth heuristic.
+    //!
+    //! Heuristic options:
+    //! - `1`: Increase only if the steady-state residual decreases
+    //!   (`ssnorm(x_after) < ssnorm(x_before)`).
+    //! - `2`: Increase only if the transient residual decreases
+    //!   (`tsnorm(x_after) < tsnorm(x_before)`). This is the default.
+    //! - `3`: Scale growth by residual improvement ratio based on transient residual:
+    //!   @f$ \min(m\_tstep\_growth, \max(1, (ts\_before/ts\_after)^{0.2})) @f$.
+    //! - `4`: Increase only if the most recent Newton solve used at most 3 iterations.
+    //!
+    //! @param heuristic  Integer in [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 (`1`-`4`).
+    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;
@@ -255,6 +331,12 @@ class SteadyStateSystem
     //! @param[in] x  Current state vector, length size()
     void evalSSJacobian(span<const double> x);
 
+    //! Determine the timestep growth factor after a successful step.
+    //!
+    //! Called only when a successful step reuses the current Jacobian.
+    double timeStepIncreaseFactor(span<const double> x_before,
+                                  span<const double> x_after);
+
     //! Array of number of steps to take after each unsuccessful steady-state solve
     //! before re-attempting the steady-state solution. For subsequent time stepping
     //! calls, the final value is reused. See setTimeStep().
@@ -267,6 +349,21 @@ class SteadyStateSystem
     //! Factor time step is multiplied by if time stepping fails ( < 1 )
     double m_tfactor = 0.5;
 
+    //! Growth factor for successful steps that reuse the Jacobian.
+    //!
+    //! Used directly when adaptive growth is disabled, and as the base / cap value
+    //! when adaptive growth is enabled.
+    double m_tstep_growth = 1.5;
+
+    //! If `true`, use heuristic gating for successful-step growth; otherwise use the
+    //! fixed growth factor.
+    bool m_adaptive_tstep_growth = false;
+
+    //! Selected adaptive growth heuristic:
+    //! `1`: steady residual gate; `2`: transient residual gate;
+    //! `3`: residual-ratio scaling; `4`: Newton-iteration gate.
+    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
@@ -314,6 +411,7 @@ class SteadyStateSystem
     double m_jacobianRelPerturb = 1e-5;
     //! Absolute perturbation of each component in finite difference Jacobian
     double m_jacobianAbsPerturb = 1e-10;
+
 };
 
 }

include/cantera/oneD/MultiNewton.h

@@ -35,6 +35,11 @@ class MultiNewton
         return m_n;
     }
 
+    //! Number of Newton iterations taken in the most recent solve() call.
+    int lastIterations() const {
+        return m_lastIterations;
+    }
+
     //! Compute the undamped Newton step.  The residual function is evaluated
     //! at `x`, but the Jacobian is not recomputed.
     //! @since Starting in %Cantera 3.2, the Jacobian is accessed via the OneDim object.
@@ -167,13 +172,16 @@ class MultiNewton
     double m_dampFactor = sqrt(2.0);
 
     //! Maximum number of damping iterations
-    size_t m_maxDampIter = 7;
+    size_t m_maxDampIter = 14;
 
     //! number of variables
     size_t m_n;
 
     //! Elapsed CPU time spent computing the Jacobian.
     double m_elapsed = 0.0;
+
+    //! Number of Newton iterations taken in the last solve() call.
+    int m_lastIterations = 0;
 };
 }
 

include/cantera/oneD/Sim1D.h

@@ -337,6 +337,25 @@ class Sim1D : public OneDim
         m_steady_callback = callback;
     }
 
+    //! Set the maximum number of regrid attempts after a timestep failure.
+    //!
+    //! This fallback is used during solve(loglevel, refine_grid=true). Set to `0`
+    //! to disable regrid-on-timestep-failure retries.
+    //!
+    //! @param nmax  Maximum retry attempts; must be >= 0.
+    void setTimeStepRegridMax(int nmax) {
+        if (nmax < 0) {
+            throw CanteraError("Sim1D::setTimeStepRegridMax",
+                "Time step regrid retry count must be >= 0. Got {}.", nmax);
+        }
+        m_ts_regrid_max = nmax;
+    }
+
+    //! Get the maximum number of regrid attempts after a timestep failure.
+    int timeStepRegridMax() const {
+        return m_ts_regrid_max;
+    }
+
 protected:
     //! the solution vector after the last successful steady-state solve (stored
     //! before grid refinement)
@@ -349,6 +368,9 @@ class Sim1D : public OneDim
     //! User-supplied function called after a successful steady-state solve.
     Func1* m_steady_callback;
 
+    //! 0 disables regrid-on-timestep-failure retries
+    int m_ts_regrid_max = 10;
+
 private:
     //! Calls method _finalize in each domain.
     void finalize();

interfaces/cython/cantera/_onedim.pxd

@@ -149,8 +149,16 @@ cdef extern from "cantera/oneD/Sim1D.h":
         void getResidual(double, span[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) except +translate_exception
+        int timeStepRegridMax()
         void setMaxGridPoints(int, size_t) except +translate_exception
         size_t maxGridPoints(size_t) except +translate_exception
         void setGridMin(int, double) except +translate_exception

interfaces/cython/cantera/_onedim.pyi

@@ -297,6 +297,7 @@ class Sim1D:
     def max_time_step_count(self) -> int: ...
     @max_time_step_count.setter
     def max_time_step_count(self, nmax: int) -> None: ...
+    def set_time_step_regrid(self, max_tries: int = 10) -> None: ...
     def set_jacobian_perturbation(
         self, relative: float, absolute: float, threshold: float
     ) -> None: ...
@@ -327,6 +328,13 @@ class Sim1D:
     ) -> None: ...
     def set_max_jac_age(self, ss_age: int, ts_age: int) -> None: ...
     def set_time_step_factor(self, tfactor: float) -> None: ...
+    def set_time_step_growth_factor(self, tfactor: float) -> None: ...
+    def time_step_growth_factor(self) -> float: ...
+    def set_adaptive_time_step_growth(self, enabled: bool = True) -> None: ...
+    def adaptive_time_step_growth(self) -> bool: ...
+    def set_time_step_growth_heuristic(self, heuristic: int) -> None: ...
+    def time_step_growth_heuristic(self) -> int: ...
+    def time_step_regrid(self) -> int: ...
     def set_min_time_step(self, tsmin: float) -> None: ...
     def set_max_time_step(self, tsmax: float) -> None: ...
     @property

interfaces/cython/cantera/_onedim.pyx

@@ -1425,11 +1425,107 @@ 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.
+
+        :param tfactor:
+            Multiplicative reduction factor applied after failed steps.
         """
         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.
+
+        If adaptive growth is disabled (default), this fixed factor is applied
+        directly after successful steps that reuse the Jacobian.
+
+        If adaptive growth is enabled, this factor is used as the accepted
+        growth value for heuristics 1, 2, and 4, and as the maximum growth for
+        heuristic 3.
+
+        :param tfactor:
+            Finite growth factor >= 1.0. The default value is 1.5.
+        """
+        self.sim.setTimeStepGrowthFactor(tfactor)
+
+    def time_step_growth_factor(self):
+        """
+        Get the configured growth factor for successful steps that reuse the
+        Jacobian.
+        """
+        return self.sim.timeStepGrowthFactor()
+
+    def set_adaptive_time_step_growth(self, enabled=True):
+        """
+        Enable or disable adaptive time-step growth heuristics.
+
+        Adaptive heuristics are only used after successful steps that reuse the
+        Jacobian. Disabled by default to preserve legacy behavior.
+
+        :param enabled:
+            `True` to enable adaptive growth heuristics, `False` to use the
+            fixed growth factor.
+        """
+        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.
+
+        Available options are:
+
+        ``1``:
+            Steady-state residual gate. Increase the timestep only if the
+            steady-state residual decreases.
+        ``2``:
+            Transient residual gate (default). Increase the timestep only if
+            the transient residual decreases.
+        ``3``:
+            Residual-ratio scaling. Scale growth using the transient residual
+            improvement ratio, capped by ``time_step_growth_factor()``.
+        ``4``:
+            Newton-iteration gate. Increase only if the most recent Newton
+            solve took at most 3 iterations.
+
+        :param heuristic:
+            Integer in the range [1, 4].
+        """
+        self.sim.setTimeStepGrowthHeuristic(heuristic)
+
+    def time_step_growth_heuristic(self):
+        """
+        Get the selected adaptive time-step growth heuristic (1-4).
+        """
+        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.
+
+        This fallback is used by :meth:`Sim1D.solve` when ``refine_grid=True``.
+
+        :param max_tries:
+            Maximum retry attempts. Set to zero to disable regrid-on-failure
+            retries. Values less than zero are invalid.
+        """
+        self.sim.setTimeStepRegridMax(max_tries)
+
+    def time_step_regrid(self):
+        """
+        Get the maximum number of regrid attempts after a time step failure.
+
+        The default value is 10.
+        """
+        return self.sim.timeStepRegridMax()
+
     def set_min_time_step(self, tsmin):
         """ Set the minimum time step. """
         self.sim.setMinTimeStep(tsmin)