diff --git a/include/cantera/zeroD/ReactorNet.h b/include/cantera/zeroD/ReactorNet.h index 25bce85535b..e1a8d8c3fde 100644 --- a/include/cantera/zeroD/ReactorNet.h +++ b/include/cantera/zeroD/ReactorNet.h @@ -363,6 +363,13 @@ class ReactorNet : public FuncEval //! @param settings the settings map propagated to all reactors and kinetics objects virtual void setDerivativeSettings(AnyMap& settings); + //! Root function used by the integrator to detect advance limit crossings + //! during calls to advance(t, applylimit=True). Returns 0 on success. + //! The function evaluates gout[0] = 1 - max_i(|y[i]-y_base[i]|/limit[i]). + //! If advance limits are disabled or the limit check is inactive, gout[0] + //! is set to a positive value so that no root is detected. + int advanceLimitRootFunc(double t, const double* y, double* gout); + protected: //! Add the reactor *r* to this reactor network. //! @since Changed in %Cantera 3.2. Previous version used a reference. @@ -386,6 +393,7 @@ class ReactorNet : public FuncEval //! and deliberately not exposed in external interfaces. virtual int lastOrder() const; + vector m_reactors; map m_counts; //!< Map used for default name generation unique_ptr m_integ; @@ -426,6 +434,14 @@ class ReactorNet : public FuncEval vector m_ydot; vector m_yest; vector m_advancelimits; + //! Base state used for evaluating advance limits during a single advance + //! call when root-finding is enabled + vector m_ybase; + //! Base time corresponding to m_ybase + double m_ybase_time = 0.0; + //! Indicates whether the advance-limit root check is active for the + //! current call to advance(t, applylimit=True) + bool m_limit_check_active = false; //! m_LHS is a vector representing the coefficients on the //! "left hand side" of each governing equation vector m_LHS; diff --git a/src/numerics/CVodesIntegrator.cpp b/src/numerics/CVodesIntegrator.cpp index 4474074c7a5..98e0a897eb1 100644 --- a/src/numerics/CVodesIntegrator.cpp +++ b/src/numerics/CVodesIntegrator.cpp @@ -4,6 +4,7 @@ // at https://cantera.org/license.txt for license and copyright information. #include "cantera/numerics/CVodesIntegrator.h" +#include "cantera/zeroD/ReactorNet.h" #include "cantera/base/stringUtils.h" #include @@ -90,6 +91,20 @@ extern "C" { FuncEval* f = (FuncEval*) f_data; return f->preconditioner_solve_nothrow(NV_DATA_S(r),NV_DATA_S(z)); } + + // Root function to enforce ReactorNet advance limits via CVODE event detection + static int advance_limit_root(realtype t, N_Vector y, realtype *gout, void *user_data) + { + // user_data points to the FuncEval (e.g., ReactorNet) + auto* f = static_cast(user_data); + // Only ReactorNet implements the limit-check root function + if (auto* net = dynamic_cast(f)) { + return net->advanceLimitRootFunc(t, NV_DATA_S(y), gout); + } + // Default: no root + gout[0] = 1.0; + return 0; + } } CVodesIntegrator::CVodesIntegrator() @@ -323,6 +338,11 @@ void CVodesIntegrator::initialize(double t0, FuncEval& func) flag = CVodeSetUserData(m_cvode_mem, &func); checkError(flag, "initialize", "CVodeSetUserData"); + // Initialize single root function hook; the underlying evaluator decides + // whether a root is active. This is safe for non-ReactorNet use cases. + flag = CVodeRootInit(m_cvode_mem, 1, advance_limit_root); + checkError(flag, "initialize", "CVodeRootInit"); + if (func.nparams() > 0) { sensInit(t0, func); flag = CVodeSetSensParams(m_cvode_mem, func.m_sens_params.data(), @@ -490,6 +510,7 @@ void CVodesIntegrator::integrate(double tout) tout, m_time); } int nsteps = 0; + bool root_triggered = false; while (m_tInteg < tout) { if (nsteps >= m_maxsteps) { string f_errs = m_func->getErrors(); @@ -502,7 +523,7 @@ void CVodesIntegrator::integrate(double tout) nsteps, tout, m_tInteg, f_errs); } int flag = CVode(m_cvode_mem, tout, m_y, &m_tInteg, CV_ONE_STEP); - if (flag != CV_SUCCESS) { + if (flag != CV_SUCCESS && flag != CV_ROOT_RETURN) { string f_errs = m_func->getErrors(); if (!f_errs.empty()) { f_errs = "Exceptions caught during RHS evaluation:\n" + f_errs; @@ -513,18 +534,27 @@ void CVodesIntegrator::integrate(double tout) "Components with largest weighted error estimates:\n{}", flag, m_error_message, f_errs, getErrorInfo(10)); } + if (flag == CV_ROOT_RETURN) { + // Stop early at root (e.g., advance limit reached) + root_triggered = true; + break; + } nsteps++; } - int flag = CVodeGetDky(m_cvode_mem, tout, 0, m_y); + // Align solution time to the requested tout when no root event occurred; this + // avoids small positive overshoots from CV_ONE_STEP, which can otherwise cause + // callers to see non-monotonic target times. + double t_eval = root_triggered ? m_tInteg : tout; + int flag = CVodeGetDky(m_cvode_mem, t_eval, 0, m_y); checkError(flag, "integrate", "CVodeGetDky"); - m_time = tout; + m_time = t_eval; m_sens_ok = false; } double CVodesIntegrator::step(double tout) { int flag = CVode(m_cvode_mem, tout, m_y, &m_tInteg, CV_ONE_STEP); - if (flag != CV_SUCCESS) { + if (flag != CV_SUCCESS && flag != CV_ROOT_RETURN) { string f_errs = m_func->getErrors(); if (!f_errs.empty()) { f_errs = "Exceptions caught during RHS evaluation:\n" + f_errs; diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp index fc1f79c27fd..d0876abeca9 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -209,51 +209,67 @@ double ReactorNet::advance(double time, bool applylimit) reinitialize(); } - if (!applylimit) { - // take full step + if (!applylimit || !hasAdvanceLimits()) { + // No limit enforcement requested; integrate to requested time advance(time); return time; } - if (!hasAdvanceLimits()) { - // take full step - advance(time); - return time; - } + // Enable root-based limit detection for this advance call + // Set the base state to the current state + m_ybase.assign(m_nv, 0.0); + getState(m_ybase.data()); + m_ybase_time = m_time; + m_limit_check_active = true; - getAdvanceLimits(m_advancelimits.data()); + // Integrate toward the requested time; integrator will return early + // if a limit is reached (CV_ROOT_RETURN) + m_integ->integrate(time); - // ensure that gradient is available - while (lastOrder() < 1) { - step(); - } + // Update reactor states to match the integrator solution at the time + // reached (which may be earlier than 'time' if a limit was triggered) + updateState(m_integ->solution()); - int k = lastOrder(); - double t = time, delta; - double* y = m_integ->solution(); + // Disable limit checking after this call + m_limit_check_active = false; - // reduce time step if limits are exceeded - while (true) { - bool exceeded = false; - getEstimate(t, k, &m_yest[0]); + // Verbose logging analogous to prior predictive limiter: when a root event + // stopped integration before reaching the requested time, report the most + // limiting component and details about the step. + if (m_verbose && m_time < time) { + // Ensure limits are available + if (m_advancelimits.size() != m_nv) { + m_advancelimits.assign(m_nv, -1.0); + } + getAdvanceLimits(m_advancelimits.data()); + double* ycurr = m_integ->solution(); + size_t jmax = npos; + double max_ratio = -1.0; + double best_delta = 0.0; + double best_limit = 0.0; for (size_t j = 0; j < m_nv; j++) { - delta = abs(m_yest[j] - y[j]); - if ( (m_advancelimits[j] > 0.) && ( delta > m_advancelimits[j]) ) { - exceeded = true; - if (m_verbose) { - writelog(" Limiting global state vector component {:d} (dt = {:9.4g}):" - "{:11.6g} > {:9.4g}\n", - j, t - m_time, delta, m_advancelimits[j]); + double lim = m_advancelimits[j]; + if (lim > 0.0) { + double delta = abs(ycurr[j] - m_ybase[j]); + double ratio = delta / lim; + if (ratio > max_ratio) { + max_ratio = ratio; + jmax = j; + best_delta = delta; + best_limit = lim; } } } - if (!exceeded) { - break; + if (jmax != npos) { + double dt = m_time - m_ybase_time; + writelog(" Advance limit triggered for component {:d} (dt = {:9.4g}):" + " |Δ| = {:11.6g}, limit = {:9.4g}, ratio = {:9.4g}\n", + jmax, dt, best_delta, best_limit, max_ratio); } - t = .5 * (m_time + t); } - advance(t); - return t; + + // m_time is tracked via callbacks during integration + return m_time; } double ReactorNet::step() @@ -334,6 +350,36 @@ int ReactorNet::lastOrder() const } } +int ReactorNet::advanceLimitRootFunc(double t, const double* y, double* gout) +{ + // Default: no root detected + double g = 1.0; + + if (m_limit_check_active) { + // Ensure limits vector is current + if (m_advancelimits.size() != m_nv) { + m_advancelimits.assign(m_nv, -1.0); + } + getAdvanceLimits(m_advancelimits.data()); + + double max_ratio = 0.0; + for (size_t i = 0; i < m_nv; i++) { + double lim = m_advancelimits[i]; + if (lim > 0.0) { + double delta = abs(y[i] - m_ybase[i]); + double ratio = delta / lim; + if (ratio > max_ratio) { + max_ratio = ratio; + } + } + } + g = 1.0 - max_ratio; // root at g = 0 when any component reaches its limit + } + + gout[0] = g; + return 0; +} + void ReactorNet::addReactor(Reactor& r) { warn_deprecated("ReactorNet::addReactor",