Skip to content

Accept a vector abstol with a scalar reltol in in-place FIRK solvers - #4920

Open
singhharsh1708 wants to merge 5 commits into
SciML:masterfrom
singhharsh1708:firk-vector-abstol
Open

singhharsh1708 wants to merge 5 commits into
SciML:masterfrom
singhharsh1708:firk-vector-abstol

Conversation

@singhharsh1708

@singhharsh1708 singhharsh1708 commented Oct 27, 2026 •

Copy link
Copy Markdown
Member

In-place RadauIIA3, RadauIIA5, RadauIIA9, AdaptiveRadau and GaussLegendre now solve with a vector abstol and a scalar reltol instead of throwing a MethodError at init, and in-place RadauIIA3 also takes a vector reltol.

The in-place caches in firk_caches.jl typed atol from reltol, so a scalar reltol gave a scalar field that was then assigned a vector; atol is now an array filled in place. RadauIIA3 passes minimum(reltol) to the linear solver when reltol is a vector (as newton_tolerance does), where master threw a MethodError. AdaptiveRadau reads rtol and atol from the cache after updating them, which keeps empty states working with the array atol.

This changes scalar-tolerance AdaptiveRadau results on purpose: master's first step used the raw reltol for both tolerances, so on this branch Van der Pol (mu = 100, u0 = [2, 0], tspan (0, 2), reltol 1e-6, abstol 1e-7) takes 91 steps instead of 92. The other four methods take the same number of steps as master on that problem, and an equal-valued vector abstol gives the same solve as the scalar one for all five methods. The vector-tolerance and Krylov tests fail without the src change; the empty-state test guards the array atol.

Depends on #4935: the changed AdaptiveRadau step pattern exposes the dense-output bug #4925 in the sparse-vs-dense test at ode_firk_tests.jl:79 (1.33e-4 against an atol of 1e-4), which is why FIRK Core and downgrade are red here. With #4935 stacked the difference is below 1e-6. #4935 also bumps FIRK to 2.8.10, so this PR moves to 2.8.11 once it merges.

Not covered here: RadauIIA5, RadauIIA9, AdaptiveRadau and GaussLegendre never set the linear solver's abstol (only a scalar reltol per solve), so a Krylov solve can stop at LinearSolve's default and return a wrong answer with scalar tolerances too (RadauIIA5 with GMRES gives 0.0 instead of 1e-10; master as well).

Fixes #4911.

AI Disclosure: Claude assisted with this change.

@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): CHANGES, risk medium-low (ML2). Reviewed head: 1fc5c1aa78f6d9c16d0319838603e9abe8bb41c5.

Full review

VERDICT: CHANGES
PR: #4920
Round: 1
Head: 1fc5c1a
Base: 45ff46a
Slot: 1
Reviewer: Codex CLI 0.163.0 / gpt-6.1-sol
Session: 01a19752-68ac-7c02-a250-73139a4a294d
Reviewed: 2026-11-01T13:44:04Z
Summary: This PR enables vector absolute tolerances in five in-place FIRK solvers and vector relative tolerances in RadauIIA3. The regression fails before and passes after, but AdaptiveRadau combines tolerance values from different updates.

Risk

medium-low (ML2, ML8: changes non-default FIRK tolerance initialization and mutable tolerance caches; AdaptiveRadau cache coherence remains defective.)

Risk assessment

  • Risk: medium-low (ML2, ML8: peripheral numerical initialization and cache changes, assessed before reading the PR description.)
  • Blast radius: In-place RadauIIA3, RadauIIA5, RadauIIA9, AdaptiveRadau, and GaussLegendre residual normalization.
  • Evidence: Changed test file passes on head; its copied regression errors on base. Five scalar-tolerance trajectories match base by ==. A separate probe confirms inconsistent AdaptiveRadau residual scaling. FIRK CI is green; four failing jobs match master failure signatures, while one remains unclassified.
  • Verdict: CHANGES
  • Round: 1
  • Blocking priorities: P2
  • Merge: needs changes (make AdaptiveRadau tolerance reads coherent and cover order changes)

