Timestep regridding for 1D solver

include/cantera/numerics/SteadyStateSystem.h

@@ -72,6 +72,11 @@ class SteadyStateSystem
     //! x. On return, array r contains the steady-state residual values.
     double ssnorm(double* x, 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(double* x, double* r);
+
     //! Total solution vector length;
     size_t size() const {
         return m_size;
@@ -204,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;
@@ -266,6 +319,9 @@ class SteadyStateSystem
     //! @param[out] rsd  Storage for the residual, length size()
     void evalSSJacobian(double* x, double* rsd);
 
+    //! Determine the timestep growth factor after a successful step.
+    double timeStepIncreaseFactor(double* x_before, 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().
@@ -278,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
@@ -325,6 +392,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.
@@ -164,13 +169,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,17 @@ class Sim1D : public OneDim
         m_steady_callback = callback;
     }
 
+    //! Set the maximum number of regrid attempts after a timestep failure.
+    //! Set to 0 to disable regrid-on-timestep-failure retries.
+    void setTimeStepRegridMax(int 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 +360,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, 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)
+        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: ...

interfaces/cython/cantera/_onedim.pyx

@@ -1420,11 +1420,62 @@ 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.
+
+        Setting ``max_tries`` to zero disables regrid-on-timestep-failure retries.
+        """
+        self.sim.setTimeStepRegridMax(max_tries)
+
     def set_min_time_step(self, tsmin):
         """ Set the minimum time step. """
         self.sim.setMinTimeStep(tsmin)

src/numerics/SteadyStateSystem.cpp

@@ -108,6 +108,50 @@ void SteadyStateSystem::solve(int loglevel)
     }
 }
 
+double SteadyStateSystem::timeStepIncreaseFactor(double* x_before, double* x_after)
+{
+    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 = m_tstep_growth;
+
+    if (heuristic == 1) {
+        // Steady-state residual gate. If the steady-state residual decreased, allow
+        // increase in timestep.
+        double ss_before = ssnorm(x_before, m_work1.data());
+        double ss_after = ssnorm(x_after, m_work1.data());
+        return (ss_after < ss_before) ? grow_factor : 1.0;
+    } else if (heuristic == 2) {
+        // Transient residual gate. If the transient residual decreased, allow
+        // increase in timestep.
+        double ts_before = tsnorm(x_before, m_work1.data());
+        double ts_after = tsnorm(x_after, m_work1.data());
+        return (ts_after < ts_before) ? grow_factor : 1.0;
+    } else if (heuristic == 3) {
+        // Residual-ratio scaling gate (transient residual). If the transient residual
+        // decreased significantly, allow larger increase in timestep.
+        double ts_before = tsnorm(x_before, m_work1.data());
+        double ts_after = tsnorm(x_after, m_work1.data());
+        if (!(ts_after > 0.0) || !(ts_before > ts_after)) {
+            return 1.0;
+        }
+        const double max_growth = m_tstep_growth;
+        const double exponent = 0.2;
+        double ratio = ts_before / ts_after;
+        double factor = std::pow(ratio, exponent);
+        return std::min(max_growth, std::max(1.0, factor));
+    } else if (heuristic == 4) {
+        // Newton-iteration gate. If the last Newton solve required only a few
+        // iterations, allow increase in timestep.
+        const int max_iters_for_growth = 3;
+        return (newton().lastIterations() <= max_iters_for_growth) ? grow_factor : 1.0;
+    }
+    return m_tstep_growth;
+}
+
 double SteadyStateSystem::timeStep(int nsteps, double dt, double* x, double* r, int loglevel)
 {
     // set the Jacobian age parameter to the transient value
@@ -123,12 +167,14 @@ double SteadyStateSystem::timeStep(int nsteps, double dt, double* x, double* r,
         writelog("============================");
     }
     while (n < nsteps) {
-        if (loglevel == 1) { // At level 1, output concise information
-            double ss = ssnorm(x, r);
-            writelog("\n{:<5d}  {:<6.4e}   {:>7.4f}", n, dt, log10(ss));
-        } else if (loglevel > 1) {
-            double ss = ssnorm(x, r);
-            writelog("\nTimestep ({}) dt= {:<11.4e}  log(ss)= {:<7.4f}", n, dt, log10(ss));
+        if (loglevel >= 1) {
+            double ss_before = ssnorm(x, r);
+            if (loglevel == 1) { // At level 1, output concise information
+                writelog("\n{:<5d}  {:<6.4e}   {:>7.4f}", n, dt, log10(ss_before));
+            } else {
+                writelog("\nTimestep ({}) dt= {:<11.4e}  log(ss)= {:<7.4f}", n, dt,
+                         log10(ss_before));
+            }
         }
 
         // set up for time stepping with stepsize dt
@@ -148,17 +194,18 @@ double SteadyStateSystem::timeStep(int nsteps, double dt, double* x, double* r,
             successiveFailures = 0;
             m_nsteps++;
             n += 1;
-            copy(r, r + m_size, x);
-            // No Jacobian evaluations were performed, so a larger timestep can be used
+            double grow_factor = 1.0;
             if (m_jac->nEvals() == j0) {
-                dt *= 1.5;
+                grow_factor = timeStepIncreaseFactor(x, r);
             }
+            copy(r, r + m_size, x);
+            dt *= grow_factor;
             if (m_time_step_callback) {
                 m_time_step_callback->eval(dt);
             }
             dt = std::min(dt, m_tmax);
             if (m_nsteps >= m_nsteps_max) {
-                throw CanteraError("OneDim::timeStep",
+                throw CanteraError("SteadyStateSystem::timeStep",
                     "Took maximum number of timesteps allowed ({}) without "
                     "reaching steady-state solution.", m_nsteps_max);
             }
@@ -181,7 +228,7 @@ double SteadyStateSystem::timeStep(int nsteps, double dt, double* x, double* r,
                 if (dt < m_tmin) {
                     string err_msg = fmt::format(
                         "Time integration failed. Minimum timestep ({}) reached.", m_tmin);
-                    throw CanteraError("OneDim::timeStep", err_msg);
+                    throw CanteraError("SteadyStateSystem::timeStep", err_msg);
                 }
             }
         }
@@ -212,6 +259,16 @@ double SteadyStateSystem::ssnorm(double* x, double* r)
     return ss;
 }
 
+double SteadyStateSystem::tsnorm(double* x, double* r)
+{
+    eval(x, r, m_rdt, 0);
+    double ts = 0.0;
+    for (size_t i = 0; i < m_size; i++) {
+        ts = std::max(fabs(r[i]), ts);
+    }
+    return ts;
+}
+
 void SteadyStateSystem::setTimeStep(double stepsize, size_t n, const int* tsteps)
 {
     m_tstep = stepsize;
@@ -273,7 +330,9 @@ void SteadyStateSystem::setSteadyMode()
     }
 
     m_rdt = 0.0;
-    m_jac->updateTransient(m_rdt, m_mask.data());
+    if (m_jac_ok) {
+        m_jac->updateTransient(m_rdt, m_mask.data());
+    }
 }
 
 void SteadyStateSystem::setJacAge(int ss_age, int ts_age)

src/oneD/MultiNewton.cpp

@@ -206,11 +206,13 @@ int MultiNewton::solve(double* x0, double* x1, SteadyStateSystem& r, int logleve
     double s1=1.e30;
 
     copy(x0, x0 + m_n, &m_x[0]);
+    m_lastIterations = 0;
 
     double rdt = r.rdt();
     int nJacReeval = 0;
     auto jac = r.linearSolver();
     while (true) {
+        m_lastIterations++;
         // Check whether the Jacobian should be re-evaluated.
         if (jac->age() > m_maxAge) {
             if (loglevel > 1) {