Fix prob_ode_filament: mistranslated 0.6 deprecations in the projection - #219
Merged
ChrisRackauckas merged 1 commit intoAug 8, 2026
Merged
Conversation
`projection!` computed the constraint projector with `mul!(P', P0, J)`. The 0.6-era call was `At_mul_B!(P, P0, J)` (P <- P0' * J); the transpose belongs on the operand, not the destination, so this wrote a 20x63 * 20x63 product into a 63x63 destination and threw `DimensionMismatch` on every RHS evaluation. `subtract_from_identity!` had the same class of error: `scale!(-1, A)` was translated to `rmul!(-1, A)`, but `rmul!` takes (array, number). The correct replacement is `lmul!(-1, A)`. That call is unreachable until the `mul!` above is fixed, so both are needed. The fixed RHS reproduces the SciMLBenchmarks Filament.jmd model bit-for-bit. Add a Core regression test asserting the RHS evaluates and that the projection holds every squared segment length constant, which is what `P` is for. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
ChrisRackauckas
marked this pull request as ready for review
August 8, 2026 11:24
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Fixes #218.
What was wrong
Two mistranslated Julia 0.6 → 1.0 deprecations in
filament_prob.jl, both in theinextensibility projection.
prob_ode_filamenthas been unusable since 2019 —it threw on a single RHS evaluation, with no solver involved.
projection!(introduced in 2913a1e, Jun 2019). The 0.6 call wasAt_mul_B!(P, P0, J), i.e.P <- P0' * J. It becamemul!(P', P0, J)— thetranspose was moved from the operand to the destination.
P0isN x 3(N+1)andJisN x 3(N+1), so this asks for a20x63 * 20x63product written into a
63x63destination. Correct translation:mul!(P, P0', J).subtract_from_identity!(introduced in e56b7c2, Jul 2018).scale!(-1, A)became
rmul!(-1, A), butrmul!takes(array, number)— the number-firstform is
lmul!. This throwsMethodError. It is unreachable until (1) isfixed, so both changes are required.
lmul!is added to theusing LinearAlgebra:list.The bug survived seven years because the test suite only checked that the module
loads; no test ever evaluated a problem's RHS. This PR adds one for the filament.
Failing before
GROUP=Coreon this branch with onlysrc/reverted to master:With the
mul!fix applied butlmul!reverted tormul!, showing the secondchange is independently needed:
Passing after
Correctness
Not just non-erroring. The fixed RHS reproduces the working copy of the same
model in
SciMLBenchmarks
benchmarks/ComplicatedPDE/Filament.jmdbit for bit — the reference model was extracted from the
.jmdinto aseparate module and evaluated side by side:
(t = 0 is near zero because the initial configuration is straight along x and
the t = 0 magnetic force is purely axial, so the projection removes all of it.)
Whole solutions agree with the reference to 0.0 as well, at
abstol = reltol = 1e-10:The physical invariant the projection exists to enforce — constant segment
lengths — holds to integration-tolerance level:
Tolerance sweep on the default problem (
tspan = (0.0, 1.0)), after warmup:Two pre-existing limitations, unchanged by this PR
Both are inherited from the upstream benchmark model and are not touched here;
flagging them because they affect how the problem should be used.
The analytic
jacis an approximation. It isP * A, which drops the∂P/∂r · (Ar + F)term —Pdepends onr. Against central differences of thetrue RHS at
u0:It does not give wrong answers, but it does cost Rosenbrock methods heavily,
since an inexact Jacobian breaks their order conditions (
tspan=(0,0.01),1e-10):18.8x more Rodas5P steps with the analytic Jacobian; FBDF is indifferent, as
expected for a Newton-based method. This is very likely the reason Filament.jmd
drops the Rosenbrock methods as "500x slower". Supplying the exact Jacobian means
differentiating the projector, which changes the model and every published
benchmark number for it — a separate decision, not this PR.
Rodas5P()with default autodiff fails, because the model's scratch buffersare
Vector{Float64}and Rosenbrock methods ForwardDiff throughtfor the timegradient:
Rodas5P(autodiff = AutoFiniteDiff())works, which is what Filament.jmd usesthroughout.
FBDF()works with defaults (it uses the suppliedjacand needs notgrad). Fixing this properly means either an analytictgrador AD-compatiblecaches; again separate from this bug.
Not verified
The OrdinaryDiffEq.jl downstream jobs and the Downgrade jobs were not run
locally — left to CI. Runic 1.7.0 and
typosare clean on the diff. No versionbump; this repo does those in dedicated release-prep PRs.
Judgment calls a reviewer may want to push back on: the regression test's
inextensibility threshold is
rate / max(1, ‖du‖∞) < 1e-10where the measuredvalue is ~2.7e-14, i.e. ~3700x margin — tight enough to catch a wrong projector
(which gives O(1)) but loose enough to survive BLAS differences. And the two
limitations above were deliberately left alone rather than fixed here.
CI note
Every check passes except
Downgrade / Downgrade Tests - Core, which is notcaused by this PR — it reproduces on unmodified master at this PR's base commit
3a884a6. It fails inside Aqua's stale-dependency check on the rootDiffEqProblemLibrarypackage:Root cause: the
julia-actions/julia-downgrade-compat@v2moving tag advanced fromv2.6.2 to v2.7.0
on 2026-07-31, between this repo's last green Downgrade run (2026-07-29, same base
commit) and today. v2.7.0 rewrites the checkout's
Project.tomlto promoteold-style
[extras]/[targets].testnames into[deps]:Aqua then reads that rewritten
Project.tomland correctly reports two[deps]the package never loads. The SciML reusable workflow SHA and the resolved
dependency versions are identical between the green and red runs.
This is a known upstream regression affecting other Julia packages too —
julia-actions/julia-downgrade-compat#58, filed
2026-08-04 by a third party. It needs fixing in the action (or worked around
here), not in this PR.
This PR touches only
lib/ODEProblemLibrary/. The non-downgradedtests / QA (julia 1)job passes, anddowngrade-sublibraries / test (lib/ODEProblemLibrary)— the downgrade job for the sublibrary actually changedhere — passes.