Repository navigation
Restart multistep history at d_discontinuities - #4932
singhharsh1708 wants to merge 7 commits into
Conversation
Independent review (Codex CLI 0.163.0, model gpt-6.1-sol, reasoning effort high): CHANGES, risk medium-high (MH1). Reviewed head: Full reviewVERDICT: CHANGES Riskmedium-high (MH1: Changes shared discontinuity handling used by the FBDF main lane and multiple other solver families; also L2 and L5.) Risk assessment
Blocking findings
Non-blocking findingsNone. What I ranAll Julia execution used the prescribed wrapper: Julia 1.13.1, four Julia threads, one OpenBLAS thread, isolated temporary directories, and the global lock. Commands below abbreviate that wrapper as
No source files were changed. Generated manifests were removed; only the authorized copied tests remain modified on base. What I did not verify
AUTHOR_GRADE: {"task":"bugfix","score":2,"difficulty":3,"rounds":1,"author_harness":"unknown","author_model":"unknown","note":"The intended regression tests fail on base and pass on head, but the shared flag introduces severe silent Radau accuracy regressions and lacks independent accuracy assertions. The task requires coordinating restart and cache behavior across several solver families."} 🤖 Posted by an AI agent · harness: Claude Code · model: claude-opus-5-5 (coordinating gpt-6.1-sol reviewers) |
6239172 to
ef16b2e
Compare
ef16b2e to
04db522
Compare
PR: #4932 Reviewed head: 04db522 Risk assessment
Review reportVERDICT: CHANGES Riskmedium-high (MH1: shared discontinuity handling changes the FBDF main lane and history/Jacobian reset behavior across solver families). Blocking findings
Both findings reproduce with: using OrdinaryDiffEqAdamsBashforthMoulton, OrdinaryDiffEqBDF
f!(du, u, p, t) = (du[1] = -3 * u[1] + (t > 0.5); nothing)
prob = ODEProblem(f!, [0.0], (0.0, 1.0))
exact(t) = -expm1(-3 * max(t - 0.5, 0.0)) / 3
for alg in (VCAB3(), VCAB4(), VCAB5(), QNDF2())
sol = solve(prob, alg; d_discontinuities = [0.5],
reltol = 1e-9, abstol = 1e-11)
@show typeof(alg) sol.retcode
@show maximum(abs(u[1] - exact(t)) for (t, u) in zip(sol.t, sol.u))
endObserved maximum saved-point absolute errors for this problem at
The control redefines only
Non-blocking findings
What I ran
What I did not verify
AUTHOR_GRADE: {"task":"bugfix","score":2,"difficulty":3,"rounds":1,"author_harness":"unknown","author_model":"unknown","note":"The intended equivalence tests fail before and pass after, and both native Core suites pass, but independent analytic probes confirm silent accuracy regressions in VCAB3/4/5 and QNDF2."} Head link: 04db522 🤖 Posted by an AI agent — harness: codex-cli 0.164.0 · model: gpt-6.1-sol |
When the integrator lands on a
d_discontinuitiesentry,update_fsal!now setsintegrator.derivative_discontinuity, so the step after it restarts the multistep history from the post-discontinuity state, as an emptyDiscreteCallbackat the same time already does. The flag is set insideapply_step!, afterloopheader!checks it, andloopfooter_reset!clears it before callbacks run, so only the first attempt of that step'sperform_step!sees it; a rejected attempt is retried without the flag, which is why the VCAB and QNDF2 startup paths must handle rejection themselves.Turning the flag on exposed two restart bugs, which this PR also fixes. VCAB3/4/5 and VCABM3/4/5 advanced their step counter and step-size history on a rejected startup step, so after a restart the retry took an Adams step from that history with a zero error estimate; this already broke VCAB and VCABM on the callback path on master (0.019 for VCAB3 and 5.8e-6 for VCABM3 with an empty callback). A rejected startup step now leaves the history alone, as a rejected full-order step already does. The same fix covers a rejected first step of a solve: on u' = -3u with dt = 0.4 and reltol 1e-8, VCAB3 on master has an error of 1.1e-3 and now 1.2e-9, and VCABM5 goes from 4.4e-7 to 5.0e-8.
QNDF2 accepted its first step after a restart with EEst set to 1 while using the step size from before the restart. That step now gets a backward Euler error estimate, (u - uprev - dt f(uprev)) / 2. With a mass matrix the residual M (u - uprev) - dt f(uprev) is mapped back to the units of u through the Newton matrix, as (M - γdt J) \ r, so scaling M and f by 1e-9 gives the same 9.6e-8 error as M = I. That solve reuses the nonlinear solver's own factorization of W and is taken only when
can_smooth_est(nlsolver)is true, the check SDIRK uses for its smoothed estimate: NLNewton, or NonlinearSolveAlg in place with W reuse. Without a reusable W (out-of-place NonlinearSolveAlg, NLFunctional, NLAnderson), QNDF2 with a mass matrix keeps its history at a discontinuity, as master does withd_discontinuities, and now also after a callback, where master restarted and accepted the first step with EEst = 1. On u' = -3u + (t > 0.5) with a 3x3 identityMatrixas the mass matrix and out-of-place NonlinearSolveAlg, this gives 6.0e-8 with eitherd_discontinuitiesor an empty callback, where master gives 6.0e-8 and 0.059. The first step of a solve keeps its old behaviour.On u' = -3u + (t > 0.5) at reltol 1e-9, setting the flag without the solver fixes gives errors of 5.9e-4, 3.1e-5, 1.6e-5 and 0.059 for VCAB3, VCAB4, VCAB5 and QNDF2; with them the errors are 4.8e-11, 3.4e-11, 2.3e-11 and 9.6e-8. Master, which does not restart there, gives 5.6e-11, 3.4e-11, 8.6e-11 and 8.7e-8. The solver fixes leave QNDF, FBDF, QNDF1 and ABDF2 unchanged. Other methods that read the flag now see it at a
d_discontinuitiesentry too, as they already do after a callback. In-place RadauIIA5, RadauIIA9 and AdaptiveRadau lose accuracy there through a W-reuse bug that a callback triggers on master as well; #4941 fixes it and should land first.One case does not improve. On an index-1 DAE whose algebraic variable jumps at the discontinuity, QNDF2 returns Unstable on master with
d_discontinuitiesand now also with an empty callback, where master with the callback returned Success with a 0.059 error. With the callback the restart step is accepted, but the jump stays in QNDF2's difference history and the next step fails its error test. Withd_discontinuitiesthe algebraic variable is not reinitialized, so the restart step's error estimate stays near the jump size and never passes. QNDF and FBDF solve it to 1e-8 or better. Keeping the history after a callback costs one case that master solved: on u1' = -u1 + u2, 0 = u2 - (t > 0.5) with M = [1 0; 0 0] and an empty callback at 0.5, out-of-place NonlinearSolveAlg now returns Unstable at t = 0.5, where master returns Success with a 4.4e-6 error. NLNewton on the same problem goes from 4.4e-6 to 8.0e-7.The new tests check that
d_discontinuitiesalone gives the sametanduas the empty callback, plus a Success return code, the final time, and the error against the exact solution, for VCAB3/4/5, VCABM3/4/5, FBDF, QNDF and QNDF2, for VCAB3/4/5 and VCABM3/4/5 through the callback path, and for QNDF2 with a mass matrix of 1e-9, in place and out of place, with NLNewton and NonlinearSolveAlg; 14 of them fail without the solver fixes.Fixes #4913
AI Disclosure: Claude assisted with this change.