Fixes #3915.

Bug

jacobian2W!/jacobian2W special-case mass_matrix isa UniformScaling to skip the axes(mass_matrix) == axes(W) boundscheck and use the scalar directly. ScalarOperator is the same thing (λ·I) but isn't a UniformScaling, and reports axes(mm) == (), so it falls through to the boundscheck meant for full matrices and throws.

using OrdinaryDiffEqSDIRK, SciMLOperators
f = ODEFunction((du, u, p, t) -> (du .= -u); mass_matrix = ScalarOperator(1.0))
solve(ODEProblem(f, [1.0], (0.0, 1.0)), TRBDF2())
# DimensionMismatch: W: (Base.OneTo(1), Base.OneTo(1)), mass matrix: ()

IdentityOperator doesn't hit this since its axes happen to match W.

Fix

Treat ScalarOperator the same as UniformScaling in both jacobian2W! methods and jacobian2W, via _is_scalar_massmatrix/_scalar_massmatrix_λ. Follows #3814 → SciMLOperators#400 → #3837, the same operator-mass-matrix chain, one step further (this is the first place a ScalarOperator mass matrix reaches once the entry crash and get_differential_vars are both already fixed).

Tests

New scalar_operator_massmatrix_tests.jl covers jacobian2W! through both the ::Matrix specialization and the generic ::AbstractMatrix method (reached with a sparse J), the out-of-place jacobian2W, and an end-to-end TRBDF2() solve of the case from #3915 — Lorenz with mass_matrix = ScalarOperator(2.0) against mass_matrix = 2.0*I. 7 assertions; the file errors on master with the reported DimensionMismatch.

Those two solves are in fact bit-identical (max abs difference 0.0 over saveat = 0:0.1:1 at abstol = reltol = 1e-12, identical accepted-step counts) for TRBDF2, KenCarp4, FBDF and Rodas5P, and Rodas5P agrees to 7e-12 with a 1e-14 reference solve of the equivalent M = I, f/2 problem.

OrdinaryDiffEqDifferentiation Core suite: 65 passes across 10 testsets (master: 61 across 9). QA also passes (JET 1/1, Aqua/ExplicitImports 15/15).

Scope

This fixes the in-place W assembly, which is what #3915 reports. Still failing and unchanged by this PR, with different errors: out-of-place and StaticArrays ScalarOperator mass matrices (Out-of-place NLNewton does not support non-concrete W, and an AssertionError in the StaticWOperator path), and a sparse jac_prototype combined with a ScalarOperator mass matrix, which fails earlier in prepare_user_sparsity.

AI Disclosure

Claude assisted with this work.