# 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:** 1
**Showing post:** 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?

---

_[View the full topic](https://discourse.julialang.org/t/rodas5-and-modelingtoolkit-much-slower-than-sundials-ida-on-large-systems/60712)._