AKODE segfault with Sparse M matrix

[jd-lara]

Segfault message

signal (11): Segmentation fault: 11
in expression starting at REPL[25]:2
getindex at ./array.jl:801 [inlined]
_broadcast_getindex at ./broadcast.jl:614 [inlined]
_getindex at ./broadcast.jl:644 [inlined]
_broadcast_getindex at ./broadcast.jl:620 [inlined]
getindex at ./broadcast.jl:575 [inlined]
macro expansion at ./broadcast.jl:984 [inlined]
macro expansion at ./simdloop.jl:77 [inlined]
copyto! at ./broadcast.jl:983 [inlined]
copyto! at ./broadcast.jl:936 [inlined]
materialize! at ./broadcast.jl:894 [inlined]
materialize! at ./broadcast.jl:891 [inlined]
convert at /Users/jdlara/.julia/packages/Sundials/sB5CR/src/types_and_consts_additions.jl:50
massmat at /Users/jdlara/.julia/packages/Sundials/sB5CR/src/common_interface/function_types.jl:159
unknown function (ip: 0x15a6d26d2)

[rodrigomha]

I got this issue

ArgumentError: 1 == colptr[1] != 1

Stacktrace:
  [1] (::SparseArrays.var"#throwstart#3")(ckp::Int64)
    @ SparseArrays /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/sparsematrix.jl:58
  [2] sparse_check
    @ /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/sparsematrix.jl:63 [inlined]
  [3] SparseArrays.SparseMatrixCSC(m::Int64, n::Int64, colptr::Vector{Int64}, rowval::Vector{Int64}, nzval::Vector{Float64})
    @ SparseArrays /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/sparsematrix.jl:39
  [4] convert(#unused#::Type{SparseArrays.SparseMatrixCSC}, J::Ptr{Sundials._generic_SUNMatrix})
    @ Sundials ~/.julia/packages/Sundials/sB5CR/src/types_and_consts_additions.jl:54
  [5] massmat(t::Float64, _M::Ptr{Sundials._generic_SUNMatrix}, mmf::Sundials.FunJac{ODEFunction{true, typeof(PowerSimulationsDynamics.system_mass_matrix!), SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Nothing, Nothing, PowerSimulationsDynamics.SimulationInputs, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Vector{Float64}, Nothing, Nothing, Nothing}, tmp1::Ptr{Sundials._generic_N_Vector}, tmp2::Ptr{Sundials._generic_N_Vector}, tmp3::Ptr{Sundials._generic_N_Vector})
    @ Sundials ~/.julia/packages/Sundials/sB5CR/src/common_interface/function_types.jl:159
  [6] ARKStepEvolve
    @ ~/.julia/packages/Sundials/sB5CR/src/API/arkode.jl:416 [inlined]
  [7] ARKStepEvolve(arkode_mem::Sundials.Handle{Sundials.ARKStepMem}, tout::Float64, yout::Sundials.NVector, tret::Vector{Float64}, itask::Int32)
    @ Sundials ~/.julia/packages/Sundials/sB5CR/src/API/arkode.jl:421
  [8] solver_step(integrator::Sundials.ARKODEIntegrator{Sundials.NVector, PowerSimulationsDynamics.SimulationInputs, Sundials.Handle{Sundials.ARKStepMem}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Nothing, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, PowerSimulationsDynamics.SimulationInputs, ODEFunction{true, typeof(PowerSimulationsDynamics.system_mass_matrix!), SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ARKODE{:Newton, :Dense, :Dense, Sundials.Implicit, Nothing, Nothing, Nothing, Nothing}, SciMLBase.HermiteInterpolation{Vector{Float64}, Vector{Vector{Float64}}, Vector{Vector{Float64}}}, DiffEqBase.DEStats}, ARKODE{:Newton, :Dense, :Dense, Sundials.Implicit, Nothing, Nothing, Nothing, Nothing}, ODEFunction{true, typeof(PowerSimulationsDynamics.system_mass_matrix!), SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Sundials.FunJac{ODEFunction{true, typeof(PowerSimulationsDynamics.system_mass_matrix!), SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Nothing, Nothing, PowerSimulationsDynamics.SimulationInputs, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Vector{Float64}, Nothing, Nothing, Nothing}, Nothing, Sundials.DEOptions{DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, CallbackSet{Tuple{}, Tuple{DiscreteCallback{PowerSimulationsDynamics.var"#32#33"{NetworkSwitch}, PowerSimulationsDynamics.var"#7#8"{NetworkSwitch}, typeof(DiffEqBase.INITIALIZE_DEFAULT), typeof(DiffEqBase.FINALIZE_DEFAULT)}}}, Float64, Float64, typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE)}, Vector{Float64}, Tuple{Int64}, Vector{Float64}, Sundials.LinSolHandle{Sundials.Dense}, Sundials.MatrixHandle{Sundials.DenseMatrix}, Sundials.LinSolHandle{Sundials.Dense}, Sundials.MatrixHandle{Sundials.DenseMatrix}, Nothing}, tstop::Float64)
    @ Sundials ~/.julia/packages/Sundials/sB5CR/src/common_interface/solve.jl:1448
  [9] solve!(integrator::Sundials.ARKODEIntegrator{Sundials.NVector, PowerSimulationsDynamics.SimulationInputs, Sundials.Handle{Sundials.ARKStepMem}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Nothing, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, PowerSimulationsDynamics.SimulationInputs, ODEFunction{true, typeof(PowerSimulationsDynamics.system_mass_matrix!), SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ARKODE{:Newton, :Dense, :Dense, Sundials.Implicit, Nothing, Nothing, Nothing, Nothing}, SciMLBase.HermiteInterpolation{Vector{Float64}, Vector{Vector{Float64}}, Vector{Vector{Float64}}}, DiffEqBase.DEStats}, ARKODE{:Newton, :Dense, :Dense, Sundials.Implicit, Nothing, Nothing, Nothing, Nothing}, ODEFunction{true, typeof(PowerSimulationsDynamics.system_mass_matrix!), SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Sundials.FunJac{ODEFunction{true, typeof(PowerSimulationsDynamics.system_mass_matrix!), SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Nothing, Nothing, PowerSimulationsDynamics.SimulationInputs, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Vector{Float64}, Nothing, Nothing, Nothing}, Nothing, Sundials.DEOptions{DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, CallbackSet{Tuple{}, Tuple{DiscreteCallback{PowerSimulationsDynamics.var"#32#33"{NetworkSwitch}, PowerSimulationsDynamics.var"#7#8"{NetworkSwitch}, typeof(DiffEqBase.INITIALIZE_DEFAULT), typeof(DiffEqBase.FINALIZE_DEFAULT)}}}, Float64, Float64, typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE)}, Vector{Float64}, Tuple{Int64}, Vector{Float64}, Sundials.LinSolHandle{Sundials.Dense}, Sundials.MatrixHandle{Sundials.DenseMatrix}, Sundials.LinSolHandle{Sundials.Dense}, Sundials.MatrixHandle{Sundials.DenseMatrix}, Nothing})
    @ Sundials ~/.julia/packages/Sundials/sB5CR/src/common_interface/solve.jl:1525
 [10] __solve(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, PowerSimulationsDynamics.SimulationInputs, ODEFunction{true, typeof(PowerSimulationsDynamics.system_mass_matrix!), SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::ARKODE{:Newton, :Dense, :Dense, Sundials.Implicit, Nothing, Nothing, Nothing, Nothing}, timeseries::Vector{Any}, ts::Vector{Any}, ks::Vector{Any}, recompile::Type{Val{true}}; kwargs::Base.Iterators.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:callback, :tstops, :dtmax), Tuple{CallbackSet{Tuple{}, Tuple{DiscreteCallback{PowerSimulationsDynamics.var"#32#33"{NetworkSwitch}, PowerSimulationsDynamics.var"#7#8"{NetworkSwitch}, typeof(DiffEqBase.INITIALIZE_DEFAULT), typeof(DiffEqBase.FINALIZE_DEFAULT)}}}, Vector{Float64}, Float64}}})
    @ Sundials ~/.julia/packages/Sundials/sB5CR/src/common_interface/solve.jl:15
 [11] #solve_call#58
    @ ~/.julia/packages/DiffEqBase/lULzQ/src/solve.jl:61 [inlined]
 [12] solve_up(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, PowerSimulationsDynamics.SimulationInputs, ODEFunction{true, typeof(PowerSimulationsDynamics.system_mass_matrix!), SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{Float64}, p::PowerSimulationsDynamics.SimulationInputs, args::ARKODE{:Newton, :Dense, :Dense, Sundials.Implicit, Nothing, Nothing, Nothing, Nothing}; kwargs::Base.Iterators.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:callback, :tstops, :dtmax), Tuple{CallbackSet{Tuple{}, Tuple{DiscreteCallback{PowerSimulationsDynamics.var"#32#33"{NetworkSwitch}, PowerSimulationsDynamics.var"#7#8"{NetworkSwitch}, typeof(DiffEqBase.INITIALIZE_DEFAULT), typeof(DiffEqBase.FINALIZE_DEFAULT)}}}, Vector{Float64}, Float64}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/lULzQ/src/solve.jl:82
 [13] #solve#59
    @ ~/.julia/packages/DiffEqBase/lULzQ/src/solve.jl:70 [inlined]
 [14] execute!(sim::Simulation{MassMatrixModel}, solver::ARKODE{:Newton, :Dense, :Dense, Sundials.Implicit, Nothing, Nothing, Nothing, Nothing}; kwargs::Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:dtmax,), Tuple{Float64}}})
    @ PowerSimulationsDynamics ~/Repos_2021/PowerSimulationsDynamics.jl/src/base/simulation.jl:346
 [15] top-level scope
    @ In[24]:2
 [16] eval
    @ ./boot.jl:360 [inlined]
 [17] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1094

[ChrisRackauckas]

Sundials does weird things with sparse matrices that seems to lead to non-deterministic segfaults. For ARKODE, I would just recommend using the KenCarp methods from OrdinaryDiffEq.

[sjdaines]

The fix for #315 should I think also fix this one (might not be the whole story of course).

[ChrisRackauckas]

Awesome. @jd-lara and @rodrigomha follow up if it didn't.