Blocking findings

  1. [P2] confirmed by running - AdaptiveRadau uses updated atol with stale scalar rtol. At lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl:2024, the local scalar rtol is captured before the cache update. The array selected at line 2026 aliases cache.atol, which line 2042 mutates. Consequently, residual calculations combine the previous relative tolerance with the current absolute tolerance.

    For u′ = -u, u₀ = [1,2], abstol = [1e-8,1e-6], reltol = 1e-6, and initial dt = 0.01, the probe observed:

    old rtol:       1e-6
    current rtol:   1.0000000000000004e-5
    actual residual:
      [4.1287301097577914e-5, 7.569338534555949e-6]
    residual using current atol and current rtol:
      [4.496636753201553e-6, 3.0277354138223795e-6]
    

    The actual residual exactly matches calculation with the old rtol. On the PR’s Van der Pol problem, equivalent scalar and vector absolute tolerances produce respectively 91 accepted / 13 rejected and 89 accepted / 6 rejected steps.

    Read both tolerances after updating them, or preserve a coherent snapshot of both. Extend lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl:355 to cover AdaptiveRadau through stage-count changes; the equivalence loop currently omits it. The cache-coherence regression required by G3 is missing.

Non-blocking findings

None.

What I ran

Absolute locations are abbreviated below as $JRUN (supplied wrapper), $BASE (supplied merge-base worktree), and $PROBE (temporary adaptive_probe.jl).

  • Git inspection, all exit 0:

    git rev-parse HEAD
    git merge-base origin/master HEAD
    git diff origin/master...HEAD
    git diff --check origin/master...HEAD
    git log -1 --no-merges --format=%B

    Head and merge-base matched the supplied SHAs. No author harness/model trailers were present.

  • Head, exit 0:

    "$JRUN" --project=lib/OrdinaryDiffEqFIRK -e 'using TestEnv; TestEnv.activate(); using Test, SafeTestsets; include("lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl")'

    Tail:

    In-place solve with a vector abstol | 28  28
    
  • Copied the changed test file into the base worktree, then ran the same command there:

    cp lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl "$BASE/lib/OrdinaryDiffEqFIRK/test/ode_firk_tests.jl"

    Base test exit 1, after reaching the new regression:

    MethodError: Cannot convert Vector{Float64} to Float64
    In-place solve with a vector abstol | Error 1
    

    The exception originates in initialize!(::RadauIIA3Cache) at lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl:161.

  • Scalar baseline probe, from base:

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

    Exit 0:

    Julia 1.13.1
    Saved five scalar-tolerance trajectories from base
    
  • Comparison and residual probe, from head:

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

    Exit 1:

    Scalar trajectories unchanged from base | 5 passed
    AdaptiveRadau tolerance consistency     | 2 passed, 2 failed
    

    Compared saved states, times, and accepted/rejected counts for all five solvers on Van der Pol. The two failures are current-tolerance residual consistency and AdaptiveRadau scalar/vector trajectory equivalence.

  • Live PR/issue/check reads:

    gh pr view 4920 --repo SciML/OrdinaryDiffEq.jl --json state,isDraft,headRefOid,baseRefOid
    gh issue view 4911 --repo SciML/OrdinaryDiffEq.jl
    gh pr checks 4920 --repo SciML/OrdinaryDiffEq.jl

    PR remains open at the supplied head. Checks report five failed jobs.

  • Inspected supplied CI excerpts for FIRK Core/QA, downgrade tests, head CI, and head downstream tests. Compared against master CI and master downstream tests: ReverseDiff tape errors, SciMLSensitivity unexpected passes, Enzyme’s mi2 !== nothing assertion, and the GPUCompiler warning assertion match. ModelingToolkit’s reversed-equation assertion has no matching failure in the supplied master excerpts.

Source files remain unchanged in both worktrees; only the copied test file is modified on base.

What I did not verify

  • Package-wide, full downstream, GPU, MPI, or Reactant execution locally.
  • The author’s complete scalar-preservation matrix across out-of-place, fixed-step, ROBER, and multiple tolerances.
  • Remaining base regression cases after the first RadauIIA3 exception.
  • The cause of the unmatched ModelingToolkit failure. Full CI logs were not independently re-downloaded; failure comparison used the supplied excerpts.

