OrdinaryDiffEq with a mass_matrix and jac_prototype

Specific Domains Numerics
1

Hi everyone,

Short story : see the MWE below.

Using OrdinaryDiffEq, I’m trying to solve an ODE involving a mass_matrix that (may) evolve depending on the solution, the time, the parameters. So I set up a SciMLOperators.MatrixOperator and provided it to the ODEFunction. So far, so good.

The problem I’m facing is that if I provide a jac_prototype (of my rhs function) to the ODEFunction, then the solve errors here OrdinaryDiffEq.jl/lib/OrdinaryDiffEqDifferentiation/src/alg_utils.jl at b80bd3b002927c63ba1d0eec9fdcbcee3d560873 · SciML/OrdinaryDiffEq.jl · GitHub. Without jac_prototype, on a MWE, everything is fine. Am I doing something wrong?

Moreover, I don’t understand the quoted lines in OrdinaryDiffEqDifferentiation : the jac_prototype is related to my rhs function, why are mass_matrix and jac_prototype connected together here?

Here is a MWE:

using OrdinaryDiffEq
using SparseArrays
using SparseDiffTools

const M_mwe = sparse([1.0 2.0; 0.0 1.0])

function update_func(A, u, p, t)
    return M_mwe
end

function run()
    mass_matrix = MatrixOperator(M_mwe; update_func = update_func)
    f!(du, u, p, t) = du .= u
    alg = ImplicitEuler()
    tspan = (0.0, 1.0)
    u0 = ones(2)

    # Without jac_prototype
    odeFunction = ODEFunction(f!; mass_matrix)
    prob = ODEProblem(odeFunction, u0, tspan)
    solve(prob, alg) #-------> OK

    # With jac_prototype
    jac_prototype = sparse([1.0 0.0; 0.0 1.0])
    colorvec = [1, 1]
    odeFunction = ODEFunction(f!; mass_matrix, jac_prototype, colorvec)
    prob = ODEProblem(odeFunction, u0, tspan)
    solve(prob, alg) #-------> KO
end

run()

PS : passing the jac_prototype of the rhs is not an option for me because in my real case the automatic tools (such as SparseConnectivityTracer) fail to find it.

2

Hi!
Can you give the error stack trace?
Also, what do you get when you try automatic sparsity detection? Just for curiosity

3

Hi Guillaume, here is the stack trace:

MethodError: no method matching keys(::MatrixOperator{Float64, SparseMatrixCSC{Float64, Int64}, SciMLOperators.FilterKwargs{typeof(update_func), Tuple{}}, SciMLOperators.FilterKwargs{typeof(SciMLOperators.DEFAULT_UPDATE_FUNC), Tuple{}}})
The function `keys` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  keys(::Cmd)
   @ Base process.jl:716
  keys(::Base.TermInfo)
   @ Base terminfo.jl:232
  keys(::LibGit2.GitTree)
   @ Revise .julia/packages/Revise/p9Zlq/src/git.jl:52
  ...

