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

Follow the MTK tutorials.

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.