AUTHOR_GRADE: {"task":"bugfix","score":3,"difficulty":3,"rounds":1,"author_harness":"unknown","author_model":"unknown","note":"The analytic regression fails on base and passes on head, and five scalar trajectories match base. Substantive repair remains necessary because AdaptiveRadau mixes tolerance generations and is omitted from the equivalence regression; coordinating stage-dependent caches across five solvers makes this moderately difficult."}


🤖 Posted by an AI agent · harness: Claude Code · model: claude-opus-5-5 (coordinating gpt-6.1-sol reviewers)
Review: unreviewed
Session: 01a19752-68ac-7c02-a250-73139a4a294d

@singhharsh1708

Copy link
Copy Markdown
Member Author

The FIRK Core failure here (Sparse FIRK solves vs dense jac_prototype, AdaptiveRadau with the mass matrix) is the dense-output bug from #4925. Both forms interpolate badly at t = 0.01 on master too, and this PR's slightly different step sequence moves their gap over the test's 1e-4:

AdaptiveRadau, M = diag(0.8..1.2), error vs exp(M\J t) at t = 0.01
            sparse    dense     sparse-dense   node error
master      1.36e-2   1.36e-2   2.0e-5         4.2e-9
this PR     1.75e-2   1.76e-2   1.3e-4         4.2e-9
+ #4935     1.79e-7   1.79e-7   1.8e-10        4.2e-9

With #4935 merged on top the testset passes. The two PRs also touch the same line (both read rtol/atol after the update), so whichever lands second needs a trivial rebase.

@ChrisRackauckas

Copy link
Copy Markdown
Member

🤖 Automated comment from an AI agent running as @ChrisRackauckas — not written or reviewed by Chris.

PR: #4920

Reviewed head: 68d1a48
Raw SHA-256: ad5155dee60eddca759e394c9ff6b2db077b1b42f226d76066e4059938ec9e67
Source SHA-256: 9f488a33acb330bf53064f2cc5ea6c840b9a5da015624cc3efb923aad1c1afa0
Source: standard-review-v1 (harness-log)
Reviewer session: 01a1d256-24df-7002-9dd0-6d875daa44f3
Slot: 1
Confidence: unknown
Own risk fixed before author rating/discussion: yes

Risk assessment

  • Risk: high (H1)
  • Verdict: CHANGES
  • Reviewed head: 68d1a48

Review report

VERDICT: CHANGES
Round: 1
Head: 68d1a48
Slot: 1
CI: CI is not green (regress:SciMLSensitivity.jl/Core4/1 / Downstream Tests - Core4; downgrade-sublibraries / test (lib/OrdinaryDiffEqFIRK); sublibrary-ci / lib/OrdinaryDiffEqFIRK / Julia 1 / Tests - Core (Julia 1))
Session: 01a1d256-24df-7002-9dd0-6d875daa44f3

Risk

high (H1: AdaptiveRadau rejects a previously successful public solve input: an empty in-place state).

Blocking findings

  • Verdict: CHANGES
  • Round: 1
  • Blocking priorities: P1, P2
  • Deferred P2 issues: None provided.

[P1] Vector reltol can produce a completely wrong solution with Success — lib/OrdinaryDiffEqFIRK/src/firk_caches.jl:114. Confirmed by running:

using OrdinaryDiffEqFIRK, LinearSolve
f_small!(du, u, p, t) = (du[1] = 1e-10; nothing)
prob = ODEProblem(f_small!, [0.0], (0.0, 1.0))
sol = solve(prob, RadauIIA3(linsolve = KrylovJL_GMRES());
    reltol = [1e-13], abstol = 1e-15)

The head returns Success, reaches t=1.0, and returns [0.0]. The analytic answer is [1e-10]: constant forcing has no Runge–Kutta truncation error. The absolute error is 100,000 times the requested abstol. The same failure occurs with abstol=[1e-15]. Scalar reltol=1e-13 returns approximately [1e-10] on both the merge base and head; LU with vector reltol also returns the correct answer.

The new empty-keyword branch leaves both linear tolerances at 1.4901161193847656e-8, allowing GMRES to accept zero corrections for this small RHS. Derive suitable scalar linear tolerances from component tolerances instead of silently using those defaults, and add this analytic regression test.

