[Montana]

Fixes #1301.

What was wrong

The issue's two messages are one defect, and the second is not a stale handle.

1. CVODES returned an output time from a step that does not span it. After the sensitivity RHS refuses a state (GH #395's recoverable return), a CV_NORMAL call can return tout from a last step that starts after it. Traced on the issue's MPEKRIS fixture, at row 44: one call from tn = 43.984 took 3 steps, ended at 44.0311 with a last step from 44.0099, and returned flag 0, t_ret = 44. CVode then interpolates to tout with (void) CVodeGetDky(...), ignoring its result. So the refusal ("Illegal value for t") left y as the call before had it, and the row took that state, at the time asked for.

2. The row's sensitivities were left at zero. The per-row CVodeGetSens interpolates at CVODES's last return time, which is that same tout, so it failed too. The output loop had no else, so 24 rows of 101 came back with all-zero sensitivities, and no error. On the plain-run comparison, those are rows 44, 47, 49–52, 56, … 98.

3. The "solver at t = 2.6e-312" is a SUNDIALS message bug. In SUNDIALS 7.2.1, CVodeGetSensDky1 calls cvProcessError(…, MSGCV_BAD_T), a format with three %g fields, with no arguments. CVodeGetDky passes t, tn - hu, tn. The denormals are whatever was in the registers.

Either recovery switch removes it (BNGSIM_SENS_NONFINITE_RECOVER=0, BNGSIM_SENS_CLAMP_NUMERIC_ZERO=0: 0 bad reads instead of 48), which confirms the trigger.

The fix (cvode_simulator.cpp)

1. Detect it. After a CV_NORMAL return of t_ret == tout in a sensitivity run, check whether tout is inside the last step [tn − hu, tn], with the fuzz CVodeGetDky itself allows (inside_last_step).
2. Read the solver where it stands. If tout is outside, take the state at tn itself, which is always readable, and set t_ret = tn. A sensitivity run already records each row at the time it stands at ("A run with sensitivity columns keeps the time it stands at"), so the row is labelled honestly. Each such row is counted in solver_stats["n_rows_outside_last_step"].
3. Read the row's sensitivities at its own time, with CVodeGetSensDky(t_ret) instead of CVodeGetSens, which reads at CVODES's stale return time.
4. Never zero silently. A row whose sensitivities still cannot be read now raises, naming the row and the flag, instead of recording 0.

-                flag = CVode(cvode_mem, look_tout, y, &t_ret,
-                             (one_step || first_alone) ? CV_ONE_STEP : CV_NORMAL);
+                const int itask = (one_step || first_alone) ? CV_ONE_STEP : CV_NORMAL;
+                flag = CVode(cvode_mem, look_tout, y, &t_ret, itask);
+                if (flag == CV_SUCCESS && itask == CV_NORMAL && sens.n_total > 0 &&
+                    t_ret == look_tout) {
+                    ... if (!inside_last_step(look_tout, tn_now, hu_now)) {
+                        CVodeGetDky(cvode_mem, tn_now, 0, y);   // the state where the run stands
+                        t_ret = tn_now;
+                        ++rows_outside_last_step;
+                    }
+                }
 ...
-            flag = CVodeGetSens(cvode_mem, &t_ret, yS_guard.arr);
+            flag = CVodeGetSensDky(cvode_mem, t_ret, 0, yS_guard.arr);
+            if (flag != CV_SUCCESS) { throw ... }   // never the zeros the buffer starts as

CVODES still prints its own CVodeGetDky message from inside CVode, once per such row. bngsim cannot stop it attempting the interpolation, only stop trusting it. The CVodeGetSensDky1 messages are gone.

Verified

MPEKRIS (sensitivity_params=["sigma_nk", "g_reg"], rtol 1e-9):

	main	this PR
rows with all-zero sensitivities (t > 0)	24	0
n_rows_outside_last_step	—	25
SUNDIALS messages	48	25 (CVODES's own, in CVode)

The rows are right, not just non-zero. Each moved row's deviation from a spline through its neighbours is, for sensitivities, 7.2e-4 against 2.0e-6 median and 1.4e-3 max elsewhere; for the state, 1.8e-4 against a 1.4e-3 max. That is within the run's own range, where main's were zeros. By construction a moved row's state is CVODES's solution at tn; on main it was the previous call's state at 43.98, labelled 44.

MODEL2001200002 itself:

• m_T8, rtol 1e-6: on main it fails with CV_ILL_INPUT at t = 91 and finishes only on the GH Four Linux-only test failures, found by #169's first whole-suite run: the GH #176 FD fallback dies at a second crossing, and a steady-state fixture is exactly singular under reference LAPACK #176 finite-difference retry. Here it runs on its own Jacobian, with no retry and 2 rows read where it stood.
• mu, k, k_th: the 28 such rows are now read. All its sensitivities are zero, correctly: none of the three appears in any rate law.

Not this PR's: the sensitivity run's state differs from a plain run by about 3e-3 from t ≈ 41 on, including rows with no bad read. That is the clamp-retry trajectory itself, #1298's domain, which #1300 addresses.

Tests (python/tests/test_sens_rows_outside_last_step.py)

On MPEKRIS, copied from #1300's test:

• the guard: the run reaches the path (n_rows_outside_last_step > 0). sens: evaluate the sensitivity RHS at the value RHS's clamped state (#1298) #1300 stops this model reaching it, so that would make the tests below pass for nothing; this test fails instead and asks for another model;
• no row of sensitivities is left at zero, and all are finite;
• each row is read where the run stands: at or past its grid time, within one interval, in order;
• a run without sensitivities counts none and keeps the grid exactly.

On main: 3 fail, 1 passes (the in-order times).

Mutations, each rebuilt and caught:

• the window check off: the run raises at t = 44 ("CVODES could not read the sensitivities at the output point t=44 … flag -25");
• the per-row read back to CVodeGetSens: the same.

Suites: python/tests and parity_checks/tests give 17,326 passed, 172 skipped, 1 xfailed, every event and sensitivity test included. clang-format 17.0.6 and ruff 0.14.0 are clean. The core stub gains the new SolverStats field.

Not run: a sanitizer build. The issue's memory-safety question is answered by reading SUNDIALS's source: point 3 above. No second cvode_mem is involved.

[wshlavacek]

Thanks for the trace: "a last step from 44.0099" was the thread that led to the cause. I pushed review commits on top of a merge of main (7eab654, 7649533, e9e1903, 8e334e5). They replace the row move, because the root cause is in SUNDIALS, not in the output loop.

The cause

In cvStep's staggered sensitivity corrector, CVODES evaluates f at the converged state. A recoverable failure there goes back round the step loop as nflag = PREV_CONV_FAIL; continue;, with neither cvRestore nor a smaller step. The next attempt predicts again from the advanced tn, so the step covers 2h and its first h is crossed by the uncorrected predictor. Neither the error test nor the root search sees that first h. That is how a CV_NORMAL return ends up with a last step that starts after tout. On MPEKRIS a probe build counted 2,333 of these in one run, starting at t ≈ 39.9. The code is unchanged on SUNDIALS develop; the fix and a unit test are now proposed upstream as llnl/sundials#1073 (and the same pattern in IDAS as llnl/sundials#1074).

What the move left

Reading a row where the run stands fixed the zeros, but the values stayed wrong. Against independent references (a plain run at rtol 1e-12, and central differences of plain runs):

MPEKRIS, rtol 1e-9	this PR as submitted	with the SUNDIALS fix
state, max error (column-scaled)	9.0e-3	3.4e-8
d/dsigma_nk, d/dg_reg	5.2e-3, 1.4e-2	3.1e-6, 2.3e-6
rows off their grid time	25	0
steps	8,283	1,577

On a 10,001-point grid the move put 2,870 rows on repeated times, up to 0.094 past the time each asked for: once the run stood past an output, the next outputs were recorded where it stood. The "3e-3 from t ≈ 41" in the body is this defect, not #1298. The spline check compared each row with neighbours that were just as wrong.

The review commits

• third_party/sundials/carry_cvodes_staggered_restore.cmake is applied as the FetchContent PATCH_COMMAND. It sends that failure through cvHandleNFlag as RHSFUNC_RECVR, like every other recoverable failure in cvStep, and counts it in ncfS. A first revision counted it in ncf, which cvStep resets on every attempt, so a refusal no smaller step cures never reached maxncf. The configure fails if the fetched source lacks the current marker (r2), which also catches a build tree holding the older revision. VENDOR.json → local_carries records the carry and its upstream PR, and vendor_sundials.py now keeps it across a refresh. BNGSIM_SUNDIALS_STAGGERED_RESTORE=0 takes the upstream path, for an A/B only.
• A CV_NORMAL return whose last step does not span tout is refused, not read where the run stands. The check now also runs after retry_while_advancing's CV_NORMAL call, which had none. It mirrors CVodeGetDky's own window: the hu < 0 fuzz was inverted, though no reachable path takes it. The helper moved out of retry_while_advancing's doc comment. A system-SUNDIALS build warns at configure time.
• n_rows_outside_last_step is dropped, since a returned run would always read 0.
• Your per-row CVodeGetSensDky read and its refusal stay. One case needed more: a row recorded before the first step, such as a start time named twice (sample_times=[2, 2, 3]). There CVODES has no step size yet, and its interpolation divides 0 by 0. That row now records the seed row 0 records. Main recorded dA/dA0 = 0 there instead of 1, with no error.
• Tests (test_sens_staggered_corrector_restore.py): the state against a tight plain run, the sensitivities against central differences, rows exactly on the grid (101 and 2,001 points), the hatch to show the model still reaches the defect, and the repeated start time. All 6 fail on the PR as submitted and pass with the fix. The finite-difference check leaves out columns below 1e-6 of the largest, so a column the parameter does not reach (Treg for sigma_nk, whose difference is libm-dependent noise) cannot get in.

Blast radius. A probe build counted the upstream path over all 1,323 rr_parity ODE models (3 parameters each, job tolerances). Three reach it: MODEL2001200002, MODEL2001160002 and MODEL2001140002.

• MODEL2001160002: completes in 1.3 s with the carry. With the carry off, the new guard refuses it at t = 24; main returned its rows silently.
• MODEL2001140002: on main it hangs past its own timeout=15, inside one CVode call. With the carry it fails loudly (the GH CVODES: max_steps does not bound the work at an if(t>=tau) discontinuity — Simulator.run never returns #54 collapsed-step message), and the 15 s budget ends the FD retry.

Suites: python/tests + parity_checks/tests, 17,753 passed, 167 skipped, 0 failed; C++ ctest 6/6 from a fresh fetch, which ran the patch step on a clean download.