Solving DAEs with known sparsity pattern and IDA

General Usage
, ,
1

Hi everyone!

Can someone provide a minimal example of how to use sparsity patterns (jac_prototype) when solving DAEs with IDA solver? I manage to use sparsity patterns with ODEs and the julia solvers…

Here are some examples that I am trying to solve:

ODE Version

n = 1000
L = 0.5 .+ rand(n)

function odefun(yp, y, p, t)
yp[:] = -L.*y
end

jac_pattern = spdiagm(ones(n))
colorvec = matrix_colors(jac_pattern)
tspan = (0,2)
y0 = ones(n)
ode = ODEFunction(odefun)
ode_sparse = ODEFunction(odefun, jac_prototype=jac_pattern, colorvec=colorvec)

prob = ODEProblem(ode, y0, tspan)
prob_sparse = ODEProblem(ode_sparse, y0, tspan)

@btime solve(prob, TRBDF2()) # solves in 85 ms
@btime solve(prob_sparse, TRBDF2()) # solves in 1 ms

However the DAE doesn’t work:

n = 1000
L = 0.5 .+ rand(n)

function residuals(res, yp, y, p, t)
res[:] = yp + L.*y
end

jac_pattern = spdiagm(ones(n))
colorvec = matrix_colors(jac_pattern)
tspan = (0,2)
y0 = ones(n)
yp0 = -L

dae = DAEFunction(residuals)
dae_sparse = DAEFunction(residuals, jac_prototype=jac_pattern, colorvec=colorvec)

prob = DAEProblem(dae, yp0, y0, tspan, differential_vars=trues(n))
prob_sparse = DAEProblem(dae_sparse, yp0, y0, tspan, differential_vars=trues(n))

@btime solve(prob, IDA()) # solves in 52 ms
@btime solve(prob_sparse, IDA()) # solves in 49 ms

I would expect that the sparse problem would have a huge decrease in time as in the ODE case.

Thanks for the support!

2

The same problem occurs when using ODE formulation with CVODE_BDF also from sundials:

n = 1000
L = 0.5 .+ rand(n)

function odefun(yp, y, p, t)
yp[:] = -L.*y
end

jac_pattern = spdiagm(ones(n))
colorvec = matrix_colors(jac_pattern)
tspan = (0,2)
y0 = ones(n)

ode = ODEFunction(odefun)
ode_sparse = ODEFunction(odefun, jac_prototype=jac_pattern, colorvec=colorvec)
prob = ODEProblem(ode, y0, tspan)
prob_sparse = ODEProblem(ode_sparse, y0, tspan)

@btime solve(prob, CVODE_BDF()) # solves in 27 ms
@btime solve(prob_sparse, CVODE_BDF()) # solves in 27ms

3

With Sundials you have to tell it to use a sparse linear solver, i.e. CVODE_BDF(linear_solver=:KLU)

4

Thank you for the answer, Cris!

However, for the KLU you have to pass the actual jacobian…

5

Oh yes, Sundials cannot use sparse differentiation.