Hi everyone,
I am trying to use the KLU method as linear solver for solving an ODE problem using Sundials.
I have tried to apply it to the CVODE_BDF solver for a simple problem, the rober equation. However, when I provide the sparsity pattern determined from the Jacobian function, that I also provide to the ODEFunction, it says:
“Sparsity Pattern in receiving SUNMatrix doesn’t match sending SparseMatrix”
The code is:
using Sundials
using IncompleteLU
using SparseArrays
using Plots
using LinearAlgebra################## Solving ODE ##################
ode repressentation
function rober!(du, u, p, t)
y₁, y₂, y₃ = u
k₁, k₂, k₃ = p
du[1] = -k₁ * y₁ + k₃ * y₂ * y₃
du[2] = k₁ * y₁ - k₂ * y₂^2 - k₃ * y₂ * y₃
du[3] = k₂ * y₂^2
nothing
endu0 = [1.0; 0.0; 0.0]
tspan = (0.0, 1e4)
p = [0.04, 3e7, 1e4]
prob = ODEProblem(rober!, u0, tspan,p)
sol = solve(prob,CVODE_BDF())
p1 = plot();
plot!(p1,sol,label=“ODE”);################## Solving ODE With Jacobian ##################
#Jacobian
function rober_jac!(J, u, p, t)
y₁, y₂, y₃ = u
k₁, k₂, k₃ = p
J[1, 1] = k₁ * -1
J[2, 1] = k₁
J[3, 1] = 0
J[1, 2] = y₃ * k₃
J[2, 2] = y₂ * k₂ * -2 + y₃ * k₃ * -1
J[3, 2] = y₂ * 2 * k₂
J[1, 3] = k₃ * y₂
J[2, 3] = k₃ * y₂ * -1
J[3, 3] = 0
nothing
endu0 = [1.0, 0, 0]
tspan = (0.0, 1e4)
p = [0.04, 3e7, 1e4]
jac_proto = zeros(length(u0),length(u0))
rober_jac!(jac_proto, u0 .+ 1e-6, p, 0.0)
fun = ODEFunction(rober!; jac = rober_jac!, jac_prototype = sparse(jac_proto))
prob = ODEProblem(fun, u0, tspan,p)
sol = solve(prob,CVODE_BDF(; linear_solver = :KLU))
p1 = plot();
plot!(p1,sol,label=“ODE”);
As seen from the ODE definition and the Jacobian, the sparsity patters is indeed:
x x x
x x x
⋅ x ⋅
However, if I instead set the sparsity patters as a dense 3x3 matrix, it solves the problem without any warning.
fun = ODEFunction(rober!; jac = rober_jac!, jac_prototype = sparse(ones(3,3)))
prob = ODEProblem(fun, u0, tspan,p)
sol = solve(prob,CVODE_BDF(; linear_solver = :KLU))
p1 = plot();
plot!(p1,sol,label=“ODE”);
To me, that seems a bit odd as the elements [3,1] and [3,3] in the jacobian prototype are indeed zeros as also seen from the ODE definition (Rober!) and the Jacobian function (Rober_jac!).
Did I specify anything incorrectly?
I also looked at the examples from Sundials.jl/test/common_interface/jacobians.jl at master · SciML/Sundials.jl (github.com), however, the Jacobian for the Lotka problem is in this case a dense 2x2 matrix.
Did some of you guys try to use the KLU linear solver in a Sundials solver, either ODE or DAE, and specify a sparsity pattern that is not a dense matrix?