GreenPolyProducts/cantera
5
0
Fork 0
mirror of https://github.com/Cantera/cantera.git synced 2025-02-25 18:55:29 -06:00
Code Issues Packages Projects Releases Wiki Activity

[Reactor] Always use "one step" mode with SUNDIALS integrators

Browse Source 
Do the interpolation required for "advance" explicitly, rather
than relying on SUNDIALS to do it as part of the call to CVode() or
IDASolve().

This circumvents a change in behavior introduced in SUNDIALS 6.6.

Fixes #1554.
This commit is contained in:
Ray Speth 2023-08-03 18:45:45 -04:00 committed by Ingmar Schoegl
parent c9a6f818a1
commit ce797baffb
5 changed files with 77 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
include/cantera/numerics/CVodesIntegrator.h
View File

@ -97,8 +97,16 @@ private:
void* m_linsol_matrix = nullptr; //!< matrix used by Sundials
FuncEval* m_func = nullptr;
double m_t0 = 0.0;
double m_time; //!< The current integrator time
//! The current system time, corresponding to #m_y
double m_time;
//! The latest time reached by the integrator. May be greater than #m_time.
double m_tInteg;
//! The system state at #m_time
N_Vector m_y = nullptr;
N_Vector m_abstol = nullptr;
N_Vector m_dky = nullptr;
string m_type = "DENSE";

6
include/cantera/numerics/IdasIntegrator.h
View File

@ -107,7 +107,11 @@ private:
FuncEval* m_func = nullptr;
double m_t0 = 0.0; //!< The start time for the integrator
double m_time; //!< The current integrator time
double m_time; //!< The current integrator time, corresponding to #m_y
//! The latest time reached by the integrator. May be greater than #m_time
double m_tInteg;
N_Vector m_y = nullptr; //!< The current system state
N_Vector m_ydot = nullptr; //!< The time derivatives of the system state
N_Vector m_abstol = nullptr; //!< Absolute tolerances for each state variable

41
src/numerics/CVodesIntegrator.cpp
View File

