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

**URL:** <https://discourse.julialang.org/t/rodas5-and-modelingtoolkit-much-slower-than-sundials-ida-on-large-systems/60712>\
**Category:** Modelling & Simulations\
**Tags:** modelingtoolkit\
**Created:** [May 7, 2021, 1:15pm UTC](https://discourse.julialang.org/t/rodas5-and-modelingtoolkit-much-slower-than-sundials-ida-on-large-systems/60712 "2021-05-07T13:15:33Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![bmetherall](https://avatars.discourse-cdn.com/v4/letter/b/b38774/32.png) [@bmetherall](https://discourse.julialang.org/u/bmetherall)\
**Post date:** [May 7, 2021, 1:15pm UTC](https://discourse.julialang.org/t/rodas5-and-modelingtoolkit-much-slower-than-sundials-ida-on-large-systems/60712/1 "2021-05-07T13:15:33Z")

</div>

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?

```julia
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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 7, 2021, 1:54pm UTC](https://discourse.julialang.org/t/rodas5-and-modelingtoolkit-much-slower-than-sundials-ida-on-large-systems/60712/2 "2021-05-07T13:54:02Z")

</div>

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.

> [@bmetherall](#):
>
> 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.

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

---

<div class="post-metadata">

**Author:** ![bmetherall](https://avatars.discourse-cdn.com/v4/letter/b/b38774/32.png) [@bmetherall](https://discourse.julialang.org/u/bmetherall)\
**Post date:** [May 7, 2021, 3:42pm UTC](https://discourse.julialang.org/t/rodas5-and-modelingtoolkit-much-slower-than-sundials-ida-on-large-systems/60712/3 "2021-05-07T15:42:42Z")

</div>

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

```julia
	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**
>
> ```julia
> 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
> in expression starting at /home/brady/Desktop/Silicon_Processes/Numerics/Example.jl:78
> 
> ```

If not, would following something like  
[https://discourse.julialang.org/t/error-evaluating-sparse-analytical-jacobian-from-modelingtoolkit-in-dae-system/45764/2](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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 9, 2021, 10:51am UTC](https://discourse.julialang.org/t/rodas5-and-modelingtoolkit-much-slower-than-sundials-ida-on-large-systems/60712/4 "2021-05-09T10:51:57Z")

</div>

> [@bmetherall](#):
>
> 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

Just do it the easy way.

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

```

---

<div class="post-metadata">

**Author:** ![bmetherall](https://avatars.discourse-cdn.com/v4/letter/b/b38774/32.png) [@bmetherall](https://discourse.julialang.org/u/bmetherall)\
**Post date:** [May 10, 2021, 5:27pm UTC](https://discourse.julialang.org/t/rodas5-and-modelingtoolkit-much-slower-than-sundials-ida-on-large-systems/60712/5 "2021-05-10T17:27:20Z")

</div>

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:

```julia
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

```julia
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:

```julia
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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 10, 2021, 5:35pm UTC](https://discourse.julialang.org/t/rodas5-and-modelingtoolkit-much-slower-than-sundials-ida-on-large-systems/60712/6 "2021-05-10T17:35:03Z")

</div>

> [@bmetherall](#):
>
> Do you know of a way I could fix this? Or how I would restructure the code since `n` is a parameter?

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

> [@bmetherall](#):
>
> 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?

Follow the MTK tutorials.

[https://mtk.sciml.ai/dev/tutorials/ode\_modeling/#Symbolic-and-sparse-derivatives-1](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.

> [@bmetherall](#):
>
> I’ve also found `SparseDiffTools.jl` , but I’m not sure how to use it given the structure of `model!` .

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.

---

<div class="post-metadata">

**Author:** ![bmetherall](https://avatars.discourse-cdn.com/v4/letter/b/b38774/32.png) [@bmetherall](https://discourse.julialang.org/u/bmetherall)\
**Post date:** [May 10, 2021, 5:56pm UTC](https://discourse.julialang.org/t/rodas5-and-modelingtoolkit-much-slower-than-sundials-ida-on-large-systems/60712/7 "2021-05-10T17:56:24Z")

</div>

Thanks

> [@ChrisRackauckas](#):
>
> Follow the MTK tutorials.
> 
> [https://mtk.sciml.ai/dev/tutorials/ode\_modeling/#Symbolic-and-sparse-derivatives-1](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.

Does this work with DAEs and Sundials? If I do

```julia
	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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 13, 2021, 1:21am UTC](https://discourse.julialang.org/t/rodas5-and-modelingtoolkit-much-slower-than-sundials-ida-on-large-systems/60712/8 "2021-05-13T01:21:38Z")

</div>

> [@bmetherall](#):
>
> 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`.