I'm trying to solve a set of stiff ordinary differential equations using DifferentialEquations.jl. It works fine when I'm just providing the ODE functions with or without the Jacobian, however when I'm trying to provide the sparsity pattern using the jac_prototype keyword it always results in a BoundsError.
ERROR: LoadError: BoundsError
Stacktrace:
[1] _copyto_impl!(dest::Vector{Float64}, doffs::Int64, src::Vector{Float64}, soffs::Int64, n::Int64)
@ Base ./array.jl:329
[2] copyto!
@ ./array.jl:322 [inlined]
[3] copyto!(dest::Vector{Float64}, src::Vector{Float64})
@ Base ./array.jl:346
[4] solve(cache::LinearSolve.LinearCache{SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, KLUFactorization, KLU.KLUFactorization{Float64, Int64}, LinearSolve.InvPreconditioner{Diagonal{Float64, Vector{Float64}}}, Diagonal{Float64, Vector{Float64}}, Float64}, alg::KLUFactorization; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:reltol,), Tuple{Float64}}})
@ LinearSolve ~/.julia/packages/LinearSolve/FDu3j/src/factorization.jl:296
[5] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{CompositeAlgorithm{Tuple{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{8, false, KLUFactorization, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}, true, Vector{Float64}, Nothing, Float64, SciMLBase.NullParameters, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, OrdinaryDiffEq.ODECompositeSolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(swing!), UniformScaling{Bool}, Nothing, Nothing, typeof(jacobian!), Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, CompositeAlgorithm{Tuple{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{8, false, KLUFactorization, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}, OrdinaryDiffEq.CompositeInterpolationData{ODEFunction{true, typeof(swing!), UniformScaling{Bool}, Nothing, Nothing, typeof(jacobian!), Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.Rosenbrock23Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, OrdinaryDiffEq.Rosenbrock23Tableau{Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, typeof(swing!), UniformScaling{Bool}, Nothing, Nothing, typeof(jacobian!), Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Float64}, SciMLBase.NullParameters}, SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(swing!), UniformScaling{Bool}, Nothing, Nothing, typeof(jacobian!), Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, SciMLBase.NullParameters}, LinearSolve.LinearCache{SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, KLUFactorization, KLU.KLUFactorization{Float64, Int64}, LinearSolve.InvPreconditioner{Diagonal{Float64, Vector{Float64}}}, Diagonal{Float64, Vector{Float64}}, Float64}, Nothing, FiniteDiff.GradientCache{Nothing, Vector{Float64}, Vector{Float64}, Float64, Val{:forward}(), Float64, Val{true}()}, Float64, Rosenbrock23{8, false, KLUFactorization, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Nothing}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}}, DiffEqBase.DEStats}, ODEFunction{true, typeof(swing!), UniformScaling{Bool}, Nothing, Nothing, typeof(jacobian!), Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.Rosenbrock23Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, OrdinaryDiffEq.Rosenbrock23Tableau{Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, typeof(swing!), UniformScaling{Bool}, Nothing, Nothing, typeof(jacobian!), Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Float64}, SciMLBase.NullParameters}, SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(swing!), UniformScaling{Bool}, Nothing, Nothing, typeof(jacobian!), Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, SciMLBase.NullParameters}, LinearSolve.LinearCache{SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, KLUFactorization, KLU.KLUFactorization{Float64, Int64}, LinearSolve.InvPreconditioner{Diagonal{Float64, Vector{Float64}}}, Diagonal{Float64, Vector{Float64}}, Float64}, Nothing, FiniteDiff.GradientCache{Nothing, Vector{Float64}, Vector{Float64}, Float64, Val{:forward}(), Float64, Val{true}()}, Float64, Rosenbrock23{8, false, KLUFactorization, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Nothing}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.Rosenbrock23Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, OrdinaryDiffEq.Rosenbrock23Tableau{Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, typeof(swing!), UniformScaling{Bool}, Nothing, Nothing, typeof(jacobian!), Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Float64}, SciMLBase.NullParameters}, SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(swing!), UniformScaling{Bool}, Nothing, Nothing, typeof(jacobian!), Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, SciMLBase.NullParameters}, LinearSolve.LinearCache{SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, KLUFactorization, KLU.KLUFactorization{Float64, Int64}, LinearSolve.InvPreconditioner{Diagonal{Float64, Vector{Float64}}}, Diagonal{Float64, Vector{Float64}}, Float64}, Nothing, FiniteDiff.GradientCache{Nothing, Vector{Float64}, Vector{Float64}, Float64, Val{:forward}(), Float64, Val{true}()}, Float64, Rosenbrock23{8, false, KLUFactorization, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Nothing}, repeat_step::Bool)
@ OrdinaryDiffEq ~/.julia/packages/LinearSolve/FDu3j/src/factorization.jl:0
[6] perform_step!
@ ~/.julia/packages/OrdinaryDiffEq/ccwlJ/src/perform_step/composite_perform_step.jl:52 [inlined]
[7] perform_step!
@ ~/.julia/packages/OrdinaryDiffEq/ccwlJ/src/perform_step/composite_perform_step.jl:49 [inlined]
[8] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{CompositeAlgorithm{Tuple{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{8, false, KLUFactorization, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}, true, Vector{Float64}, Nothing, Float64, SciMLBase.NullParameters, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, OrdinaryDiffEq.ODECompositeSolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(swing!), UniformScaling{Bool}, Nothing, Nothing, typeof(jacobian!), Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, CompositeAlgorithm{Tuple{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{8, false, KLUFactorization, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}, OrdinaryDiffEq.CompositeInterpolationData{ODEFunction{true, typeof(swing!), UniformScaling{Bool}, Nothing, Nothing, typeof(jacobian!), Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.Rosenbrock23Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, OrdinaryDiffEq.Rosenbrock23Tableau{Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, typeof(swing!), UniformScaling{Bool}, Nothing, Nothing, typeof(jacobian!), Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Float64}, SciMLBase.NullParameters}, SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(swing!), UniformScaling{Bool}, Nothing, Nothing, typeof(jacobian!), Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, SciMLBase.NullParameters}, LinearSolve.LinearCache{SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, KLUFactorization, KLU.KLUFactorization{Float64, Int64}, LinearSolve.InvPreconditioner{Diagonal{Float64, Vector{Float64}}}, Diagonal{Float64, Vector{Float64}}, Float64}, Nothing, FiniteDiff.GradientCache{Nothing, Vector{Float64}, Vector{Float64}, Float64, Val{:forward}(), Float64, Val{true}()}, Float64, Rosenbrock23{8, false, KLUFactorization, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Nothing}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}}, DiffEqBase.DEStats}, ODEFunction{true, typeof(swing!), UniformScaling{Bool}, Nothing, Nothing, typeof(jacobian!), Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.Rosenbrock23Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, OrdinaryDiffEq.Rosenbrock23Tableau{Float64}, SciMLBase.TimeGradientWrapper{ODEFunction{true, typeof(swing!), UniformScaling{Bool}, Nothing, Nothing, typeof(jacobian!), Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Float64}, SciMLBase.NullParameters}, SciMLBase.UJacobianWrapper{ODEFunction{true, typeof(swing!), UniformScaling{Bool}, Nothing, Nothing, typeof(jacobian!), Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, SciMLBase.NullParameters}, LinearSolve.LinearCache{SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, KLUFactorization, KLU.KLUFactorization{Float64, Int64}, LinearSolve.InvPreconditioner{Diagonal{Float64, Vector{Float64}}}, Diagonal{Float64, Vector{Float64}}, Float64}, Nothing, FiniteDiff.GradientCache{Nothing, Vector{Float64}, Vector{Float64}, Float64, Val{:forward}(), Float64, Val{true}()}, Float64, Rosenbrock23{8, false, KLUFactorization, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Nothing}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/ccwlJ/src/solve.jl:475
[9] #__solve#562
@ ~/.julia/packages/OrdinaryDiffEq/ccwlJ/src/solve.jl:5 [inlined]
[10] __solve(::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(swing!), UniformScaling{Bool}, Nothing, Nothing, typeof(jacobian!), Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ::Nothing; default_set::Bool, kwargs::Base.Pairs{Symbol, Integer, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:second_time, :progress, :progress_steps), Tuple{Bool, Bool, Int64}}})
@ DifferentialEquations ~/.julia/packages/DifferentialEquations/TZ7Qg/src/default_solve.jl:9
[11] #__solve#48
@ ~/.julia/packages/DiffEqBase/GHOAd/src/solve.jl:1049 [inlined]
[12] #solve_call#24
@ ~/.julia/packages/DiffEqBase/GHOAd/src/solve.jl:451 [inlined]
[13] solve_up(::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, typeof(swing!), UniformScaling{Bool}, Nothing, Nothing, typeof(jacobian!), Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ::Nothing, ::Vector{Float64}, ::SciMLBase.NullParameters; kwargs::Base.Pairs{Symbol, Integer, Tuple{Symbol, Symbol}, NamedTuple{(:progress, :progress_steps), Tuple{Bool, Int64}}})
@ DiffEqBase ~/.julia/packages/DiffEqBase/GHOAd/src/solve.jl:818
[14] #solve#29
@ ~/.julia/packages/DiffEqBase/GHOAd/src/solve.jl:784 [inlined]
[15] top-level scope
@ ~/Documents/juliatest/disctest2.jl:59
in expression starting at ~/Documents/juliatest/disctest2.jl:59
This is happening for both the case when providing an example Jacobian for the sparsity pattern and when using the Symbolics.jacobian_sparsity function as described in the documentation. The problem can be reproduced using the following code. I am using Julia v1.8.0 on MacOS with DifferentialEquations.jl v7.2.0.
using DifferentialEquations
using LinearAlgebra
using SparseArrays
edgelist = [1 5 10.; 1 3 10.; 5 3 10.; 3 4 3.; 4 2 12.; 4 6 12.; 2 6 12.]
m = [1; 1; 0; 0; 0; 0]
d = [0.1; 0.1; 0.1; 0.1; 0.1; 0.1]
P = [1.7; 1.5; -0.6; -0.8; -0.8; -1.0]
N = 6
Niner = 2
function swing!(du, u, p, t)
du[1:Niner] = u[N + 1:end]
du[Niner + 1 : N] = P[Niner+1:N]
du[N + 1 : end] = P[1 : Niner] - d[1:Niner] .* u[N+1:end]
for row in eachrow(edgelist)
_i = Int(row[1])
_j = Int(row[2])
_k = row[3]
_i1 = _i > Niner ? _i : _i + N
_i2 = _j > Niner ? _j : _j + N
du[_i1] -= _k * sin(u[_i] - u[_j])
du[_i2] -= _k * sin(u[_j] - u[_i])
end
du[Niner + 1 : N] ./= d[Niner+1:N]
du[N + 1 : end] ./= m[1:Niner]
end
function jacobian!(J, u, p, t)
J[1:Niner, N+1:end] = sparse(I, Niner, Niner)
J[N+1:end, N+1:end] = -spdiagm(d[1:Niner]./m[1:Niner])
L = spzeros(N, N)
for row in eachrow(edgelist)
_i = Int(row[1])
_j = Int(row[2])
_k = row[3] * cos(u[_i] - u[_j])
L[_i, _j] = _k
L[_j, _i] = _k
end
for i=1:N
L[i, i] = -sum(L[i, :])
end
J[Niner+1:N, 1:N] = diagm(1 ./ d[Niner+1:N]) * L[Niner+1:end, :]
J[N+1:end, 1:N] = diagm(1 ./ m[1:Niner]) * L[1:Niner, :]
end
u0 = [0.06866403; 0.02864531; -0.00473835; -0.03807784; -0.00807496; -0.0464182; zeros(Niner)]
# u0 = zeros(12)
jacproto = spzeros(N+Niner, N+Niner)
jacobian!(jacproto, ones(N+Niner), 0, 0)
tspan = (0., 25.)
func = ODEFunction(swing!, jac=jacobian!, jac_prototype=jacproto)
prob = ODEProblem(func, u0, tspan)
solve(prob, progress=true, progress_steps=1)
The jacproto doesn't is an 8x8 matrix,
julia> jacproto
8×8 SparseMatrixCSC{Float64, Int64} with 24 stored entries:
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.0 ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.0
100.0 ⋅ -230.0 30.0 100.0 ⋅ ⋅ ⋅
⋅ 120.0 30.0 -270.0 ⋅ 120.0 ⋅ ⋅
100.0 ⋅ 100.0 ⋅ -200.0 ⋅ ⋅ ⋅
⋅ 120.0 ⋅ 120.0 ⋅ -240.0 ⋅ ⋅
-20.0 ⋅ 10.0 ⋅ 10.0 ⋅ -0.1 ⋅
⋅ -24.0 ⋅ 12.0 ⋅ 12.0 ⋅ -0.1
Which does not fully represent the diagonal. Remember that the sparsity pattern is the prototype sparsity of I - gamma*J for J = f', and so you need to account for all of the diagonal. Adding the diagonal terms like:
jacobian!(jacproto, ones(N+Niner), 0, 0)
jacproto[1,1] = 1
jacproto[2,2] = 1
tspan = (0., 25.)
func = ODEFunction(swing!, jac=jacobian!, jac_prototype=jacproto)
prob = ODEProblem(func, u0, tspan)
solve(prob, progress=true, progress_steps=1)
makes the solve successful with KLU.