Skip to content

Reuse in-place FIRK W only at the dt it was factored for - #4941

Open
singhharsh1708 wants to merge 1 commit into
SciML:masterfrom
singhharsh1708:firk-w-reuse-dt
Open

singhharsh1708 wants to merge 1 commit into
SciML:masterfrom
singhharsh1708:firk-w-reuse-dt

Conversation

@singhharsh1708

Copy link
Copy Markdown
Member

In-place RadauIIA3, RadauIIA5, RadauIIA9 and AdaptiveRadau now reuse a factored W only at the dt it was built for. They reused it for a step up to 20% away (new_W_γdt_cutoff = 1//5) while still accepting the first Newton iterate on the previous step's convergence rate, so a step with a stale W could be accepted after one iteration with its error unresolved. cache.W_γdt was only updated with a new Jacobian, so in practice reuse stayed off until something forced one: any callback that fires and leaves u_modified set (it marks a derivative discontinuity), or a d_discontinuities entry once #4932 lands.

W_γdt is now recorded whenever W is rebuilt, and the default new_W_γdt_cutoff for these four methods is 0, so W is kept only while the controller holds dt fixed, which is how RADAU5 reuses its decomposition. Out of place, W is built every step and nothing changes.

On u' = -3u + 1 with a no-op DiscreteCallback at t = 0.5 (reltol 1e-8), master's in-place RadauIIA5, RadauIIA9 and AdaptiveRadau end 68, 779 and 736 times reltol off at their own steps; here they match the run without the callback (0.031, 0.0008 and 0.0026). With #4932 applied, a drive switching at 0.5:0.5:2.5 given as d_discontinuities goes from up to 1.6e4 times reltol to below 0.04. On Van der Pol (mu = 1e3), OREGO and ROBER at 1e-6 and 1e-9, RadauIIA5 needs fewer factorizations for the same steps (Van der Pol 1e-9: 2030 to 720; OREGO 1e-9: 1227 to 814). Final errors are mostly smaller (Van der Pol RadauIIA5 1e-9: 8.3 to 0.23 times tolerance; OREGO 1e-6: 13 to 0.7) but not everywhere (ROBER RadauIIA9 1e-9: 8 to 19), and AdaptiveRadau on Van der Pol at 1e-9 takes 627 steps against 302.

The new testset is the callback problem; it fails on master for all three methods. FIRK Core and QA pass, as do the NonlinearSolve mass-matrix tests (which the Core group never reaches on master, see #4681) and the root Integrators_I, InterfaceI, InterfaceIII, InterfaceIV and Regression_II groups. RadauIIA3 is left out of the test: after a discontinuity its error is 2.6e-5 at every tolerance, on master too and out of place as well, which is a separate problem.

FIRK goes to 2.8.10, the same version as the other open FIRK PRs.

AI Disclosure: Claude assisted with this change.

@oscardssmith

Copy link
Copy Markdown
Member

This is going to substantially slow down the method.

@singhharsh1708

Copy link
Copy Markdown
Member Author

@oscardssmith Measured against master, in-place, reltol 1e-9, abstol 1e-12. Master refactors W on almost every step, because W_γdt is only set when J is rebuilt, so the 1/5 reuse rarely kicks in.

                        steps   factorizations (master -> PR)
VdP mu=1e3  RadauIIA5     982   938 -> 406
OREGO       RadauIIA5     737   585 -> 342
ROBER       RadauIIA5     283   263 -> 266

@singhharsh1708

Copy link
Copy Markdown
Member Author

🤖 Automated comment from an AI agent running as @singhharsh1708, not written or reviewed by the account holder.

Independent review (Codex CLI 0.163.0, model gpt-6.1-sol, reasoning effort high): MERGE-WITH-ISSUES, risk medium (M2). Reviewed head: d8c42744e3ee30590470abad988d425cc8dc75bd.

Full review

VERDICT: MERGE-WITH-ISSUES
PR: #4941
Round: 1
Head: d8c4274
Base: 45ff46a
Slot: 1
Reviewer: Codex CLI 0.163.0 / gpt-6.1-sol
Session: 01a19788-e507-7853-8c58-fd1f0bfe7eb2
Reviewed: 2026-11-01T15:03:01Z
Summary: The PR records the timestep associated with rebuilt FIRK stage matrices and changes Radau defaults to prevent inaccurate cross-timestep reuse. The analytic callback regression fails on base and passes on head; two non-blocking improvements remain.

Risk

medium (M2: Factorization reuse changes numerical results on existing in-place Radau paths, beyond rounding differences.)

Risk assessment

  • Risk: medium (M2: Changes numerical behavior on working peripheral solver paths; independent analytic evidence verifies the callback fix.)
  • Blast radius: In-place RadauIIA3, RadauIIA5, RadauIIA9 and AdaptiveRadau stage-matrix reuse. Default solver selection and other solver families are unchanged.
  • Evidence: The changed test file passes on head and fails all three added assertions on base. An additional analytic order-change probe passes 12/12 checks on both revisions. FIRK Core and QA CI are green.
  • Verdict: MERGE-WITH-ISSUES
  • Round: 1
  • Blocking priorities: none
  • Merge: merge

Blocking findings

None.

Non-blocking findings

  1. [P3] Invalidate stage matrices when AdaptiveRadau changes order - confirmed by running. At lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl:2058, recording the latest factorization timestep enables reuse when dt remains fixed. However, lib/OrdinaryDiffEqFIRK/src/controllers.jl:41 changes the stage count without invalidating those matrices, whose shifts also depend on the tableau. With dtmax = 0.02, both manufactured-solution cases required 153 steps and one nonlinear-convergence failure on head, versus 151 steps and zero failures on base. Both recovered and retained accuracy. Include the order/tableau in cache validity, or force rebuilding after an order change, and test that transition.

  2. [P3] Assert successful completion in the regression - suspected from reading. lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl:352 checks accuracy only over returned points. An accurate partial solution terminating before the callback could satisfy that assertion. Add checks for ReturnCode.Success and sol.t[end] == 1.0.

What I ran

Location aliases below represent the supplied wrapper ($JRUN), merge-base worktree ($BASE_WORKTREE), and scratch diagnostic script ($PROBE); absolute locations are omitted as requested.

Read-only git checks, all exit 0:

git rev-parse HEAD
git diff origin/master...HEAD
git merge-base origin/master HEAD
git log -1 --no-merges --format=%B
git diff --check origin/master...HEAD
git status --short
git -C "$BASE_WORKTREE" rev-parse HEAD
git -C "$BASE_WORKTREE" diff --name-only

The head and merge-base SHAs matched the supplied revisions. The head remained clean; base differed only by the authorized test-file copy. No agent trailers were present.

Copied the changed test file into base and verified identical contents; both commands exited 0:

cp lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl "$BASE_WORKTREE/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl"
cmp lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl "$BASE_WORKTREE/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl"

Ran this exact Julia command separately from head and base:

"$JRUN" --project=lib/OrdinaryDiffEqFIRK -e 'using TestEnv; TestEnv.activate(); using Test, SafeTestsets; include("lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl")'
  • Head: exit 0. Tail: In-place FIRK stays accurate after a callback | 3 3.
  • Base: exit 1. Tail: 0 passed, 3 failed, 0 errored, 0 broken. Maximum errors were 6.789798343009856e-7, 7.787375100709593e-6, and 7.356080193643422e-6, each exceeding 1e-8.

Ran the diagnostic script separately from head and base:

"$JRUN" --project=lib/OrdinaryDiffEqFIRK "$PROBE"

The script used the analytic solution (u(t)=\sin(\omega t)) for
(u'=\omega\cos(\omega t)-\lambda(u-\sin(\omega t))), with (λ, ω) equal to (3, 1) and (1000, 3). It compared fixed order 5 against adaptive orders 5–9 on [0, 3], using dt=1e-4, dtmax=0.02, reltol=1e-6, and abstol=1e-9. Assertions checked completion, analytic accuracy, and accepted-step matrix consistency.

  • Head: exit 0. Tail: Order-sensitive W reuse | 12 12. Adaptive-order cases: 153 steps, one nonlinear-convergence failure each.
  • Base: exit 0. Tail: Order-sensitive W reuse | 12 12. Adaptive-order cases: 151 steps, zero nonlinear-convergence failures.
  • Both revisions produced the same maximum analytic errors in these cases, at most 5.59e-9.

All Julia execution used the supplied wrapper: Julia 1.13.1, four threads, one BLAS thread, and serialized execution.

Inspected the supplied CI snapshots and failure excerpts:

  • FIRK Core and QA passed.
  • Compared failed AD/GPU jobs and downstream jobs against recent master CI and integration logs. The ReverseDiff tape errors, unexpected adjoint passes, Enzyme threading assertion, and GPUCompiler warning failure have matching signatures on master.

What I did not verify

  • The ModelingToolkit equation-orientation failure has no matching signature in the supplied master excerpts; I did not reproduce or classify it as pre-existing.
  • I did not rerun package-wide, downstream, dedicated GPU/AD, or cluster suites.
  • I did not reproduce the author’s broader performance claims or measure runtime/work-precision changes.

AUTHOR_GRADE: {"task":"bugfix","score":4,"difficulty":3,"rounds":1,"author_harness":"unknown","author_model":"unknown","note":"The analytic callback regression fails on base and passes on head, and targeted numerical probes retain accuracy. Minor follow-ups are order-change cache invalidation to avoid a recoverable failed attempt and explicit completion assertions; the task requires moderate solver-cache reasoning."}


🤖 Posted by an AI agent · harness: Claude Code · model: claude-opus-5-5 (coordinating gpt-6.1-sol reviewers)
Review: unreviewed
Session: 01a19788-e507-7853-8c58-fd1f0bfe7eb2

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants