I have coded an MTK system that I build a DAEProblem
from. I am now trying to find out how to solve that problem using Sundials IDA with KLU as linear solver. As far as I know, this requires that the Jacobian matrix is given. Here is what I tried:
using OrdinaryDiffEq, ModelingToolkit
using DASSL
using DASKR
using DifferentialEquations
using BenchmarkTools
using Plots
using Sundials
function build_system(n)
@parameters t # or would @variables t be better!?
@variables Tgas[1:n](t) Tsol[1:n](t) cR[1:n](t) mflux(t)
@parameters params[1:3]
D = Differential(t)
dp, Tin, h = params
dx = h / n
############################################################
# Build equations ##########################################
############################################################
eqns = Array{Equation}(undef, 3 * n + 1)
count = 1
# specification as function or intermediate variabel does not matter. I could also make this an
# actual @variable, then I could retrieve it later.
rR = 5.5e6 .* exp.(.5680 ./(273 .+ Tsol)) .* (1 . exp.(.cR./50))
q = 1e4 * (Tgas  Tsol)
# gas energy balance
eqns[count] = 0 ~ mflux * 1000 * (Tin  Tgas[1])  q[1] * dx
count += 1
for i in 2:n
eqns[count] = 0 ~ mflux * 1000 * (Tgas[i1]  Tgas[i])  q[i] * dx
count += 1
end
# solids energy balance
for i in 1:n
eqns[count] = D(Tsol[i]) ~ q[i] / (1000 * 2000)
count += 1
end
# gas momentum balance
eqns[count] = 0 ~  mflux[1] + 5e3 * sqrt(dp)
count += 1
# reactant concentration equations
for i in 1:n
eqns[count] = D(cR[i]) ~ rR[i]
count += 1
end
@named sys = ODESystem(eqns, t, vcat(Tgas, Tsol, cR, [mflux]), params)
return sys
end
dp = 100e2
Tin = 200
h = 0.5
parameters = [dp, Tin, h]
N = 100
tspan = (0.0, 1500.0)
sys = build_system(N)
# quick and dirty definition of u0, should be Pairs instead
u0 = ones(length(sys.states))
u0[1:N] .= parameters[2]
u0[N+1 : 2*N] .= 20
u0[2*N+1 : 3*N] .= 3000
u0[3*N+1] = 0.5
du0 = 0 .* copy(u0)
prob2 = DAEProblem(sys, du0, u0, tspan, parameters)
prob3 = DAEProblem(sys, du0, u0, tspan, parameters, jac=true, sparse=true)
jac = eval(ModelingToolkit.generate_jacobian(sys)[2])
prob4 = DAEProblem(sys, du0, u0, tspan, parameters, jac=jac, sparse=true)
diff_vars = vcat(zeros(N),ones(N), ones(N), [0])
prob5 = DAEProblem(sys, du0, u0, tspan, parameters, differential_vars=diff_vars, jac=jac, sparse=true)
prob6 = DAEProblem(sys, du0, u0, tspan, parameters, differential_vars=diff_vars, jac=true, sparse=true)
# 124ms, 54285 allocs
@benchmark solve(prob2, IDA())
# 120ms, 54289 allocs
@benchmark solve(prob3, IDA())
# 123ms, 54289 allocs
@benchmark solve(prob4, IDA())
# 125ms, 54289 allocs
@benchmark solve(prob5, IDA())
# 124ms, 54289 allocs
@benchmark solve(prob6, IDA())
# Errors with: MethodError: no method matching nonzeros(::Nothing)
@benchmark solve(prob5, IDA(linear_solver = :KLU))
I have two problems with that code:

First thing that is suspicious is that all solves on
prob2
throughprob6
seem to do pretty much the same with practically exact same performance and allocations, although the problems are set up in quite different ways (with/without jac and sparsity, with/without explicit jac, etc). This finding is similar to what was discussed a year ago here. 
Second thing is to solve the
DAEProblem
with IDA using:KLU
as linear solver. What do I have to do to get this working? As you see at the end of my code, my trials are erroring with the same error message as reported earlier here.
Any help is greatly appreciated.