@ -254,6 +254,7 @@ void CVodesIntegrator::initialize(double t0, FuncEval& func)
m_neq = func.neq();
m_t0 = t0;
m_time = t0;
m_tInteg = t0;
m_func = &func;
func.clearErrors();
// Initialize preconditioner if applied
@ -351,6 +352,7 @@ void CVodesIntegrator::reinitialize(double t0, FuncEval& func)
{
m_t0 = t0;
m_time = t0;
m_tInteg = t0;
func.getState(NV_DATA_S(m_y));
m_func = &func;
func.clearErrors();
@ -513,24 +515,40 @@ void CVodesIntegrator::integrate(double tout)
if (tout == m_time) {
return;
}
int flag = CVode(m_cvode_mem, tout, m_y, &m_time, CV_NORMAL);
if (flag != CV_SUCCESS) {
string f_errs = m_func->getErrors();
if (!f_errs.empty()) {
f_errs = "Exceptions caught during RHS evaluation:\n" + f_errs;
int nsteps = 0;
while (m_tInteg < tout) {
if (nsteps >= m_maxsteps) {
throw CanteraError("CVodesIntegrator::integrate",
"Maximum number of timesteps ({}) taken without reaching output "
"time ({}).\nCurrent integrator time: {}",
nsteps, tout, m_tInteg);
}
throw CanteraError("CVodesIntegrator::integrate",
"CVodes error encountered. Error code: {}\n{}\n"
"{}"
"Components with largest weighted error estimates:\n{}",
flag, m_error_message, f_errs, getErrorInfo(10));
int flag = CVode(m_cvode_mem, tout, m_y, &m_tInteg, CV_ONE_STEP);
if (flag != CV_SUCCESS) {
string f_errs = m_func->getErrors();
if (!f_errs.empty()) {
f_errs = "Exceptions caught during RHS evaluation:\n" + f_errs;
}
throw CanteraError("CVodesIntegrator::integrate",
"CVodes error encountered. Error code: {}\n{}\n"
"{}"
"Components with largest weighted error estimates:\n{}",
flag, m_error_message, f_errs, getErrorInfo(10));
}
nsteps++;
}
int flag = CVodeGetDky(m_cvode_mem, tout, 0, m_y);
if (flag != CV_SUCCESS) {
throw CanteraError("CvodesIntegrator::integrate",
"CVodeGetDky failed. result = {}", flag);
}
m_time = tout;
m_sens_ok = false;
}
double CVodesIntegrator::step(double tout)
{
int flag = CVode(m_cvode_mem, tout, m_y, &m_time, CV_ONE_STEP);
int flag = CVode(m_cvode_mem, tout, m_y, &m_tInteg, CV_ONE_STEP);
if (flag != CV_SUCCESS) {
string f_errs = m_func->getErrors();
if (!f_errs.empty()) {
@ -544,6 +562,7 @@ double CVodesIntegrator::step(double tout)
}
m_sens_ok = false;
m_time = m_tInteg;
return m_time;
}

38
src/numerics/IdasIntegrator.cpp
View File

@ -227,6 +227,7 @@ void IdasIntegrator::initialize(double t0, FuncEval& func)
m_neq = func.neq();
m_t0 = t0;
m_time = t0;
m_tInteg = t0;
m_func = &func;
func.clearErrors();
@ -317,6 +318,7 @@ void IdasIntegrator::reinitialize(double t0, FuncEval& func)
{
m_t0 = t0;
m_time = t0;
m_tInteg = t0;
func.getStateDae(NV_DATA_S(m_y), NV_DATA_S(m_ydot));
m_func = &func;
func.clearErrors();
@ -452,23 +454,36 @@ void IdasIntegrator::integrate(double tout)
if (tout == m_time) {
return;
}
int flag = IDASolve(m_ida_mem, tout, &m_time, m_y, m_ydot, IDA_NORMAL);
if (flag != IDA_SUCCESS) {
string f_errs = m_func->getErrors();
if (!f_errs.empty()) {
f_errs = "Exceptions caught during RHS evaluation:\n" + f_errs;
int nsteps = 0;
while (m_tInteg < tout) {
if (nsteps >= m_maxsteps) {
throw CanteraError("IdasIntegrator::integrate",
"Maximum number of timesteps ({}) taken without reaching output "
"time ({}).\nCurrent integrator time: {}",
nsteps, tout, m_time);
}
throw CanteraError("IdasIntegrator::integrate",
"IDA error encountered. Error code: {}\n{}\n"
"{}"
"Components with largest weighted error estimates:\n{}",
flag, m_error_message, f_errs, getErrorInfo(10));
int flag = IDASolve(m_ida_mem, tout, &m_tInteg, m_y, m_ydot, IDA_ONE_STEP);
if (flag != IDA_SUCCESS) {
string f_errs = m_func->getErrors();
if (!f_errs.empty()) {
f_errs = "Exceptions caught during RHS evaluation:\n" + f_errs;
}
throw CanteraError("IdasIntegrator::integrate",
"IDA error encountered. Error code: {}\n{}\n"
"{}"
"Components with largest weighted error estimates:\n{}",
flag, m_error_message, f_errs, getErrorInfo(10));
}
nsteps++;
}
int flag = IDAGetDky(m_ida_mem, tout, 0, m_y);
checkError(flag, "integrate", "IDAGetDky");
m_time = tout;
}
double IdasIntegrator::step(double tout)
{
int flag = IDASolve(m_ida_mem, tout, &m_time, m_y, m_ydot, IDA_ONE_STEP);
int flag = IDASolve(m_ida_mem, tout, &m_tInteg, m_y, m_ydot, IDA_ONE_STEP);
if (flag != IDA_SUCCESS) {
string f_errs = m_func->getErrors();
if (!f_errs.empty()) {
@ -481,6 +496,7 @@ double IdasIntegrator::step(double tout)
flag, f_errs, m_error_message, getErrorInfo(10));
}
m_time = m_tInteg;
return m_time;
}

11
test/python/test_reactor.py
View File

@ -221,7 +221,7 @@ class TestReactor(utilities.CanteraTest):
self.net.max_time_step = max_step_size
self.net.max_steps = max_steps
with self.assertRaisesRegex(
ct.CanteraError, 'mxstep steps taken before reaching tout'):
ct.CanteraError, 'Maximum number of timesteps'):
self.net.advance(1e-04)
self.assertLessEqual(self.net.time, max_steps * max_step_size)
self.assertEqual(self.net.max_steps, max_steps)
@ -1605,7 +1605,7 @@ class TestFlowReactor2(utilities.CanteraTest):
sim.max_steps = 13
assert sim.max_steps == 13
with pytest.raises(ct.CanteraError, match="mxstep steps taken"):
with pytest.raises(ct.CanteraError, match="Maximum number of timesteps"):
sim.advance(0.1)
assert sim.solver_stats['steps'] == 13
@ -1624,11 +1624,12 @@ class TestFlowReactor2(utilities.CanteraTest):
surf.TP = gas.TP
r, rsurf, sim = self.make_reactors(gas, surf)
sim.max_time_step = 0.2
sim.max_time_step = 0.002
sim.advance(0.01)
sim.step()
x1 = sim.distance
x2 = sim.step()
dx_limit = 0.05 * (x2-x1)
dx_limit = 0.1 * (x2-x1)
sim.max_time_step = dx_limit
assert sim.max_time_step == dx_limit
@ -1746,7 +1747,7 @@ class TestFlowReactor2(utilities.CanteraTest):
# With few steps allowed, won't be able to reach steady-state
r.inlet_surface_max_error_failures = 10
r.inlet_surface_max_steps = 200
with pytest.raises(ct.CanteraError, match="mxstep steps taken"):
with pytest.raises(ct.CanteraError, match="Maximum number of timesteps"):
sim.initialize()
# Relaxing the tolerances will allow the integrator to reach steady-state