# ModelingToolkit: Building a DAEProblem and solving it with IDA with KLU & Jacobian usage
**URL:** https://discourse.julialang.org/t/modelingtoolkit-building-a-daeproblem-and-solving-it-with-ida-with-klu-jacobian-usage/77814
**Category:** Modelling & Simulations
**Tags:** modelingtoolkit
**Created:** [March 13, 2022, 8:53am UTC](https://discourse.julialang.org/t/modelingtoolkit-building-a-daeproblem-and-solving-it-with-ida-with-klu-jacobian-usage/77814 "2022-03-13T08:53:57Z")
**Posts on this page:** 4
**Page:** 1
### Author:  [@RobertS](https://discourse.julialang.org/u/RobertS)
#### Post date: [March 13, 2022, 8:53am UTC](https://discourse.julialang.org/t/modelingtoolkit-building-a-daeproblem-and-solving-it-with-ida-with-klu-jacobian-usage/77814/1 "2022-03-13T08:53:57Z")
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:
```julia
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[i-1] - 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] + 5e-3 * 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` through `prob6` 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](https://discourse.julialang.org/t/rodas5-and-modelingtoolkit-much-slower-than-sundials-ida-on-large-systems/60712/1).
- 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](https://discourse.julialang.org/t/rodas5-and-modelingtoolkit-much-slower-than-sundials-ida-on-large-systems/60712/7).
Any help is greatly appreciated.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [March 14, 2022, 1:20am UTC](https://discourse.julialang.org/t/modelingtoolkit-building-a-daeproblem-and-solving-it-with-ida-with-klu-jacobian-usage/77814/2 "2022-03-14T01:20:21Z")
Everything beyond `prob3` is trying too hard. This is just a current issue with the `DAEProblem` constructor.
[https://github.com/SciML/ModelingToolkit.jl/issues/1482](https://github.com/SciML/ModelingToolkit.jl/issues/1482)
---
### Author:  [@bmetherall](https://discourse.julialang.org/u/bmetherall)
#### Post date: [March 15, 2022, 1:14pm UTC](https://discourse.julialang.org/t/modelingtoolkit-building-a-daeproblem-and-solving-it-with-ida-with-klu-jacobian-usage/77814/3 "2022-03-15T13:14:51Z")
I was able to get it working by specifying the Jacobian and sparsity in the undocumented(?) `DAEFunction`. My code is along the lines of
```julia
dae_f = DAEFunction(model!; jac = jac!, jac_prototype = build_sparse(n))
prob = DAEProblem(dae_f, dICs, ICs, tspan, params, differential_vars = differential_vars)
sol = solve(prob, IDA(linear_solver = :KLU))
```
In my implementation `jac!` numerically computes the Jacobian with `FiniteDiff.jl`, and `build_sparse(n)` specifies the sparsity of the Jacobian. You could probably do this automatically, but it was very slow for me when I tried it.
You will probably need to rewrite most of your code for it to be compatible with this setup, but I don’t know of any other way.
---
### Author:  [@RobertS](https://discourse.julialang.org/u/RobertS)
#### Post date: [March 20, 2022, 8:41am UTC](https://discourse.julialang.org/t/modelingtoolkit-building-a-daeproblem-and-solving-it-with-ida-with-klu-jacobian-usage/77814/4 "2022-03-20T08:41:13Z")
Thank you, I will take a deeper look on your proposed workaround.