search

SciML / Sundials.jl

Julia interface to Sundials, including a nonlinear solver (KINSOL), ODE's (CVODE and ARKODE), and DAE's (IDA) in a SciML scientific machine learning enabled manner
https://diffeq.sciml.ai
BSD 2-Clause "Simplified" License
209 stars 79 forks source link

Fix KLU and sparse matrix support for Sundials 5 #244

Closed ChrisRackauckas closed 4 years ago

ChrisRackauckas commented 4 years ago

https://github.com/SciML/Sundials.jl/pull/235#discussion_r403608473

ChrisRackauckas commented 4 years ago

The issue stems from the Julia v1.3 changes

https://github.com/JuliaLang/julia/blob/v1.3.0/NEWS.md#sparsearrays

MWE:

using Sundials, Test, SparseArrays, DiffEqOperators

# Test for Jacobian usage
function Lotka(du, u, p, t)
  du[1] = u[1] - u[1] * u[2] # REPL[7], line 3:
  du[2] = -3 * u[2] + 1 * u[1] * u[2]
  nothing
end

jac_called = false
function Lotka_jac(J,u,p,t)
  global jac_called
  jac_called = true
  J[1,1] = 1.0 - u[2]
  J[1,2] = -u[1]
  J[2,1] = 1 * u[2]
  J[2,2] = -3 + u[1]
  nothing
end

Lotka_f = ODEFunction(Lotka,jac=Lotka_jac,
                      jac_prototype = sparse([1,2,1,2],[1,1,2,2],zeros(4)))

prob = ODEProblem(Lotka_f,ones(2),(0.0,10.0))
jac_called = false
sol9 = solve(prob,CVODE_BDF(linear_solver=:KLU))
jac_called == true

Res:

