# Rodas5 and ModelingToolkit much slower than Sundials' IDA on large systems

I’m trying to solve a 1-D advection-diffusion-reaction equation which has been discretized into a set of ODEs. My code below is about 30x slower for n~1000 than a similar implementation as a DAEProblem solved with IDA(). Are there any obvious reasons why this should be so much slower?

``````using OrdinaryDiffEq, ModelingToolkit

function build_system(n)
###
# Declarations
###

@parameters t
@variables carbon[1:n](t) CO[1:n](t) alpha(t) xi(t)
@parameters params[1:7]
D = Differential(t)

por, r1, cald, calD, calU, calR, calS = params # Give useful names to the parameters
Pe = calU / calD

r = LinRange(r1, 1, n)
dr = (1 - r1) / (n - 1)

lambda_plus(ri) = cald * calD * por / dr * (ri + dr / 2)^2 - calU * alpha / 2
lambda_naught(ri, carbon) = cald * calD * por / dr * ((ri + dr / 2)^2 + (ri - dr / 2)^2) + por * (1 - por) * carbon * ri^2 * dr
lambda_minus(ri) = cald * calD * por / dr * (ri - dr / 2)^2 + calU * alpha / 2

###
# Build equations
###

eqns = Array{Equation}(undef, 2 * n + 2)
count = 1

# Carbon layer time update
for i in 1:n
eqns[count] = D(carbon[i]) ~ - por * carbon[i] * CO[i]
count += 1
end

# Inner concentration of CO
eqns[count] = 0 ~ CO[1] + calR / Pe - (alpha / xi^2 + calR / Pe) * exp(alpha * Pe * (1 / xi - 1 / r1))
count += 1

# Inner flux of CO
eqns[count] = 0 ~ (lambda_plus(r1) + lambda_minus(r1)) * CO[2] - lambda_minus(r1) * (CO[1] + calR / Pe) * alpha * Pe / r1^2 * 2 * dr / (cald * por) - lambda_naught(r1, carbon[1]) * CO[1]
count += 1

# Interior CO
for i in 2:n-1
eqns[count] = 0 ~ lambda_plus(r[i]) * CO[i+1] - lambda_naught(r[i], carbon[i]) * CO[i] + lambda_minus(r[i]) * CO[i-1]
count += 1
end

# Matching condition
eqns[count] = 0 ~ (lambda_plus(r[n]) + lambda_minus(r[n])) * CO[n-1] + lambda_plus(r[n]) * 2.0 * alpha * Pe * dr / (cald * por) * (CO[n] - 1) / (1 - exp(alpha * Pe)) - lambda_naught(r[n], carbon[n]) * CO[n]
count += 1

# Quartz boundary
eqns[count] = D(xi) ~ - calU * calS * alpha / xi^2

return ODESystem(eqns, t, vcat(carbon, CO, [alpha, xi]), params)
end

function solve_system(sys; tspan = (0.0, 1.5), param_values = [0.3, 5.0 / 6.0, 1.0, 1.0, 1.0, 1.0, 1.0])
ICs = ones(length(sys.states))
ICs[end] = param_values[2]

#ModelingToolkit.generate_jacobian(sys)[2]
#jac = eval(ModelingToolkit.generate_jacobian(sys)[2])

prob = ODEProblem(ODEFunction(sys), ICs, tspan, param_values)

@time sol = solve(prob, Rodas5(autodiff = false), reltol = 1e-10, abstol = 1e-6);

return sol
end

###
# Compile
###

solve_system(build_system(5))

###
# Proper
###

sys = build_system(101)
sol = solve_system(sys)
``````

I know that supplying the Jacobian and sparsity pattern should make it a lot faster, but nothing I do seems to actually be using the analytic Jacobian. I’ve tried uncommenting the `generate_jacobian` lines in `solve_system`, which I think should add the Jacobian to `sys`, but the time is the same. I’ve also tried passing it directly to `ODEFunction`, but it expects a boolean, and to `ODEProblem`, but it doesn’t seem to make a difference. I’ve also tried with `autodiff = true`, but it’s just much slower.

Would it be better to just use Sundials to solve this, or am I missing something that is making my code much slower?

As mentioned in the documentation, Rosenbrock methods do not scale to large systems of equations very well.

At medium tolerances ( `>1e-8` ?) it is recommended you use `Rodas5` , `Rodas4P` (the former is more efficient but the later is more reliable), `Kvaerno5` , or `KenCarp4` .

For asymptotically large systems of ODEs ( `N>1000` ?) where `f` is very costly and the complex eigenvalues are minimal (low oscillations), in that case `CVODE_BDF` will be the most efficient but requires `Vector{Float64}` . `CVODE_BDF` will also do surprisingly well if the solution is smooth.

replace CVODE_BDF with IDA since it has the same properties. We should probably update the DAE solver page to say that. The cutoffs aren’t very precise for many reasons (it’s system and problem dependent), but around 1000 ODEs or so is around the time to switch away from Rosenbrock methods. At your 2,000, I would definitely not be surprised to see IDA outperforming.

