don't apply MM twice, simpler test

Author: jClugstor

lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl

@@ -177,13 +177,14 @@ function calc_J!(J, integrator, cache, next_step::Bool = false)
 
             # we need to set all nzval to a non-zero number
             # otherwise in the following line any zero gets interpreted as a structural zero
-            integrator.f.jac_prototype.nzval .= 1.0
-            J .= 1.0 .* integrator.f.jac_prototype
-            J.nzval .= 0.0
-            f.jac(J, uprev, p, t)
-            MM = integrator.f.mass_matrix isa UniformScaling ?
-                 integrator.f.mass_matrix(length(integrator.u)) : integrator.f.mass_matrix
-            J .= J .+ MM
+            if !isnothing(integrator.f.jac_prototype)
+                integrator.f.jac_prototype.nzval .= 1.0
+                J .= 1.0 .* integrator.f.jac_prototype
+                J.nzval .= 0.0
+                f.jac(J, duprev, uprev, p, uf.α * uf.invγdt, t)
+            else
+                f.jac(J, duprev, uprev, p, uf.α * uf.invγdt, t)
+            end
         else
             @unpack du1, uf, jac_config = cache
             # using `dz` as temporary array
@@ -199,12 +200,14 @@ function calc_J!(J, integrator, cache, next_step::Bool = false)
 
             # we need to set all nzval to a non-zero number
             # otherwise in the following line any zero gets interpreted as a structural zero
-            integrator.f.jac_prototype.nzval .= 1.0
-            J .= 1.0 .* integrator.f.jac_prototype
-            J.nzval .= 0.0
-            f.jac(J, uprev, p, t)
-            MM = integrator.f.mass_matrix isa UniformScaling ? integrator.f.mass_matrix(length(integrator.u)) : integrator.f.mass_matrix
-            J .= J .+ MM
+            if !isnothing(integrator.f.jac_prototype)
+                integrator.f.jac_prototype.nzval .= 1.0
+                J .= 1.0 .* integrator.f.jac_prototype
+                J.nzval .= 0.0
+                f.jac(J, uprev, p, t)
+            else 
+                f.jac(J, uprev, p, t)
+            end
         else
             @unpack du1, uf, jac_config = cache
             uf.f = nlsolve_f(f, alg)

test/interface/sparsediff_tests.jl

@@ -87,35 +87,24 @@ end
 
 # test for https://github.com/SciML/OrdinaryDiffEq.jl/issues/2653#issuecomment-2778430025
 
-function sparse_f!(du, u, p, t)
-    du[1] = u[1] + u[2]
-    du[2] = u[3]^2
-    return du[3] = u[1]^2
-end
-
-backend = AutoSparse(
-    AutoForwardDiff();
-    sparsity_detector = TracerSparsityDetector(),
-    coloring_algorithm = GreedyColoringAlgorithm()
-)
-
-u = ones(3)
-du = zero(u)
-p = t = nothing
-
-prep = DI.prepare_jacobian(
-    sparse_f!, du, backend, u, DI.Constant(p), DI.Constant(t))
-# this is what the user may typically provide to the ODE problem
+using LinearAlgebra, SparseArrays
+using OrdinaryDiffEq
 
-function inplace_jac!(J, u, p, t)
-    return DI.jacobian!(
-        sparse_f!, zeros(3), J, prep, backend, u, DI.Constant(p), DI.Constant(t))
+function f(du, u, p, t)
+    du[1] = u[1]
+    return du
 end
 
-jac_prototype = similar(sparsity_pattern(prep), eltype(u))
+function jac(J::SparseMatrixCSC, u, p, t)
+    @assert nnz(J) == 1  # mirrors the strict behavior of SparseMatrixColorings
+    nonzeros(J)[1] = 1
+    return J
+end
 
-ode_f = ODEFunction(sparse_f!, jac = inplace_jac!, jac_prototype = jac_prototype)
-prob = ODEProblem(ode_f, [1, 1, 1], (0.0, 1.0))
+u0 = ones(10)
+jac_prototype = sparse(Diagonal(vcat(1, zeros(9))))
 
-@test_no_warn sol = solve(prob, Rodas5())
+fun = ODEFunction(f; jac, jac_prototype)
+prob = ODEProblem(fun, u0, (0.0, 1.0))
+@test_nowarn sol = solve(prob, Rodas4(); reltol = 1e-8, abstol = 1e-8)