Skip to content

advance limit handling to improve performance (#1453) #1976

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

advance limit handling to improve performance (#1453) #1976

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 16 additions & 0 deletions include/cantera/zeroD/ReactorNet.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -386,6 +393,7 @@ class ReactorNet : public FuncEval
//! and deliberately not exposed in external interfaces.
virtual int lastOrder() const;


vector<Reactor*> m_reactors;
map<string, int> m_counts; //!< Map used for default name generation
unique_ptr<Integrator> m_integ;
Expand Down Expand Up @@ -426,6 +434,14 @@ class ReactorNet : public FuncEval
vector<double> m_ydot;
vector<double> m_yest;
vector<double> m_advancelimits;
//! Base state used for evaluating advance limits during a single advance
//! call when root-finding is enabled
vector<double> 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<double> m_LHS;
Expand Down
38 changes: 34 additions & 4 deletions src/numerics/CVodesIntegrator.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 <iostream>
Expand Down Expand Up @@ -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<Cantera::FuncEval*>(user_data);
// Only ReactorNet implements the limit-check root function
if (auto* net = dynamic_cast<Cantera::ReactorNet*>(f)) {
return net->advanceLimitRootFunc(t, NV_DATA_S(y), gout);
}
// Default: no root
gout[0] = 1.0;
return 0;
Comment on lines +98 to +106
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I suspect you can avoid breaking the abstraction between CVodesIntegrator and ReactorNet by making the call to CVodeRootInit optional, i.e. triggered only if a call to advance with limits applied is made. The CVodes documentation seems to suggest that repeated calls to CVodeRootInit changing the number of root functions to evaluate (including setting it to zero) are allowed.

}
}

CVodesIntegrator::CVodesIntegrator()
Expand Down Expand Up @@ -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(),
Expand Down Expand Up @@ -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();
Expand All @@ -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;
Expand All @@ -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;
Expand Down
108 changes: 77 additions & 31 deletions src/zeroD/ReactorNet.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Comment on lines +271 to +272
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the current test failure (integrator hanging) stems from this. There is a tricky issue with both CVodesIntegrator and ReactorNet having a m_time member, where the former is shadowed by the latter. What I think is happening here is that only the former is being correctly updated to the user-specified time with the interpolated solution while the latter ends up having the time reached by the integrator internally on the last step.

}

double ReactorNet::step()
Expand Down Expand Up @@ -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",
Expand Down
Loading