From 2f668f28c0f5aa0b9bdd380726d504b7747e0e2d Mon Sep 17 00:00:00 2001 From: Christopher Neal Date: Sun, 11 Jan 2026 03:07:53 -0500 Subject: [PATCH 1/2] first commit for feature --- include/cantera/numerics/SteadyStateSystem.h | 9 +++ include/cantera/oneD/MultiNewton.h | 10 ++- include/cantera/oneD/Sim1D.h | 14 ++++ interfaces/cython/cantera/_onedim.pxd | 2 + interfaces/cython/cantera/_onedim.pyi | 1 + interfaces/cython/cantera/_onedim.pyx | 8 ++ src/numerics/SteadyStateSystem.cpp | 81 ++++++++++++++++---- src/oneD/MultiNewton.cpp | 2 + src/oneD/Sim1D.cpp | 26 ++++++- 9 files changed, 138 insertions(+), 15 deletions(-) diff --git a/include/cantera/numerics/SteadyStateSystem.h b/include/cantera/numerics/SteadyStateSystem.h index 8579578202d..872254870a7 100644 --- a/include/cantera/numerics/SteadyStateSystem.h +++ b/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; @@ -266,6 +271,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(). @@ -325,6 +333,7 @@ class SteadyStateSystem double m_jacobianRelPerturb = 1e-5; //! Absolute perturbation of each component in finite difference Jacobian double m_jacobianAbsPerturb = 1e-10; + }; } diff --git a/include/cantera/oneD/MultiNewton.h b/include/cantera/oneD/MultiNewton.h index 028d16434c1..e0f41d4bb6d 100644 --- a/include/cantera/oneD/MultiNewton.h +++ b/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; }; } diff --git a/include/cantera/oneD/Sim1D.h b/include/cantera/oneD/Sim1D.h index b4999375f0e..a6890ac1fb3 100644 --- a/include/cantera/oneD/Sim1D.h +++ b/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(); diff --git a/interfaces/cython/cantera/_onedim.pxd b/interfaces/cython/cantera/_onedim.pxd index a452d629e85..12b47a15f11 100644 --- a/interfaces/cython/cantera/_onedim.pxd +++ b/interfaces/cython/cantera/_onedim.pxd @@ -151,6 +151,8 @@ cdef extern from "cantera/oneD/Sim1D.h": void setTimeStepFactor(double) 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 diff --git a/interfaces/cython/cantera/_onedim.pyi b/interfaces/cython/cantera/_onedim.pyi index aea93f5b8da..9202ec47f43 100644 --- a/interfaces/cython/cantera/_onedim.pyi +++ b/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: ... diff --git a/interfaces/cython/cantera/_onedim.pyx b/interfaces/cython/cantera/_onedim.pyx index 6f905ae6b2e..64b61fd8ecf 100644 --- a/interfaces/cython/cantera/_onedim.pyx +++ b/interfaces/cython/cantera/_onedim.pyx @@ -1425,6 +1425,14 @@ cdef class Sim1D: """ self.sim.setTimeStepFactor(tfactor) + 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) diff --git a/src/numerics/SteadyStateSystem.cpp b/src/numerics/SteadyStateSystem.cpp index e57f262d8bd..41a9914c659 100644 --- a/src/numerics/SteadyStateSystem.cpp +++ b/src/numerics/SteadyStateSystem.cpp @@ -108,6 +108,46 @@ void SteadyStateSystem::solve(int loglevel) } } +double SteadyStateSystem::timeStepIncreaseFactor(double* x_before, double* x_after) +{ + int heuristic = 2; // Can be 1-4 + m_work1.resize(m_size); + const double grow_factor = 1.25; + + 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 = 1.5; + 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 1.0; +} + double SteadyStateSystem::timeStep(int nsteps, double dt, double* x, double* r, int loglevel) { // set the Jacobian age parameter to the transient value @@ -123,12 +163,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 +190,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 +224,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 +255,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; @@ -261,7 +314,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) { + if (fabs(rdt_old - m_rdt) > Tiny && m_jac_ok) { m_jac->updateTransient(m_rdt, m_mask.data()); } } @@ -273,7 +326,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) diff --git a/src/oneD/MultiNewton.cpp b/src/oneD/MultiNewton.cpp index eb0dd2a2f02..fadf57c6911 100644 --- a/src/oneD/MultiNewton.cpp +++ b/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) { diff --git a/src/oneD/Sim1D.cpp b/src/oneD/Sim1D.cpp index b82560716d1..da4872724c0 100644 --- a/src/oneD/Sim1D.cpp +++ b/src/oneD/Sim1D.cpp @@ -344,8 +344,32 @@ void Sim1D::solve(int loglevel, bool refine_grid) clearDebugFile(); } + int retries = 0; while (new_points > 0) { - SteadyStateSystem::solve(loglevel); + if (refine_grid == true) { + try { + SteadyStateSystem::solve(loglevel); + } catch (CanteraError& err) { + if (err.getMethod() == "SteadyStateSystem::timeStep" && retries < m_ts_regrid_max) { + writelog("\nTime stepping failed; attempting to refine the grid and retry " + "({}/{})...\n", retries+1, m_ts_regrid_max); + new_points = refine(loglevel); + if (new_points > 0) { + retries++; + continue; + } + if (loglevel > 0) { + if (new_points == 0) { + writelog("Regrid retry aborted: no new points added.\n"); + } + } + } + throw; + } + } else { + SteadyStateSystem::solve(loglevel); + } + if (loglevel > 0) { writelog("\nNewton steady-state solve succeeded.\n\n"); writelog("Problem solved on ["); From bbd68d400236a7f1e2704b7ded38ccbcd4daf2a7 Mon Sep 17 00:00:00 2001 From: Christopher Neal Date: Thu, 12 Feb 2026 02:57:53 -0500 Subject: [PATCH 2/2] address broken steady state tests, add options for time step regridding controls, add tests for new options --- include/cantera/numerics/SteadyStateSystem.h | 59 ++++++++++++++++++++ interfaces/cython/cantera/_onedim.pxd | 6 ++ interfaces/cython/cantera/_onedim.pyx | 47 +++++++++++++++- src/numerics/SteadyStateSystem.cpp | 14 +++-- test/python/test_onedim.py | 26 +++++++++ 5 files changed, 145 insertions(+), 7 deletions(-) diff --git a/include/cantera/numerics/SteadyStateSystem.h b/include/cantera/numerics/SteadyStateSystem.h index 872254870a7..2b6ed248558 100644 --- a/include/cantera/numerics/SteadyStateSystem.h +++ b/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> m_state; //!< Solution vector //! Work array used to hold the residual or the new solution diff --git a/interfaces/cython/cantera/_onedim.pxd b/interfaces/cython/cantera/_onedim.pxd index 12b47a15f11..e621a46b95e 100644 --- a/interfaces/cython/cantera/_onedim.pxd +++ b/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) diff --git a/interfaces/cython/cantera/_onedim.pyx b/interfaces/cython/cantera/_onedim.pyx index 64b61fd8ecf..c255feaec64 100644 --- a/interfaces/cython/cantera/_onedim.pyx +++ b/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. diff --git a/src/numerics/SteadyStateSystem.cpp b/src/numerics/SteadyStateSystem.cpp index 41a9914c659..7bfb0357bdb 100644 --- a/src/numerics/SteadyStateSystem.cpp +++ b/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()); } } diff --git a/test/python/test_onedim.py b/test/python/test_onedim.py index 07fbadb1b3e..cd286de0599 100644 --- a/test/python/test_onedim.py +++ b/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')