One method to try could be RadauIIA5, though IDa is probably what you want in this instance.

With Rosenbrock methods you’d probably want to `tgrad=true` since it needs a separate derivative for the computation.

6 Likes

Thanks for the comments, I guess IDA is the way to go. Is there a way to define a DAEProblem in MTK to use with Sundials (I can’t really find any documentation about this)? When I naively try

``````	prob = DAEProblem(DAEFunction(sys), ICs, tspan, param_values)

@time sol = solve(prob, IDA(), reltol = 1e-10, abstol = 1e-6);
``````

I get the error `ERROR: LoadError: MethodError: no method matching size(::Tuple{Float64, Float64})`.

Full Trace
``````ERROR: LoadError: MethodError: no method matching size(::Tuple{Float64, Float64})
Closest candidates are:
size(::Tuple, ::Integer) at tuple.jl:27
size(::Union{LinearAlgebra.QR, LinearAlgebra.QRCompactWY, LinearAlgebra.QRPivoted}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/qr.jl:524
size(::Union{LinearAlgebra.QR, LinearAlgebra.QRCompactWY, LinearAlgebra.QRPivoted}, ::Integer) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/qr.jl:523
...
Stacktrace:
[1] (DAEProblem{true, duType, tType, isinplace, P, F, K, D} where {duType, tType, isinplace, P, F, K, D})(f::DAEFunction{true, ModelingToolkit.var"#f#160"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#576"), Symbol("##arg#577"), Symbol("##arg#578"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x77493f90, 0x4499b563, 0xe171eec5, 0xf4ea53b4, 0x48cc079b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#579"), Symbol("##arg#576"), Symbol("##arg#577"), Symbol("##arg#578"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x3967f6b7, 0xd19b559f, 0xc7f890a0, 0xde6739f5, 0x3c2fea5e)}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing}, du0::Vector{Float64}, u0::Tuple{Float64, Float64}, tspan::Vector{Float64}, p::SciMLBase.NullParameters; differential_vars::Nothing, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ SciMLBase ~/.julia/packages/SciMLBase/9EjAY/src/problems/dae_problems.jl:18
[2] (DAEProblem{true, duType, tType, isinplace, P, F, K, D} where {duType, tType, isinplace, P, F, K, D})(f::DAEFunction{true, ModelingToolkit.var"#f#160"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#576"), Symbol("##arg#577"), Symbol("##arg#578"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x77493f90, 0x4499b563, 0xe171eec5, 0xf4ea53b4, 0x48cc079b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#579"), Symbol("##arg#576"), Symbol("##arg#577"), Symbol("##arg#578"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x3967f6b7, 0xd19b559f, 0xc7f890a0, 0xde6739f5, 0x3c2fea5e)}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing}, du0::Vector{Float64}, u0::Tuple{Float64, Float64}, tspan::Vector{Float64}, p::SciMLBase.NullParameters)
@ SciMLBase ~/.julia/packages/SciMLBase/9EjAY/src/problems/dae_problems.jl:18
[3] DAEProblem(f::DAEFunction{true, ModelingToolkit.var"#f#160"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#576"), Symbol("##arg#577"), Symbol("##arg#578"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x77493f90, 0x4499b563, 0xe171eec5, 0xf4ea53b4, 0x48cc079b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#579"), Symbol("##arg#576"), Symbol("##arg#577"), Symbol("##arg#578"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x3967f6b7, 0xd19b559f, 0xc7f890a0, 0xde6739f5, 0x3c2fea5e)}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing}, du0::Vector{Float64}, u0::Tuple{Float64, Float64}, tspan::Vector{Float64}, p::SciMLBase.NullParameters; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ SciMLBase ~/.julia/packages/SciMLBase/9EjAY/src/problems/dae_problems.jl:37
[4] DAEProblem(f::DAEFunction{true, ModelingToolkit.var"#f#160"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#576"), Symbol("##arg#577"), Symbol("##arg#578"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x77493f90, 0x4499b563, 0xe171eec5, 0xf4ea53b4, 0x48cc079b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#579"), Symbol("##arg#576"), Symbol("##arg#577"), Symbol("##arg#578"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x3967f6b7, 0xd19b559f, 0xc7f890a0, 0xde6739f5, 0x3c2fea5e)}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing}, du0::Vector{Float64}, u0::Tuple{Float64, Float64}, tspan::Vector{Float64}, p::SciMLBase.NullParameters) (repeats 2 times)
@ SciMLBase ~/.julia/packages/SciMLBase/9EjAY/src/problems/dae_problems.jl:37
[5] solve_system(sys::ODESystem; tspan::Tuple{Float64, Float64}, param_values::Vector{Float64})
@ Main ~/Desktop/Silicon_Processes/Numerics/Example.jl:67
[6] solve_system(sys::ODESystem)
@ Main ~/Desktop/Silicon_Processes/Numerics/Example.jl:61
[7] top-level scope
@ ~/Desktop/Silicon_Processes/Numerics/Example.jl:78
[8] include(fname::String)
@ Base.MainInclude ./client.jl:444
[9] top-level scope
@ REPL[25]:1
``````

If not, would following something like
https://discourse.julialang.org/t/error-evaluating-sparse-analytical-jacobian-from-modelingtoolkit-in-dae-system/45764/2
be best for optimizing the Jacobian calculation and exploiting the sparsity?

Just do it the easy way.

``````prob = DAEProblem(sys, u0, du0, tspan, param_values)
``````

Right, I was forgetting the `du0` argument. I have a few more questions @ChrisRackauckas, if you don’t mind.

I’m trying to automatically calculate the Jacobian and use `:KLU` to speed up the solution:

``````function model!(res, du, u, p, t)
...
end

function solve_system(n::Int64, tspan = (0.0, 1.5), params = Params(n); jac = nothing)
ICs = ones(2 * n + 2)
ICs[end - 1] = 0.25
ICs[end] = params.r1
dICs = ones(2 * n + 2)

differential_vars = vcat(ones(n), zeros(n + 1), [1])

dae_f = DAEFunction(model!; jac = jac, jac_prototype = build_sparse(n))

prob = DAEProblem(dae_f, dICs, ICs, tspan, params, differential_vars = differential_vars)

@time sol = solve(prob, IDA(linear_solver = :KLU), reltol = 1e-10, abstol = 1e-6)

return sol
end
``````

I’m calculating the Jacobian using MTK by

``````function build_jac(n)
@variables u[1:2n+2] du[1:2n+2] γ

res = similar(u)
t = 0.0

model!(res, du, u, Params(n), t)

# need to combine d(res)/du + γ.*d(res)/d(du)
J_u  = ModelingToolkit.jacobian(res, u)
J_du = γ .* ModelingToolkit.jacobian(res, du)
J = similar(J_u)

for i in 1:length(J_u)
J[i] = simplify(J_u[i] + J_du[i])
end

J_func! = eval(build_function(J, u, γ)[2])

return (J, du, u, p, gamma, t) -> J_func!(J, u, gamma)
end
``````

This works, but `solve` seems to be recompiled for each `n` value (because the type of the Jacobian has changed?), which makes scaling up the problem impractical:

``````n = 5
@time jac = build_jac(n)
#  0.328283 seconds (676.99 k allocations: 35.423 MiB, 88.17% compilation time)
sol = solve_system(n; jac = jac);
#  1.610014 seconds (2.93 M allocations: 170.309 MiB, 3.11% gc time, 99.86% compilation time)

n = 21
@time jac = build_jac(n)
#  0.198586 seconds (598.08 k allocations: 22.969 MiB)
sol1 = solve_system(n; jac = jac);
#  8.773903 seconds (7.35 M allocations: 325.680 MiB, 1.54% gc time, 99.95% compilation time)
sol2 = solve_system(n; jac = jac);
#  0.002189 seconds (1.58 k allocations: 212.414 KiB)

julia> typeof(build_jac(n))
var"#508#511"{var"#516#517"}

julia> typeof(build_jac(n))
var"#508#511"{var"#518#519"}
``````

Do you know of a way I could fix this? Or how I would restructure the code since `n` is a parameter?

My other question is that computing the Jacobian is slower than I would expect (and if it gets too big I get `ERROR: syntax: expression too large`). Since the sparsity pattern is fairly straightforward, and I have a function to give the sparsity pattern, is there a way to tell `ModelingToolkit.jacobian` the sparsity pattern to be more efficient? I’ve also found `SparseDiffTools.jl`, but I’m not sure how to use it given the structure of `model!`.

Thanks

This needs array symbolics which is coming soon. Otherwise the structure is dependent on `n`.

https://mtk.sciml.ai/dev/tutorials/ode_modeling/#Symbolic-and-sparse-derivatives-1

`jac=true,sparse=true` creates the sparse form, and the sparsity detection is so quick that giving the sparsity pattern wouldn’t help. But if you don’t do `sparse=true`, it’ll create the dense Jacobian in the naïve way which won’t scale well.

If you give `jac_prototype=A` where `A` is the right sparse matrix structure and type, it’ll perform the coloring and all of that and make use of it in the background. At least not with Sundials. With Sundials, you’d have to do that manually. Just follow the first example and use an anonymous function to enclose the other values.

2 Likes

Thanks

Does this work with DAEs and Sundials? If I do

``````	prob = DAEProblem(model!, dICs, ICs, tspan, params, differential_vars = differential_vars; jac = true, sparse = true)
@time sol = solve(prob, IDA(linear_solver = :KLU), reltol = 1e-10, abstol = 1e-6)
``````

I get `ERROR: LoadError: MethodError: no method matching nonzeros(::Nothing)`; it looks like the Jacobian and/or the sparsity aren’t actually being added to the problem and used. If I use the default solver, adding `jac = true, sparse = true` doesn’t change the timing, or allocations.

It looks like `model!` is a function, not an `ODESystem`. All of the documentation for ModelingToolkit shows this on `ODESystem`.