[P2] Empty-state AdaptiveRadau solves now throw — lib/OrdinaryDiffEqFIRK/src/firk_perform_step.jl:2052. Confirmed by running:

using OrdinaryDiffEqFIRK
f_empty!(du, u, p, t) = nothing
solve(ODEProblem(f_empty!, Float64[], (0.0, 1.0)), AdaptiveRadau();
    dt = 0.1, adaptive = false)

This succeeds at the merge base. The head throws BoundsError: attempt to access 0-element Vector{Float64} at index [1] from first(cache.atol). The cache has the state's length, so scalar tolerances do not imply a nonempty cache. Handle the empty cache without indexing its first element and add an empty-state regression test.

Non-blocking findings

None.

What I ran

AI review by Codex (version unknown), model gpt-6.1-sol. Conversation: local session 01a1d256-24df-7002-9dd0-6d875daa44f3 (no session URL). No PR body, comments, reviews, or author rating were accessed.

Here $repo denotes the checkout and $scratch the authorized scratch directory. TMPDIR was set to the requested job directory, and the runtime JULIA_DEPOT_PATH was preserved.

  • git diff origin/master...HEAD, git diff --check origin/master...HEAD, and final git status --porcelain=v1: diff inspected, whitespace check passed, working tree clean, head unchanged.
  • julia --version: 1.13.2. Scratch environments used FIRK's declared dependencies, test extras, compat and local sources; their generated manifests fixed the resolution. Relevant versions: LinearSolve 5.18.3, FastBroadcast 1.5.1, SciMLBase 3.58.0. No repository manifest or Julia version pin was present.
  • GROUP=Core JULIA_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 julia --startup-file=no --project="$scratch/env" "$scratch/core.jl" "$repo": the complete FIRK Core runner, 357/357 passed, 6m28s.
  • julia --startup-file=no --project="$scratch/env" "$scratch/probes.jl" "$repo": corrected adversarial probe subtree 170 passed, 1 errored, with only the AdaptiveRadau empty-state error. It covered nonuniform stiff component tolerances, matrix states/tolerances, Float32 and mixed tolerance precision, ForwardDiff gradients, empty states, and six Krylov accuracy cases against a matrix exponential. An initial scratch run had an unqualified helper name; those harness errors were corrected and the affected probes rerun.
  • The same job ran merge-base controls using FIRK from the merge-base commit, with other dependencies held fixed: 10/10 empty-state assertions passed; the PR's added vector-abstol testset failed with the expected Vector{Float64}-to-Float64 conversion error. That added testset passes within the head's Core run.
  • julia --startup-file=no --project="$scratch/env" "$scratch/krylov-probe.jl" "$repo": scaled constant-forcing controls 14 passed, 2 failed; both failures are vector-reltol GMRES cases returning zero with Success. Its merge-base scalar control 2/2 passed. The final baseline driver reused the recorded negative-control result; it did not rerun tests.
  • uv run --locked python --version: Python 3.15.0. No Python project or uv.lock exists; uv warned that --locked has no effect outside a project. Python was not used to run the Julia tests.

What I did not verify

The supplied failing CI jobs were not reproduced in their CI environments, and their failures remain unclassified. Minimum dependency versions, other Julia versions, QA, concurrent execution beyond one Julia thread, GPU, Enzyme and Zygote paths were not verified. No project-locked interpreter was available; Julia verification used the scratch manifests described above.

AUTHOR_GRADE: {"task": "bugfix", "score": 2, "difficulty": 3, "rounds": 1, "author_harness": "unknown", "author_model": "unknown", "note": "Core passes 357/357 and the vector-abstol test fails at the merge base, but probes confirm a successful zero solution for small constant forcing with vector reltol and GMRES, plus a new empty-state BoundsError."}

Head link: 68d1a48


🤖 Posted by an AI agent — harness: codex-cli 0.164.0 · model: gpt-6.1-sol
Conversation: local session 01a1d256-24df-7002-9dd0-6d875daa44f3 (no session URL)

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.

In-place FIRK methods throw a MethodError with a vector abstol and a scalar reltol

2 participants