Stacktrace:
  [1] pairs(collection::MatrixOperator{Float64, SparseMatrixCSC{Float64, Int64}, SciMLOperators.FilterKwargs{typeof(update_func), Tuple{}}, SciMLOperators.FilterKwargs{typeof(SciMLOperators.DEFAULT_UPDATE_FUNC), Tuple{}}})
    @ Base ./abstractdict.jl:178
  [2] findall(testf::ComposedFunction{typeof(!), typeof(iszero)}, A::MatrixOperator{Float64, SparseMatrixCSC{Float64, Int64}, SciMLOperators.FilterKwargs{typeof(update_func), Tuple{}}, SciMLOperators.FilterKwargs{typeof(SciMLOperators.DEFAULT_UPDATE_FUNC), Tuple{}}})
    @ Base ./array.jl:2700
  [3] prepare_user_sparsity(ad_alg::AutoForwardDiff{nothing, ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}}, prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, var"#f!#run##1", MatrixOperator{Float64, SparseMatrixCSC{Float64, Int64}, SciMLOperators.FilterKwargs{typeof(update_func), Tuple{}}, SciMLOperators.FilterKwargs{typeof(SciMLOperators.DEFAULT_UPDATE_FUNC), Tuple{}}}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem})
    @ OrdinaryDiffEqDifferentiation .julia/packages/OrdinaryDiffEqDifferentiation/DE4z5/src/alg_utils.jl:120
  [4] prepare_alg(alg::ImplicitEuler{0, AutoForwardDiff{nothing, Nothing}, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Nothing}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!)}, u0::Vector{Float64}, p::SciMLBase.NullParameters, prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, var"#f!#run##1", MatrixOperator{Float64, SparseMatrixCSC{Float64, Int64}, SciMLOperators.FilterKwargs{typeof(update_func), Tuple{}}, SciMLOperators.FilterKwargs{typeof(SciMLOperators.DEFAULT_UPDATE_FUNC), Tuple{}}}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem})
    @ OrdinaryDiffEqDifferentiation .julia/packages/OrdinaryDiffEqDifferentiation/DE4z5/src/alg_utils.jl:53
  [5] solve_up(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, var"#f!#run##1", MatrixOperator{Float64, SparseMatrixCSC{Float64, Int64}, SciMLOperators.FilterKwargs{typeof(update_func), Tuple{}}, SciMLOperators.FilterKwargs{typeof(SciMLOperators.DEFAULT_UPDATE_FUNC), Tuple{}}}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{Float64}, p::SciMLBase.NullParameters, args::ImplicitEuler{0, AutoForwardDiff{nothing, Nothing}, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Nothing}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!)}; kwargs::@Kwargs{})
    @ DiffEqBase .julia/packages/DiffEqBase/qvEPa/src/solve.jl:1200
  [6] solve_up(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, var"#f!#run##1", MatrixOperator{Float64, SparseMatrixCSC{Float64, Int64}, SciMLOperators.FilterKwargs{typeof(update_func), Tuple{}}, SciMLOperators.FilterKwargs{typeof(SciMLOperators.DEFAULT_UPDATE_FUNC), Tuple{}}}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{Float64}, p::SciMLBase.NullParameters, args::ImplicitEuler{0, AutoForwardDiff{nothing, Nothing}, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Nothing}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!)})
    @ DiffEqBase .julia/packages/DiffEqBase/qvEPa/src/solve.jl:1183
  [7] solve(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, var"#f!#run##1", MatrixOperator{Float64, SparseMatrixCSC{Float64, Int64}, SciMLOperators.FilterKwargs{typeof(update_func), Tuple{}}, SciMLOperators.FilterKwargs{typeof(SciMLOperators.DEFAULT_UPDATE_FUNC), Tuple{}}}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem}, args::ImplicitEuler{0, AutoForwardDiff{nothing, Nothing}, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Nothing}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!)}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{true}, kwargs::@Kwargs{})
    @ DiffEqBase .julia/packages/DiffEqBase/qvEPa/src/solve.jl:1096
  [8] solve(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, var"#f!#run##1", MatrixOperator{Float64, SparseMatrixCSC{Float64, Int64}, SciMLOperators.FilterKwargs{typeof(update_func), Tuple{}}, SciMLOperators.FilterKwargs{typeof(SciMLOperators.DEFAULT_UPDATE_FUNC), Tuple{}}}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Vector{Int64}, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem}, args::ImplicitEuler{0, AutoForwardDiff{nothing, Nothing}, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Nothing}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!)})
    @ DiffEqBase .julia/packages/DiffEqBase/qvEPa/src/solve.jl:1086
  [9] run()
    @ Main ./REPL[11]:18
 [10] top-level scope
    @ REPL[12]:1

Regarding the automatic sparsity detection, the problem is this one : Bcube and DifferentiationInterface · Issue #189 · bcube-project/Bcube.jl · GitHub. I could use your trick, but as a first step I’d prefer not to.

1 Like
4

If I recall correctly, this probably has something to do with the fact that the jac_prototype and sparsity pattern need to have certain entries that correspond to the non-zeroes of the mass matrix set to be potentially non-zero. It has something to do with the way the WOperator is calculated.

Probably a related issue: Strict in-place sparse Jacobians now cause errors when they don't anticipate the mass matrix · Issue #2653 · SciML/OrdinaryDiffEq.jl · GitHub

Can you file an issue on the github repo?

5

Yes, I read again the fonction and indeed the jacobian of the whole system M dq/dt=f(x) is something like M - dt J xhere J=df/dq.

Ok I’ll open an issue on the repo with this MWE and in the meantime I will write the specialization for findall(!iszero, M::MatrixOperator in my script.

1 Like