julia> sol9 = solve(prob,CVODE_BDF(linear_solver=:KLU))
1 = 1
2 = 2
ERROR: ArgumentError: 1 == colptr[1] != 1
Stacktrace:
 [1] sparse_check(::Int64, ::Array{Int64,1}, ::Array{Int64,1}, ::Array{Float64,1}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/SparseArrays/src/sparsematrix.jl:58
 [2] SparseMatrixCSC(::Int64, ::Int64, ::Array{Int64,1}, ::Array{Int64,1}, ::Array{Float64,1}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/SparseArrays/src/sparsematrix.jl:39
 [3] convert(::Type{SparseMatrixCSC}, ::Ptr{Sundials._generic_SUNMatrix}) at /home/scheme/.julia/dev/Sundials/src/types_and_consts_additions.jl:49
 [4] cvodejac(::Float64, ::Ptr{Sundials._generic_N_Vector}, ::Ptr{Sundials._generic_N_Vector}, ::Ptr{Sundials._generic_SUNMatrix}, ::Sundials.FunJac{ODEFunction{true,typeof(Lotka),UniformScaling{Bool},Nothing,Nothing,typeof(Lotka_jac),Nothing,Nothing,SparseMatrixCSC{Float64,Int64},Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,typeof(Lotka_jac),DiffEqBase.NullParameters,Nothing,SparseMatrixCSC{Float64,Int64},Array{Float64,1},Nothing,Nothing,Nothing}, ::Ptr{Sundials._generic_N_Vector}, ::Ptr{Sundials._generic_N_Vector}, ::Ptr{Sundials._generic_N_Vector}) at /home/scheme/.julia/dev/Sundials/src/common_interface/function_types.jl:66
 [5] CVode at /home/scheme/.julia/dev/Sundials/src/API/cvodes.jl:166 [inlined]
 [6] CVode(::Sundials.Handle{Sundials.CVODEMem}, ::Float64, ::Array{Float64,1}, ::Array{Float64,1}, ::Int32) at /home/scheme/.julia/dev/Sundials/src/API/cvodes.jl:171
 [7] solver_step(::Sundials.CVODEIntegrator{Array{Float64,1},DiffEqBase.NullParameters,Sundials.Handle{Sundials.CVODEMem},ODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Nothing,ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(Lotka),UniformScaling{Bool},Nothing,Nothing,typeof(Lotka_jac),Nothing,Nothing,SparseMatrixCSC{Float64,Int64},Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},CVODE_BDF{:Newton,:KLU,Nothing,Nothing},DiffEqBase.HermiteInterpolation{Array{Float64,1},Array{Array{Float64,1},1},Array{Array{Float64,1},1}},DiffEqBase.DEStats},CVODE_BDF{:Newton,:KLU,Nothing,Nothing},ODEFunction{true,typeof(Lotka),UniformScaling{Bool},Nothing,Nothing,typeof(Lotka_jac),Nothing,Nothing,SparseMatrixCSC{Float64,Int64},Nothing,Nothing,Nothing,Nothing,Nothing},Sundials.FunJac{ODEFunction{true,typeof(Lotka),UniformScaling{Bool},Nothing,Nothing,typeof(Lotka_jac),Nothing,Nothing,SparseMatrixCSC{Float64,Int64},Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,typeof(Lotka_jac),DiffEqBase.NullParameters,Nothing,SparseMatrixCSC{Float64,Int64},Array{Float64,1},Nothing,Nothing,Nothing},Ptr{Nothing},Sundials.DEOptions{DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},CallbackSet{Tuple{},Tuple{}},Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE)},Array{Float64,1},Tuple{Int64},Array{Float64,1},Sundials.LinSolHandle{Sundials.KLU},Sundials.MatrixHandle{Sundials.SparseMatrix},Nothing}, ::Float64) at /home/scheme/.julia/dev/Sundials/src/common_interface/solve.jl:1097
 [8] solve!(::Sundials.CVODEIntegrator{Array{Float64,1},DiffEqBase.NullParameters,Sundials.Handle{Sundials.CVODEMem},ODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Nothing,ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(Lotka),UniformScaling{Bool},Nothing,Nothing,typeof(Lotka_jac),Nothing,Nothing,SparseMatrixCSC{Float64,Int64},Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},CVODE_BDF{:Newton,:KLU,Nothing,Nothing},DiffEqBase.HermiteInterpolation{Array{Float64,1},Array{Array{Float64,1},1},Array{Array{Float64,1},1}},DiffEqBase.DEStats},CVODE_BDF{:Newton,:KLU,Nothing,Nothing},ODEFunction{true,typeof(Lotka),UniformScaling{Bool},Nothing,Nothing,typeof(Lotka_jac),Nothing,Nothing,SparseMatrixCSC{Float64,Int64},Nothing,Nothing,Nothing,Nothing,Nothing},Sundials.FunJac{ODEFunction{true,typeof(Lotka),UniformScaling{Bool},Nothing,Nothing,typeof(Lotka_jac),Nothing,Nothing,SparseMatrixCSC{Float64,Int64},Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,typeof(Lotka_jac),DiffEqBase.NullParameters,Nothing,SparseMatrixCSC{Float64,Int64},Array{Float64,1},Nothing,Nothing,Nothing},Ptr{Nothing},Sundials.DEOptions{DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},CallbackSet{Tuple{},Tuple{}},Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE)},Array{Float64,1},Tuple{Int64},Array{Float64,1},Sundials.LinSolHandle{Sundials.KLU},Sundials.MatrixHandle{Sundials.SparseMatrix},Nothing}) at /home/scheme/.julia/dev/Sundials/src/common_interface/solve.jl:1151
 [9] __solve(::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(Lotka),UniformScaling{Bool},Nothing,Nothing,typeof(Lotka_jac),Nothing,Nothing,SparseMatrixCSC{Float64,Int64},Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::CVODE_BDF{:Newton,:KLU,Nothing,Nothing}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/scheme/.julia/dev/Sundials/src/common_interface/solve.jl:12
 [10] __solve at /home/scheme/.julia/dev/Sundials/src/common_interface/solve.jl:10 [inlined] (repeats 5 times)
 [11] (::DiffEqBase.var"#461#462"{ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(Lotka),UniformScaling{Bool},Nothing,Nothing,typeof(Lotka_jac),Nothing,Nothing,SparseMatrixCSC{Float64,Int64},Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},Tuple{CVODE_BDF{:Newton,:KLU,Nothing,Nothing}}})() at /home/scheme/.julia/packages/DiffEqBase/k3AhB/src/solve.jl:49
 [12] maybe_with_logger at /home/scheme/.julia/packages/DiffEqBase/k3AhB/src/utils.jl:259 [inlined]
 [13] solve_call(::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(Lotka),UniformScaling{Bool},Nothing,Nothing,typeof(Lotka_jac),Nothing,Nothing,SparseMatrixCSC{Float64,Int64},Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::CVODE_BDF{:Newton,:KLU,Nothing,Nothing}; merge_callbacks::Bool, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/scheme/.julia/packages/DiffEqBase/k3AhB/src/solve.jl:48
 [14] solve_call at /home/scheme/.julia/packages/DiffEqBase/k3AhB/src/solve.jl:37 [inlined]
 [15] #solve#463 at /home/scheme/.julia/packages/DiffEqBase/k3AhB/src/solve.jl:66 [inlined]
 [16] solve(::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(Lotka),UniformScaling{Bool},Nothing,Nothing,typeof(Lotka_jac),Nothing,Nothing,SparseMatrixCSC{Float64,Int64},Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::CVODE_BDF{:Newton,:KLU,Nothing,Nothing}) at /home/scheme/.julia/packages/DiffEqBase/k3AhB/src/solve.jl:54
 [17] top-level scope at REPL[206]:1
 [18] eval(::Module, ::Any) at ./boot.jl:331
 [19] eval_user_input(::Any, ::REPL.REPLBackend) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/REPL/src/REPL.jl:86
 [20] run_backend(::REPL.REPLBackend) at /home/scheme/.julia/packages/Revise/Pcs5V/src/Revise.jl:1073
 [21] top-level scope at none:0

julia> sol9 = solve(prob,CVODE_BDF(linear_solver=:KLU))
julia> sol9 = solve(prob,CVODE_BDF(linear_solver=:KLU))
ERROR: ArgumentError: 1 == colptr[1] != 1
Stacktrace:
 [1] sparse_check(::Int64, ::Array{Int64,1}, ::Array{Int64,1}, ::Array{Float64,1}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/SparseArrays/src/sparsematrix.jl:58
 [2] SparseMatrixCSC(::Int64, ::Int64, ::Array{Int64,1}, ::Array{Int64,1}, ::Array{Float64,1}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/SparseArrays/src/sparsematrix.jl:39
 [3] convert(::Type{SparseMatrixCSC}, ::Ptr{Sundials._generic_SUNMatrix}) at /home/scheme/.julia/dev/Sundials/src/types_and_consts_additions.jl:49
 [4] cvodejac(::Float64, ::Ptr{Sundials._generic_N_Vector}, ::Ptr{Sundials._generic_N_Vector}, ::Ptr{Sundials._generic_SUNMatrix}, ::Sundials.FunJac{ODEFunction{true,typeof(Lotka),UniformScaling{Bool},Nothing,Nothing,typeof(Lotka_jac),Nothing,Nothing,SparseMatrixCSC{Float64,Int64},Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,typeof(Lotka_jac),DiffEqBase.NullParameters,Nothing,SparseMatrixCSC{Float64,Int64},Array{Float64,1},Nothing,Nothing,Nothing}, ::Ptr{Sundials._generic_N_Vector}, ::Ptr{Sundials._generic_N_Vector}, ::Ptr{Sundials._generic_N_Vector}) at /home/scheme/.julia/dev/Sundials/src/common_interface/function_types.jl:66
 [5] CVode at /home/scheme/.julia/dev/Sundials/src/API/cvodes.jl:166 [inlined]
 [6] CVode(::Sundials.Handle{Sundials.CVODEMem}, ::Float64, ::Array{Float64,1}, ::Array{Float64,1}, ::Int32) at /home/scheme/.julia/dev/Sundials/src/API/cvodes.jl:171
 [7] solver_step(::Sundials.CVODEIntegrator{Array{Float64,1},DiffEqBase.NullParameters,Sundials.Handle{Sundials.CVODEMem},ODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Nothing,ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(Lotka),UniformScaling{Bool},Nothing,Nothing,typeof(Lotka_jac),Nothing,Nothing,SparseMatrixCSC{Float64,Int64},Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},CVODE_BDF{:Newton,:KLU,Nothing,Nothing},DiffEqBase.HermiteInterpolation{Array{Float64,1},Array{Array{Float64,1},1},Array{Array{Float64,1},1}},DiffEqBase.DEStats},CVODE_BDF{:Newton,:KLU,Nothing,Nothing},ODEFunction{true,typeof(Lotka),UniformScaling{Bool},Nothing,Nothing,typeof(Lotka_jac),Nothing,Nothing,SparseMatrixCSC{Float64,Int64},Nothing,Nothing,Nothing,Nothing,Nothing},Sundials.FunJac{ODEFunction{true,typeof(Lotka),UniformScaling{Bool},Nothing,Nothing,typeof(Lotka_jac),Nothing,Nothing,SparseMatrixCSC{Float64,Int64},Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,typeof(Lotka_jac),DiffEqBase.NullParameters,Nothing,SparseMatrixCSC{Float64,Int64},Array{Float64,1},Nothing,Nothing,Nothing},Ptr{Nothing},Sundials.DEOptions{DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},CallbackSet{Tuple{},Tuple{}},Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE)},Array{Float64,1},Tuple{Int64},Array{Float64,1},Sundials.LinSolHandle{Sundials.KLU},Sundials.MatrixHandle{Sundials.SparseMatrix},Nothing}, ::Float64) at /home/scheme/.julia/dev/Sundials/src/common_interface/solve.jl:1095
 [8] solve!(::Sundials.CVODEIntegrator{Array{Float64,1},DiffEqBase.NullParameters,Sundials.Handle{Sundials.CVODEMem},ODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Nothing,ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(Lotka),UniformScaling{Bool},Nothing,Nothing,typeof(Lotka_jac),Nothing,Nothing,SparseMatrixCSC{Float64,Int64},Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},CVODE_BDF{:Newton,:KLU,Nothing,Nothing},DiffEqBase.HermiteInterpolation{Array{Float64,1},Array{Array{Float64,1},1},Array{Array{Float64,1},1}},DiffEqBase.DEStats},CVODE_BDF{:Newton,:KLU,Nothing,Nothing},ODEFunction{true,typeof(Lotka),UniformScaling{Bool},Nothing,Nothing,typeof(Lotka_jac),Nothing,Nothing,SparseMatrixCSC{Float64,Int64},Nothing,Nothing,Nothing,Nothing,Nothing},Sundials.FunJac{ODEFunction{true,typeof(Lotka),UniformScaling{Bool},Nothing,Nothing,typeof(Lotka_jac),Nothing,Nothing,SparseMatrixCSC{Float64,Int64},Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,typeof(Lotka_jac),DiffEqBase.NullParameters,Nothing,SparseMatrixCSC{Float64,Int64},Array{Float64,1},Nothing,Nothing,Nothing},Ptr{Nothing},Sundials.DEOptions{DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},CallbackSet{Tuple{},Tuple{}},Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE)},Array{Float64,1},Tuple{Int64},Array{Float64,1},Sundials.LinSolHandle{Sundials.KLU},Sundials.MatrixHandle{Sundials.SparseMatrix},Nothing}) at /home/scheme/.julia/dev/Sundials/src/common_interface/solve.jl:1149
 [9] __solve(::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(Lotka),UniformScaling{Bool},Nothing,Nothing,typeof(Lotka_jac),Nothing,Nothing,SparseMatrixCSC{Float64,Int64},Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::CVODE_BDF{:Newton,:KLU,Nothing,Nothing}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/scheme/.julia/dev/Sundials/src/common_interface/solve.jl:12
 [10] __solve at /home/scheme/.julia/dev/Sundials/src/common_interface/solve.jl:10 [inlined] (repeats 5 times)
 [11] (::DiffEqBase.var"#461#462"{ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(Lotka),UniformScaling{Bool},Nothing,Nothing,typeof(Lotka_jac),Nothing,Nothing,SparseMatrixCSC{Float64,Int64},Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},Tuple{CVODE_BDF{:Newton,:KLU,Nothing,Nothing}}})() at /home/scheme/.julia/packages/DiffEqBase/k3AhB/src/solve.jl:49
 [12] maybe_with_logger at /home/scheme/.julia/packages/DiffEqBase/k3AhB/src/utils.jl:259 [inlined]
 [13] solve_call(::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(Lotka),UniformScaling{Bool},Nothing,Nothing,typeof(Lotka_jac),Nothing,Nothing,SparseMatrixCSC{Float64,Int64},Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::CVODE_BDF{:Newton,:KLU,Nothing,Nothing}; merge_callbacks::Bool, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/scheme/.julia/packages/DiffEqBase/k3AhB/src/solve.jl:48
 [14] solve_call at /home/scheme/.julia/packages/DiffEqBase/k3AhB/src/solve.jl:37 [inlined]
 [15] #solve#463 at /home/scheme/.julia/packages/DiffEqBase/k3AhB/src/solve.jl:66 [inlined]
 [16] solve(::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(Lotka),UniformScaling{Bool},Nothing,Nothing,typeof(Lotka_jac),Nothing,Nothing,SparseMatrixCSC{Float64,Int64},Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::CVODE_BDF{:Newton,:KLU,Nothing,Nothing}) at /home/scheme/.julia/packages/DiffEqBase/k3AhB/src/solve.jl:54
 [17] top-level scope at REPL[206]:1
 [18] eval(::Module, ::Any) at ./boot.jl:331
 [19] eval_user_input(::Any, ::REPL.REPLBackend) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/REPL/src/REPL.jl:86
 [20] run_backend(::REPL.REPLBackend) at /home/scheme/.julia/packages/Revise/Pcs5V/src/Revise.jl:1073
 [21] top-level scope at none:0
ChrisRackauckas commented 4 years ago

It seems like we need to move the colptr when converting, and that was done before using this code:

https://github.com/SciML/Sundials.jl/blob/master/src/common_interface/function_types.jl#L57-L77 https://github.com/SciML/Sundials.jl/blob/master/src/types_and_consts_additions.jl#L38-L50