Unsupported interpolator chosen when mass matrix is not diagonal

[termi-official]

Describe the bug 🐞

Solver chooses unsupported interpolation (Hermite) when mass matrix is not diagonal.

Expected behavior

The solver should either fall back to a supported interpolation (e.g. linear interpolation)

Also, users should be able to specify the used interpolator manually.

Minimal Reproducible Example 👇

using OrdinaryDiffEq, SparseArrays, LinearAlgebra

Δt₀ = 0.1
Δt_save = 0.025
T = 1.0

M = sparse([1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 3, 4, 5, 6, 3, 4, 5, 6], [1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6], [0.012527777777777773, 0.0062638888888888865, 0.0031319444444444437, 0.0062638888888888865, 0.0062638888888888865, 0.012527777777777773, 0.0062638888888888865, 0.0031319444444444437, 0.0031319444444444437, 0.0062638888888888865, 0.025055555555555546, 0.012527777777777773, 0.006263888888888887, 0.003131944444444444, 0.0062638888888888865, 0.0031319444444444437, 0.012527777777777773, 0.025055555555555546, 0.003131944444444444, 0.006263888888888887, 0.006263888888888887, 0.003131944444444444, 0.012527777777777775, 0.006263888888888888, 0.003131944444444444, 0.006263888888888887, 0.006263888888888888, 0.012527777777777775], 6, 6)
K = sparse([1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 3, 4, 5, 6, 3, 4, 5, 6], [1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6], [-1.018551367331855, -0.3229120473022911, 0.5092756836659275, 0.8321877309682186, -0.3229120473022911, -1.018551367331855, 0.8321877309682185, 0.5092756836659275, 0.5092756836659275, 0.8321877309682185, -2.037102734663709, -0.6458240946045821, 0.8321877309682181, 0.5092756836659272, 0.8321877309682186, 0.5092756836659275, -0.6458240946045821, -2.0371027346637094, 0.5092756836659272, 0.8321877309682181, 0.8321877309682181, 0.5092756836659272, -1.0185513673318543, -0.32291204730229095, 0.5092756836659272, 0.8321877309682181, -0.32291204730229095, -1.0185513673318543], 6, 6)
u₀ = rand(6)

@show cond(Matrix(M - Δt₀*K))

struct RHSparams
    K::SparseMatrixCSC
end
p = RHSparams(K)

function heat!(du,u,p,t)
    @show "rhs",t
    
    du .= p.K * u
end;

function heat_jac!(J,u,p,t)
    @show "jac",t
    J .= p.K
end;

rhs = ODEFunction(heat!, mass_matrix=M; jac=heat_jac!, jac_prototype=copy(K))
problem = ODEProblem(rhs, u₀, (0.0,T), p);
timestepper = ImplicitEuler()

solve(problem, timestepper, dt=Δt₀, saveat=Δt_save);

Error & Stacktrace ⚠️

ERROR: Hermite interpolation is not defined in this case. The Hermite interpolation
fallback only supports diagonal mass matrices. If you have a DAE with a
non-diagonal mass matrix, then the dense output is not supported with this
ODE solver. Either use a method which has a specialized interpolation,
such as Rodas5P, or use `dense=false`

You can find the list of available DAE solvers with their documented interpolations at:
https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae_solve/


Some of the types have been truncated in the stacktrace for improved reading. To emit complete information
in the stack trace, evaluate `TruncatedStacktraces.VERBOSE[] = true` and re-run the code.

Stacktrace:
  [1] _ode_interpolant!(out::Vector{…}, Θ::Float64, dt::Float64, y₀::Vector{…}, y₁::Vector{…}, k::Vector{…}, cache::OrdinaryDiffEq.ImplicitEulerCache{…}, idxs::Nothing, T::Type{…}, differential_vars::OrdinaryDiffEq.DifferentialVarsUndefined)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/AfyyB/src/dense/generic_dense.jl:672
  [2] ode_interpolant(Θ::Float64, dt::Float64, y₀::Vector{…}, y₁::Vector{…}, k::Vector{…}, cache::OrdinaryDiffEq.ImplicitEulerCache{…}, idxs::Nothing, T::Type{…}, differential_vars::OrdinaryDiffEq.DifferentialVarsUndefined)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/AfyyB/src/dense/generic_dense.jl:602
  [3] ode_interpolant
    @ ~/.julia/packages/OrdinaryDiffEq/AfyyB/src/dense/generic_dense.jl:114 [inlined]
  [4] _savevalues!(integrator::OrdinaryDiffEq.ODEIntegrator{…}, force_save::Bool, reduce_size::Bool)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/AfyyB/src/integrators/integrator_utils.jl:73
  [5] savevalues! (repeats 2 times)
    @ ~/.julia/packages/OrdinaryDiffEq/AfyyB/src/integrators/integrator_utils.jl:58 [inlined]
  [6] handle_callbacks!
    @ ~/.julia/packages/OrdinaryDiffEq/AfyyB/src/integrators/integrator_utils.jl:353 [inlined]
  [7] _loopfooter!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/AfyyB/src/integrators/integrator_utils.jl:254
  [8] loopfooter!
    @ ~/.julia/packages/OrdinaryDiffEq/AfyyB/src/integrators/integrator_utils.jl:207 [inlined]
  [9] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/AfyyB/src/solve.jl:545
 [10] __solve(::ODEProblem{…}, ::ImplicitEuler{…}; kwargs::@Kwargs{…})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/AfyyB/src/solve.jl:7
 [11] __solve
    @ ~/.julia/packages/OrdinaryDiffEq/AfyyB/src/solve.jl:1 [inlined]
 [12] solve_call(_prob::ODEProblem{…}, args::ImplicitEuler{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/WyGjp/src/solve.jl:612
 [13] solve_up(prob::ODEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::RHSparams, args::ImplicitEuler{…}; kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/WyGjp/src/solve.jl:1080
 [14] solve_up
    @ ~/.julia/packages/DiffEqBase/WyGjp/src/solve.jl:1066 [inlined]
 [15] solve(prob::ODEProblem{…}, args::ImplicitEuler{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/WyGjp/src/solve.jl:1003
 [16] top-level scope
    @ REPL[99]:1

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

Status `~/Repos/Ferrite.jl/docs/Project.toml`
  [8e7c35d0] BlockArrays v0.16.43
  [e30172f5] Documenter v1.4.1
  [daee34ce] DocumenterCitations v1.3.3
  [c061ca5d] Ferrite v0.3.14 `..`
  [4f95f4f8] FerriteGmsh v1.1.0
  [0f8c756f] FerriteMeshParser v0.1.8
  [f6369f11] ForwardDiff v0.10.36
  [705231aa] Gmsh v0.3.1
  [42fd0dbc] IterativeSolvers v0.9.4
  [d3d80556] LineSearches v7.2.0
  [7ed4a6bd] LinearSolve v2.30.0
  [98b081ad] Literate v2.18.0
  [16fef848] LiveServer v1.3.1
  [429524aa] Optim v1.9.4
  [1dea7af3] OrdinaryDiffEq v6.76.0 
  [91a5bcdd] Plots v1.40.4
  [92933f4c] ProgressMeter v1.10.0
  [90137ffa] StaticArrays v1.9.3
  [48a634ad] Tensors v1.16.1
  [a759f4b9] TimerOutputs v0.5.23
  [3a884ed6] UnPack v1.0.2