do_newJW: refactorize W on the linear path when γdt drifts (#3934)

For islinearfunction problems do_newJW returned (false, false) after the
first step, assuming a constant J means W never changes — but
W = J - M/(γdt) depends on the step size, so concrete-A linear solvers
kept the factorization from the first step's dt forever. Combined with
WOperator._concrete_form being stale for operator Jacobians
(SciML/SciMLOperators.jl#408), implicit solvers on
SplitODEProblem(MatrixOperator(A), f) with LU/KLU returned garbage with
retcode Success in ~2 steps (#3933).

Keep new_jac = false on the linear path (J really is constant) but apply
the standard W_γdt cutoff test for new_W, matching the DAE branch below.

Adds an end-to-end regression test (KenCarp3/4/5 + LU on a 2x2 split
operator problem, non-split TRBDF2, and the always-working Krylov path);
requires SciMLOperators >= 1.24.4 (compat bumped) for the companion fix,
so CI will resolve only once that releases.

Fixes #3933 together with SciML/SciMLOperators.jl#408.

Co-authored-by: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>

Author: ChrisRackauckas-Claude

lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl

@@ -509,7 +509,17 @@ function do_newJW(integrator, alg, nlsolver, repeat_step)::NTuple{2, Bool}
     integrator.iter <= 1 && return true, true # at least one JW eval at the start
     repeat_step && return false, false
     islin, _ = islinearfunction(integrator)
-    islin && return false, false # no further JW eval when it's linear
+    if islin
+        # J is constant for a linear function, but W = J - M/(γdt) still depends
+        # on γdt: W must be rebuilt/refactorized when the step size has drifted
+        # past the cutoff, otherwise concrete-A linear solvers keep a stale
+        # factorization from the first step's (possibly tiny) dt.
+        isnewton(nlsolver) || return false, true
+        W_iγdt = inv(nlsolver.cache.W_γdt)
+        iγdt = inv(nlsolver.γ * integrator.dt)
+        smallstepchange = abs(iγdt / W_iγdt - 1) <= get_new_W_γdt_cutoff(nlsolver)
+        return false, !smallstepchange
+    end
     !integrator.opts.adaptive && return true, true # Not adaptive will always refactorize
     errorfail = OrdinaryDiffEqCore.get_EEst(integrator) > one(OrdinaryDiffEqCore.get_EEst(integrator))
     # TODO: add `isJcurrent` support for Rosenbrock solvers

lib/OrdinaryDiffEqDifferentiation/test/runtests.jl

@@ -35,6 +35,7 @@ if TEST_GROUP ∉ ("QA", "Sparse", "ModelingToolkit")
     @time @safetestset "Differentiation Trait Tests" include("differentiation_traits_tests.jl")
     @time @safetestset "Autodiff Error Tests" include("autodiff_error_tests.jl")
     @time @safetestset "No Jac Tests" include("nojac_tests.jl")
+    @time @safetestset "Stale W Linear Operator Tests" include("stale_w_linear_operator_tests.jl")
 end
 
 # Run sparse tests (separate environment due to ComponentArrays dep conflicts)

lib/OrdinaryDiffEqDifferentiation/test/stale_w_linear_operator_tests.jl

@@ -0,0 +1,44 @@
+# Regression tests for SciML/OrdinaryDiffEq.jl#3933: implicit solvers with an
+# AbstractSciMLOperator Jacobian (islinearfunction path) and a concrete-matrix
+# linear solver silently returned garbage, because do_newJW never requested a W
+# refactorization for linear operators (W = J - M/(γdt) changes with dt even
+# when J is constant). Requires SciMLOperators >= 1.24.4 for the companion fix
+# that rebuilds WOperator._concrete_form on convert for operator Jacobians.
+using OrdinaryDiffEqSDIRK, OrdinaryDiffEqDifferentiation, LinearSolve, SciMLOperators
+using SciMLBase, LinearAlgebra, Test
+
+@testset "stale W with linear operator f and concrete-A linsolve" begin
+    A = [-2.0 1.0; 1.0 -2.0]
+    B = [0.0 0.5; -0.5 0.0]
+    u0 = [1.0, 0.5]
+    tspan = (0.0, 1.0)
+    uref = exp((A + B) * tspan[2]) * u0
+
+    f2! = (du, u, p, t) -> mul!(du, B, u)
+    prob_split = SplitODEProblem(MatrixOperator(A), f2!, u0, tspan)
+
+    @testset "$(nameof(typeof(alg))) split, LU" for alg in (
+            KenCarp3(linsolve = LUFactorization()),
+            KenCarp4(linsolve = LUFactorization()),
+            KenCarp5(linsolve = LUFactorization()),
+        )
+        sol = solve(prob_split, alg; abstol = 1.0e-8, reltol = 1.0e-8)
+        @test SciMLBase.successful_retcode(sol)
+        @test norm(sol.u[end] - uref) / norm(uref) < 1.0e-3
+        # the broken path accepted the whole interval in ~2 giant steps
+        @test length(sol.t) > 5
+    end
+
+    # Non-split linear-operator problem hits the same do_newJW islin branch
+    prob_lin = ODEProblem(ODEFunction(MatrixOperator(A + B)), u0, tspan)
+    sol = solve(prob_lin, TRBDF2(linsolve = LUFactorization()); abstol = 1.0e-8, reltol = 1.0e-8)
+    @test SciMLBase.successful_retcode(sol)
+    @test norm(sol.u[end] - uref) / norm(uref) < 1.0e-3
+
+    # The Krylov path was always correct; make sure it stays that way
+    sol_krylov = solve(
+        prob_split, KenCarp4(linsolve = KrylovJL_GMRES());
+        abstol = 1.0e-8, reltol = 1.0e-8
+    )
+    @test norm(sol_krylov.u[end] - uref) / norm(uref) < 1.0e-3
+end