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?