Skip to content

Fix prob_ode_filament: mistranslated 0.6 deprecations in the projection - #219

Merged
ChrisRackauckas merged 1 commit into
SciML:masterfrom
ChrisRackauckas-Claude:fix-filament-projection
Aug 8, 2026
Merged

Fix prob_ode_filament: mistranslated 0.6 deprecations in the projection#219
ChrisRackauckas merged 1 commit into
SciML:masterfrom
ChrisRackauckas-Claude:fix-filament-projection

Conversation

@ChrisRackauckas-Claude

@ChrisRackauckas-Claude ChrisRackauckas-Claude commented Aug 8, 2026

Copy link
Copy Markdown
Member

Please ignore this PR until it has been reviewed by @ChrisRackauckas.

Fixes #218.

What was wrong

Two mistranslated Julia 0.6 → 1.0 deprecations in filament_prob.jl, both in the
inextensibility projection. prob_ode_filament has been unusable since 2019 —
it threw on a single RHS evaluation, with no solver involved.

  1. projection! (introduced in 2913a1e, Jun 2019). The 0.6 call was
    At_mul_B!(P, P0, J), i.e. P <- P0' * J. It became mul!(P', P0, J) — the
    transpose was moved from the operand to the destination. P0 is
    N x 3(N+1) and J is N x 3(N+1), so this asks for a 20x63 * 20x63
    product written into a 63x63 destination. Correct translation: mul!(P, P0', J).

  2. subtract_from_identity! (introduced in e56b7c2, Jul 2018). scale!(-1, A)
    became rmul!(-1, A), but rmul! takes (array, number) — the number-first
    form is lmul!. This throws MethodError. It is unreachable until (1) is
    fixed, so both changes are required.

lmul! is added to the using 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=Core on this branch with only src/ reverted to master:

prob_ode_filament: Error During Test at lib/ODEProblemLibrary/test/filament.jl:16
  Got exception outside of a @test
  DimensionMismatch: matrix A has axes (Base.OneTo(20),Base.OneTo(63)), matrix C has axes (Base.OneTo(63),Base.OneTo(63))
   [6] projection!(f::ODEProblemLibrary.FilamentCache{...})
   [7] (::ODEProblemLibrary.FilamentCache{...})(dr::Vector{Float64}, r::Vector{Float64}, p::SciMLBase.NullParameters, t::Float64)

Test Summary:       | Pass  Error  Total  Time
Filament            |    1      1      2  5.8s
ERROR: LoadError: Some tests did not pass: 1 passed, 0 failed, 1 errored, 0 broken.

With the mul! fix applied but lmul! reverted to rmul!, showing the second
change is independently needed:

  MethodError: no method matching rmul!(::Int64, ::Matrix{Float64})
Test Summary:       | Pass  Error  Total  Time
Filament            |    1      1      2  8.6s

Passing after

$ GROUP=Core julia --project=lib/ODEProblemLibrary -e 'using Pkg; Pkg.test()'
Test Summary: | Pass  Total  Time
Load Tests    |    1      1  0.0s
Test Summary: | Pass  Total  Time
Filament      |    8      8  4.5s
     Testing ODEProblemLibrary tests passed

$ GROUP=QA julia --project=lib/ODEProblemLibrary -e 'using Pkg; Pkg.test()'
Test Summary: | Pass  Total   Time
Aqua          |   20     20  38.8s
     Testing ODEProblemLibrary tests passed

Correctness

Not just non-erroring. The fixed RHS reproduces the working copy of the same
model in
SciMLBenchmarks benchmarks/ComplicatedPDE/Filament.jmd
bit for bit — the reference model was extracted from the .jmd into a
separate module and evaluated side by side:

#### RHS: library vs SciMLBenchmarks reference
  t=0.000  ‖du‖∞=1.386624e-11  max|Δ|=0.000e+00  rel=0.000e+00
  t=0.013  ‖du‖∞=2.075123e+04  max|Δ|=0.000e+00  rel=0.000e+00
  t=0.026  ‖du‖∞=1.403473e+04  max|Δ|=0.000e+00  rel=0.000e+00
  t=0.039  ‖du‖∞=2.244388e+04  max|Δ|=0.000e+00  rel=0.000e+00
  t=0.052  ‖du‖∞=2.108862e+04  max|Δ|=0.000e+00  rel=0.000e+00
  t=0.065  ‖du‖∞=1.519610e+04  max|Δ|=0.000e+00  rel=0.000e+00

(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:

tspan(0,0.01), Rodas5P(autodiff=AutoFiniteDiff())  lib=Success ref=Success  max|Δu(T)|=0.000e+00
tspan(0,1.0),  FBDF()                              lib=Success ref=Success  max|Δu(T)|=0.000e+00

The physical invariant the projection exists to enforce — constant segment
lengths — holds to integration-tolerance level:

initial: min=0.050000000000 max=0.050000000000 (1/N=0.050000000000)
Rodas5P t≤0.01   max rel drift over solution=2.774e-08   final min=0.049999998613 max=0.050000000036
FBDF t≤1.0       max rel drift over solution=1.937e-07   final min=0.049999990362 max=0.049999994270

Tolerance sweep on the default problem (tspan = (0.0, 1.0)), after warmup:

  tol=1e-4  FBDF()       retcode=Success  steps=1321   3.73s
  tol=1e-4  CVODE_BDF()  retcode=Success  steps=8655   8.27s
  tol=1e-6  FBDF()       retcode=Success  steps=1950   3.39s
  tol=1e-6  CVODE_BDF()  retcode=Success  steps=4372   4.89s
  tol=1e-8  FBDF()       retcode=Success  steps=4008   12.26s
  tol=1e-8  CVODE_BDF()  retcode=Success  steps=15288  11.24s
  tol=1e-10 FBDF()       retcode=Success  steps=9437   19.77s
  tol=1e-10 CVODE_BDF()  retcode=Success  steps=48129  24.38s
  tol=1e-12 FBDF()       retcode=Success  steps=29844  37.72s
  tol=1e-12 CVODE_BDF()  retcode=Success  steps=64531  39.66s

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 jac is an approximation. It is P * A, which drops the
∂P/∂r · (Ar + F) term — P depends on r. Against central differences of the
true RHS at u0:

rel Frobenius err = 2.0456e-02 ; rel max-entry err = 2.2517e-02

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):

Rodas5P  analytic jac  Success  steps=19670  68.5s  u(T)[1:2]=[0.4279443915, -0.4927683452]
Rodas5P  FD jac        Success  steps=1048   31.6s  u(T)[1:2]=[0.4279443928, -0.4927683469]
FBDF     analytic jac  Success  steps=230    3.9s   u(T)[1:2]=[0.4279443933, -0.4927683465]
FBDF     FD jac        Success  steps=228    5.4s   u(T)[1:2]=[0.4279443933, -0.4927683463]

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 buffers
are Vector{Float64} and Rosenbrock methods ForwardDiff through t for the time
gradient:

Rodas5P()  FAILED: First call to automatic differentiation for time gradient

Rodas5P(autodiff = AutoFiniteDiff()) works, which is what Filament.jmd uses
throughout. FBDF() works with defaults (it uses the supplied jac and needs no
tgrad). Fixing this properly means either an analytic tgrad or AD-compatible
caches; 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 typos are clean on the diff. No version
bump; 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-10 where the measured
value 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 not
caused 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 root
DiffEqProblemLibrary package:

Stale dependencies: Test Failed at Aqua/src/stale_deps.jl:32
  Expression: isempty(stale_deps)
   Evaluated: isempty(Base.PkgId[SafeTestsets [1bc83da4-...], SciMLTesting [09d9d899-...]])

Root cause: the julia-actions/julia-downgrade-compat@v2 moving tag advanced from
v2.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.toml to promote
old-style [extras]/[targets].test names into [deps]:

[ Info: Promoted old-style test dependencies into [deps] for locked Pkg.test
│   promoted = ["SafeTestsets", "SciMLTesting"]

Aqua then reads that rewritten Project.toml and 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-downgraded
tests / QA (julia 1) job passes, and downgrade-sublibraries / test (lib/ODEProblemLibrary) — the downgrade job for the sublibrary actually changed
here — passes.

`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
ChrisRackauckas marked this pull request as ready for review August 8, 2026 11:24
@ChrisRackauckas
ChrisRackauckas merged commit a5551f0 into SciML:master Aug 8, 2026
26 of 27 checks passed
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.

prob_ode_filament throws DimensionMismatch on a single RHS evaluation

2 participants