Brusselator Scaling PDE Differentiation Benchmarks

From the paper A Comparison of Automatic Differentiation and Continuous Sensitivity Analysis for Derivatives of Differential Equation Solutions

using OrdinaryDiffEq, ReverseDiff, ForwardDiff, FiniteDiff, SciMLSensitivity
using OrdinaryDiffEqRosenbrock
using LinearAlgebra, Tracker, Mooncake, Plots
function makebrusselator(N = 8)
    xyd_brusselator = range(0, stop = 1, length = N)
    function limit(a, N)
        if a == N+1
            return 1
        elseif a == 0
            return N
        else
            return a
        end
    end
    brusselator_f(x, y, t) = ifelse(
        (((x-0.3)^2 + (y-0.6)^2) <= 0.1^2) &&
        (t >= 1.1), 5.0, 0.0)
    brusselator_2d_loop = let N=N, xyd=xyd_brusselator, dx=step(xyd_brusselator)
        function brusselator_2d_loop(du, u, p, t)
            @inbounds begin
                ii1 = N^2
                ii2 = ii1+N^2
                ii3 = ii2+2(N^2)
                A = @view p[1:ii1]
                B = @view p[(ii1 + 1):ii2]
                α = @view p[(ii2 + 1):ii3]
                II = LinearIndices((N, N, 2))
                for I in CartesianIndices((N, N))
                    x = xyd[I[1]]
                    y = xyd[I[2]]
                    i = I[1]
                    j = I[2]
                    ip1 = limit(i+1, N);
                    im1 = limit(i-1, N)
                    jp1 = limit(j+1, N);
                    jm1 = limit(j-1, N)
                    du[II[i, j, 1]] = α[II[
                                          i, j, 1]]*(u[II[im1, j, 1]] + u[II[ip1, j, 1]] +
                                                     u[II[i, jp1, 1]] + u[II[i, jm1, 1]] -
                                                     4u[II[i, j, 1]])/dx^2 +
                                      B[II[i, j, 1]] + u[II[i, j, 1]]^2*u[II[i, j, 2]] -
                                      (A[II[i, j, 1]] + 1)*u[II[i, j, 1]] +
                                      brusselator_f(x, y, t)
                end
                for I in CartesianIndices((N, N))
                    i = I[1]
                    j = I[2]
                    ip1 = limit(i+1, N)
                    im1 = limit(i-1, N)
                    jp1 = limit(j+1, N)
                    jm1 = limit(j-1, N)
                    du[II[i, j, 2]] = α[II[
                        i, j, 2]]*(u[II[im1, j, 2]] + u[II[ip1, j, 2]] + u[II[i, jp1, 2]] +
                                   u[II[i, jm1, 2]] - 4u[II[i, j, 2]])/dx^2 +
                                      A[II[i, j, 1]]*u[II[i, j, 1]] -
                                      u[II[i, j, 1]]^2*u[II[i, j, 2]]
                end
                return nothing
            end
        end
    end
    function init_brusselator_2d(xyd)
        N = length(xyd)
        u = zeros(N, N, 2)
        for I in CartesianIndices((N, N))
            x = xyd[I[1]]
            y = xyd[I[2]]
            u[I, 1] = 22*(y*(1-y))^(3/2)
            u[I, 2] = 27*(x*(1-x))^(3/2)
        end
        vec(u)
    end
    dx = step(xyd_brusselator)
    e1 = ones(N-1)
    off = N-1
    e4 = ones(N-off)
    T = diagm(0=>-2ones(N), -1=>e1, 1=>e1, off=>e4, -off=>e4) ./ dx^2
    Ie = Matrix{Float64}(I, N, N)
    # A + df/du
    Op = kron(Ie, T) + kron(T, Ie)
    brusselator_jac = let N=N
        (J, a, p, t) -> begin
            ii1 = N^2
            ii2 = ii1+N^2
            ii3 = ii2+2(N^2)
            A = @view p[1:ii1]
            B = @view p[(ii1 + 1):ii2]
            α = @view p[(ii2 + 1):ii3]
            u = @view a[1:(end ÷ 2)]
            v = @view a[(end ÷ 2 + 1):end]
            N2 = length(a)÷2
            α1 = @view α[1:(end ÷ 2)]
            α2 = @view α[(end ÷ 2 + 1):end]
            fill!(J, 0)

            J[1:N2, 1:N2] .= α1 .* Op
            J[(N2 + 1):end, (N2 + 1):end] .= α2 .* Op

            J1 = @view J[1:N2, 1:N2]
            J2 = @view J[(N2 + 1):end, 1:N2]
            J3 = @view J[1:N2, (N2 + 1):end]
            J4 = @view J[(N2 + 1):end, (N2 + 1):end]
            J1[diagind(J1)] .+= @. 2u*v-(A+1)
            J2[diagind(J2)] .= @. A-2u*v
            J3[diagind(J3)] .= @. u^2
            J4[diagind(J4)] .+= @. -u^2
            nothing
        end
    end
    Jmat = zeros(2N*N, 2N*N)
    dp = zeros(2N*N, 4N*N)
    brusselator_comp = let N=N, xyd=xyd_brusselator, dx=step(xyd_brusselator), Jmat=Jmat,
        dp=dp, brusselator_jac=brusselator_jac

        function brusselator_comp(dus, us, p, t)
            @inbounds begin
                ii1 = N^2
                ii2 = ii1+N^2
                ii3 = ii2+2(N^2)
                @views u, s = us[1:ii2], us[(ii2 + 1):end]
                du = @view dus[1:ii2]
                ds = @view dus[(ii2 + 1):end]
                fill!(dp, 0)
                A = @view p[1:ii1]
                B = @view p[(ii1 + 1):ii2]
                α = @view p[(ii2 + 1):ii3]
                dfdα = @view dp[:, (ii2 + 1):ii3]
                diagind(dfdα)
                for i in 1:ii1
                    dp[i, ii1 + i] = 1
                end
                II = LinearIndices((N, N, 2))
                uu = @view u[1:(end ÷ 2)]
                for i in eachindex(uu)
                    dp[i, i] = -uu[i]
                    dp[i + ii1, i] = uu[i]
                end
                for I in CartesianIndices((N, N))
                    x = xyd[I[1]]
                    y = xyd[I[2]]
                    i = I[1]
                    j = I[2]
                    ip1 = limit(i+1, N);
                    im1 = limit(i-1, N)
                    jp1 = limit(j+1, N);
                    jm1 = limit(j-1, N)
                    au = dfdα[II[i, j, 1], II[i, j, 1]] = (u[II[im1, j, 1]] +
                                                           u[II[ip1, j, 1]] +
                                                           u[II[i, jp1, 1]] +
                                                           u[II[i, jm1, 1]] -
                                                           4u[II[i, j, 1]])/dx^2
                    du[II[i, j, 1]] = α[II[i, j, 1]]*(au) + B[II[i, j, 1]] +
                                      u[II[i, j, 1]]^2*u[II[i, j, 2]] -
                                      (A[II[i, j, 1]] + 1)*u[II[i, j, 1]] +
                                      brusselator_f(x, y, t)
                end
                for I in CartesianIndices((N, N))
                    i = I[1]
                    j = I[2]
                    ip1 = limit(i+1, N)
                    im1 = limit(i-1, N)
                    jp1 = limit(j+1, N)
                    jm1 = limit(j-1, N)
                    av = dfdα[II[i, j, 2], II[i, j, 2]] = (u[II[im1, j, 2]] +
                                                           u[II[ip1, j, 2]] +
                                                           u[II[i, jp1, 2]] +
                                                           u[II[i, jm1, 2]] -
                                                           4u[II[i, j, 2]])/dx^2
                    du[II[i, j, 2]] = α[II[i, j, 2]]*(av) + A[II[i, j, 1]]*u[II[i, j, 1]] -
                                      u[II[i, j, 1]]^2*u[II[i, j, 2]]
                end
                brusselator_jac(Jmat, u, p, t)
                BLAS.gemm!('N', 'N', 1.0, Jmat, reshape(s, 2N*N, 4N*N), 1.0, dp)
                copyto!(ds, vec(dp))
                return nothing
            end
        end
    end
    u0 = init_brusselator_2d(xyd_brusselator)
    p = [fill(3.4, N^2); fill(1.0, N^2); fill(10.0, 2*N^2)]
    brusselator_2d_loop, u0,
    p,
    brusselator_jac,
    ODEProblem(brusselator_comp, copy([u0; zeros((N^2*2)*(N^2*4))]), (0.0, 10.0), p)
end

Base.eps(::Type{Tracker.TrackedReal{T}}) where {T} = eps(T)
Base.vec(v::Adjoint{<:Real, <:AbstractVector}) = vec(v') # bad bad hack

Setup AutoDiff

bt = 0:0.1:1
tspan = (0.0, 1.0)
forwarddiffn = vcat(2:10, 12, 15)
reversediffn = 2:10
numdiffn = vcat(2:10, 12)
csan = vcat(2:10, 12, 15, 17)
#csaseedn = 2:10
tols = (abstol = 1e-5, reltol = 1e-7)

@isdefined(PROBS) || (const PROBS = Dict{Int, Any}())
makebrusselator!(dict, n) = get!(()->makebrusselator(n), dict, n)

_adjoint_methods_iq = ntuple(2) do ii
    Alg = (InterpolatingAdjoint, QuadratureAdjoint)[ii]
    (
        user = Alg(autodiff = false, autojacvec = false), # user Jacobian
        adjc = Alg(autodiff = true, autojacvec = false), # AD Jacobian
        advj = Alg(autodiff = true, autojacvec = EnzymeVJP()) # AD vJ
    )
end |> NamedTuple{(:interp, :quad)}
# GaussAdjoint/GaussKronrodAdjoint do not support user-provided Jacobians (autodiff=false)
_adjoint_methods_g = ntuple(2) do ii
    Alg = (GaussAdjoint, GaussKronrodAdjoint)[ii]
    (
        adjc = Alg(autodiff = true, autojacvec = false), # AD Jacobian
        advj = Alg(autodiff = true, autojacvec = EnzymeVJP()) # AD vJ
    )
end |> NamedTuple{(:gauss, :gausskronrod)}
@isdefined(ADJOINT_METHODS_IQ) ||
    (const ADJOINT_METHODS_IQ = mapreduce(collect, vcat, _adjoint_methods_iq))
@isdefined(ADJOINT_METHODS_G) ||
    (const ADJOINT_METHODS_G = mapreduce(collect, vcat, _adjoint_methods_g))

function auto_sen_l2(
        f, u0, tspan, p, t, alg = Tsit5(); diffalg = ReverseDiff.gradient, kwargs...)
    test_f(p) = begin
        prob = ODEProblem{true, SciMLBase.FullSpecialize}(f, convert.(eltype(p), u0), tspan, p)
        sol = solve(prob, alg, saveat = t; kwargs...)
        sum(sol.u) do x
            sum(z->(1-z)^2/2, x)
        end
    end
    diffalg(test_f, p)
end
@inline function diffeq_sen_l2(df, u0, tspan, p, t, alg = Tsit5();
        abstol = 1e-5, reltol = 1e-7, iabstol = abstol, ireltol = reltol,
        sensalg = SensitivityAlg(), kwargs...)
    prob = ODEProblem{true, SciMLBase.FullSpecialize}(df, u0, tspan, p)
    saveat = tspan[1] != t[1] && tspan[end] != t[end] ? vcat(tspan[1], t, tspan[end]) : t
    sol = solve(prob, alg, abstol = abstol, reltol = reltol, saveat = saveat; kwargs...)
    dg(out, u, p, t, i) = (out.=u .- 1.0)
    adjoint_sensitivities(sol, alg; t, abstol = abstol, dgdu_discrete = dg,
        reltol = reltol, sensealg = sensalg)
end
diffeq_sen_l2 (generic function with 2 methods)

AD Choice Benchmarks

forwarddiff = map(forwarddiffn) do n
    bfun, b_u0, b_p, brusselator_jac, brusselator_comp = makebrusselator!(PROBS, n)
    @elapsed auto_sen_l2(
        bfun, b_u0, tspan, b_p, bt, (Rodas5()); diffalg = (ForwardDiff.gradient), tols...)
    t = @elapsed auto_sen_l2(
        bfun, b_u0, tspan, b_p, bt, (Rodas5()); diffalg = (ForwardDiff.gradient), tols...)
    @show n, t
    t
end
(n, t) = (2, 0.001156163)
(n, t) = (3, 0.009697823)
(n, t) = (4, 0.038978882)
(n, t) = (5, 0.142719323)
(n, t) = (6, 0.403370744)
(n, t) = (7, 0.961728277)
(n, t) = (8, 2.234831128)
(n, t) = (9, 4.052730812)
(n, t) = (10, 9.359019992)
(n, t) = (12, 30.720095731)
(n, t) = (15, 120.941890062)
11-element Vector{Float64}:
   0.001156163
   0.009697823
   0.038978882
   0.142719323
   0.403370744
   0.961728277
   2.234831128
   4.052730812
   9.359019992
  30.720095731
 120.941890062
#=
reversediff = map(reversediffn) do n
  bfun, b_u0, b_p, brusselator_jac, brusselator_comp = makebrusselator!(PROBS, n)
  @elapsed auto_sen_l2(bfun, b_u0, tspan, b_p, bt, (Rodas5(autodiff=AutoFiniteDiff())); diffalg=(ReverseDiff.gradient), tols...)
  t = @elapsed auto_sen_l2(bfun, b_u0, tspan, b_p, bt, (Rodas5(autodiff=AutoFiniteDiff())); diffalg=(ReverseDiff.gradient), tols...)
  @show n,t
  t
end
=#
numdiff = map(numdiffn) do n
    bfun, b_u0, b_p, brusselator_jac, brusselator_comp = makebrusselator!(PROBS, n)
    @elapsed auto_sen_l2(bfun, b_u0, tspan, b_p, bt, (Rodas5());
        diffalg = (FiniteDiff.finite_difference_gradient), tols...)
    t = @elapsed auto_sen_l2(bfun, b_u0, tspan, b_p, bt, (Rodas5());
        diffalg = (FiniteDiff.finite_difference_gradient), tols...)
    @show n, t
    t
end
(n, t) = (2, 0.003685878)
(n, t) = (3, 0.030100603)
(n, t) = (4, 0.101166087)
(n, t) = (5, 0.320808188)
(n, t) = (6, 0.848674562)
(n, t) = (7, 2.193527103)
(n, t) = (8, 4.37603293)
(n, t) = (9, 10.292124476)
(n, t) = (10, 18.450115912)
(n, t) = (12, 89.21977739)
10-element Vector{Float64}:
  0.003685878
  0.030100603
  0.101166087
  0.320808188
  0.848674562
  2.193527103
  4.37603293
 10.292124476
 18.450115912
 89.21977739

Warmup: run each adjoint method once at the smallest size to ensure all compilation is complete before we start timing.

let n = first(csan)
    bfun, b_u0, b_p, brusselator_jac, brusselator_comp = makebrusselator!(PROBS, n)
    solver = Rodas5(autodiff = AutoFiniteDiff())
    for alg in ADJOINT_METHODS_IQ
        f = SciMLSensitivity.alg_autodiff(alg) ? bfun :
            ODEFunction(bfun, jac = brusselator_jac)
        diffeq_sen_l2(f, b_u0, tspan, b_p, bt, solver; sensalg = alg, tols...)
    end
    for alg in ADJOINT_METHODS_G
        diffeq_sen_l2(bfun, b_u0, tspan, b_p, bt, solver; sensalg = alg, tols...)
    end
end
csa_iq = map(csan) do n
    bfun, b_u0, b_p, brusselator_jac, brusselator_comp = makebrusselator!(PROBS, n)
    @time ts = map(ADJOINT_METHODS_IQ) do alg
        @info "Running $alg"
        f = SciMLSensitivity.alg_autodiff(alg) ? bfun :
            ODEFunction(bfun, jac = brusselator_jac)
        solver = Rodas5(autodiff = AutoFiniteDiff())
        @time diffeq_sen_l2(f, b_u0, tspan, b_p, bt, solver; sensalg = alg, tols...)
        t = @elapsed diffeq_sen_l2(f, b_u0, tspan, b_p, bt, solver; sensalg = alg, tols...)
        return t
    end
    @show n, ts
    ts
end
0.006164 seconds (9.61 k allocations: 1.117 MiB)
  0.003583 seconds (6.32 k allocations: 2.050 MiB)
  0.003578 seconds (7.34 k allocations: 569.445 KiB)
  0.001953 seconds (4.58 k allocations: 304.148 KiB)
  0.001717 seconds (2.30 k allocations: 230.867 KiB)
  0.002788 seconds (7.90 k allocations: 499.922 KiB)
  0.976041 seconds (1.21 M allocations: 69.317 MiB, 94.29% compilation time
: 5% of which was recompilation)
(n, ts) = (2, [0.006010405, 0.002860783, 0.003242051, 0.00168653, 0.0013033
42, 0.002396186])
  0.025931 seconds (14.74 k allocations: 2.059 MiB)
 10.131532 seconds (8.09 M allocations: 388.735 MiB, 0.75% gc time, 99.86% 
compilation time)
  0.010108 seconds (10.96 k allocations: 994.891 KiB)
  0.003899 seconds (5.60 k allocations: 550.633 KiB)
 10.517351 seconds (6.85 M allocations: 326.931 MiB, 0.91% gc time, 99.94% 
compilation time)
  0.005808 seconds (9.10 k allocations: 695.078 KiB)
 20.766679 seconds (15.03 M allocations: 726.160 MiB, 0.83% gc time, 99.33%
 compilation time)
(n, ts) = (3, [0.024999653, 0.012127769, 0.009983322, 0.003728068, 0.003833
067, 0.005496518])
  0.107499 seconds (22.87 k allocations: 4.004 MiB)
 10.226935 seconds (8.10 M allocations: 389.627 MiB, 2.59% gc time, 99.54% 
compilation time)
  0.021784 seconds (16.77 k allocations: 1.758 MiB)
  0.008954 seconds (7.44 k allocations: 1.057 MiB)
 10.370348 seconds (6.85 M allocations: 327.160 MiB, 0.45% gc time, 99.88% 
compilation time)
  0.010933 seconds (11.39 k allocations: 1005.797 KiB)
 20.960383 seconds (15.08 M allocations: 735.293 MiB, 1.48% gc time, 97.98%
 compilation time)
(n, ts) = (4, [0.106446596, 0.043759773, 0.021849043, 0.008914557, 0.010683
967, 0.011065245])
  0.360743 seconds (32.72 k allocations: 6.980 MiB)
  9.586964 seconds (7.12 M allocations: 343.371 MiB, 0.92% gc time, 98.18% 
compilation time)
  0.047880 seconds (23.80 k allocations: 3.015 MiB)
  0.020269 seconds (9.86 k allocations: 1.938 MiB)
 10.548319 seconds (6.84 M allocations: 327.735 MiB, 0.87% gc time, 99.66% 
compilation time)
  0.019507 seconds (14.40 k allocations: 1.447 MiB)
 21.249844 seconds (14.14 M allocations: 702.312 MiB, 0.85% gc time, 93.77%
 compilation time)
(n, ts) = (5, [0.35972572, 0.17215279, 0.048393135, 0.020228451, 0.03381940
2, 0.020005683])
  0.994032 seconds (44.78 k allocations: 12.106 MiB)
  9.967069 seconds (7.13 M allocations: 345.994 MiB, 0.72% gc time, 95.05% 
compilation time)
  0.091176 seconds (32.42 k allocations: 5.000 MiB)
  0.040011 seconds (12.60 k allocations: 3.539 MiB)
 10.541472 seconds (6.85 M allocations: 328.008 MiB, 1.48% gc time, 99.33% 
compilation time)
  0.034958 seconds (17.82 k allocations: 2.003 MiB)
 23.489257 seconds (14.22 M allocations: 726.501 MiB, 1.37% gc time, 84.91%
 compilation time)
(n, ts) = (6, [0.990946328, 0.490307554, 0.182655568, 0.03917235, 0.0685727
08, 0.035033565])
  2.541370 seconds (61.29 k allocations: 19.863 MiB)
  9.927483 seconds (6.05 M allocations: 297.280 MiB, 0.35% gc time, 89.04% 
compilation time)
  0.173765 seconds (44.24 k allocations: 7.965 MiB)
  0.074872 seconds (15.82 k allocations: 5.740 MiB)
  2.273085 seconds (789.87 k allocations: 41.046 MiB, 92.87% compilation ti
me)
  0.054741 seconds (21.85 k allocations: 2.693 MiB)
 19.285144 seconds (7.16 M allocations: 421.874 MiB, 0.85% gc time, 56.78% 
compilation time)
(n, ts) = (7, [2.594097966, 1.12190019, 0.17391593, 0.075661206, 0.20383925
5, 0.055536354])
  5.736195 seconds (78.33 k allocations: 30.015 MiB)
  2.764484 seconds (24.18 k allocations: 11.810 MiB)
  0.296238 seconds (56.40 k allocations: 12.095 MiB)
  0.126767 seconds (20.61 k allocations: 9.053 MiB)
  0.344395 seconds (7.27 k allocations: 4.423 MiB)
  0.085445 seconds (27.71 k allocations: 3.636 MiB)
 18.647406 seconds (431.11 k allocations: 142.931 MiB, 0.61% gc time)
(n, ts) = (8, [5.58743157, 2.811020104, 0.297288077, 0.155598257, 0.3433718
76, 0.084827333])
 12.642975 seconds (111.52 k allocations: 48.410 MiB, 0.44% gc time)
  6.568246 seconds (34.65 k allocations: 18.375 MiB, 0.43% gc time)
  0.771444 seconds (107.37 k allocations: 20.228 MiB, 6.43% gc time)
  0.260190 seconds (24.97 k allocations: 13.632 MiB, 14.26% gc time)
  0.677881 seconds (8.36 k allocations: 6.499 MiB)
  0.138889 seconds (33.16 k allocations: 4.793 MiB)
 42.277904 seconds (642.17 k allocations: 224.740 MiB, 0.88% gc time)
(n, ts) = (9, [12.687357132, 6.569018145, 0.722367394, 0.222826403, 0.85754
295, 0.13817358])
 21.380205 seconds (123.39 k allocations: 67.505 MiB, 0.78% gc time)
 10.836489 seconds (37.06 k allocations: 26.776 MiB)
  0.839252 seconds (88.59 k allocations: 25.620 MiB, 2.10% gc time)
  0.351589 seconds (29.83 k allocations: 20.838 MiB)
  1.245915 seconds (9.58 k allocations: 9.225 MiB)
  0.194954 seconds (39.23 k allocations: 6.158 MiB)
 69.721582 seconds (657.46 k allocations: 313.108 MiB, 0.58% gc time)
(n, ts) = (10, [21.26981936, 10.87202168, 0.863249517, 0.381403743, 1.26805
7573, 0.195809412])
 64.533957 seconds (181.92 k allocations: 130.713 MiB, 0.81% gc time)
 32.898040 seconds (54.28 k allocations: 53.241 MiB, 0.11% gc time)
  2.364259 seconds (130.62 k allocations: 48.696 MiB, 11.30% gc time)
  0.958237 seconds (41.39 k allocations: 40.368 MiB, 7.58% gc time)
  3.966074 seconds (12.68 k allocations: 17.459 MiB, 0.33% gc time)
  0.446718 seconds (53.61 k allocations: 9.886 MiB, 3.46% gc time)
209.583028 seconds (951.13 k allocations: 601.596 MiB, 0.59% gc time)
(n, ts) = (12, [64.11596754, 32.81239846, 2.22189177, 0.889011557, 3.920252
231, 0.429867029])
232.368557 seconds (271.09 k allocations: 298.637 MiB, 0.18% gc time)
124.448534 seconds (82.50 k allocations: 125.349 MiB, 0.04% gc time)
  5.357463 seconds (207.92 k allocations: 109.802 MiB, 5.23% gc time)
  3.042593 seconds (62.13 k allocations: 96.879 MiB, 8.60% gc time)
 14.956209 seconds (17.87 k allocations: 39.433 MiB, 0.10% gc time)
  1.066714 seconds (79.53 k allocations: 18.623 MiB)
761.084272 seconds (1.44 M allocations: 1.346 GiB, 0.30% gc time)
(n, ts) = (15, [230.927492573, 124.882455901, 5.251848781, 2.7871184, 14.91
2401246, 1.055379152])
585.457308 seconds (475.80 k allocations: 479.871 MiB, 0.13% gc time)
641.719892 seconds (101.29 k allocations: 203.646 MiB, 0.09% gc time)
  8.445158 seconds (248.09 k allocations: 173.692 MiB, 3.59% gc time)
  5.812806 seconds (116.40 k allocations: 149.973 MiB, 3.78% gc time)
 32.346769 seconds (21.97 k allocations: 63.004 MiB, 0.04% gc time)
  1.824952 seconds (100.02 k allocations: 27.201 MiB, 0.94% gc time)
2550.175775 seconds (2.13 M allocations: 2.144 GiB, 0.14% gc time)
(n, ts) = (17, [586.599723097, 639.725617118, 8.260552394, 5.771834001, 32.
290931626, 1.893561637])
12-element Vector{Vector{Float64}}:
 [0.006010405, 0.002860783, 0.003242051, 0.00168653, 0.001303342, 0.0023961
86]
 [0.024999653, 0.012127769, 0.009983322, 0.003728068, 0.003833067, 0.005496
518]
 [0.106446596, 0.043759773, 0.021849043, 0.008914557, 0.010683967, 0.011065
245]
 [0.35972572, 0.17215279, 0.048393135, 0.020228451, 0.033819402, 0.02000568
3]
 [0.990946328, 0.490307554, 0.182655568, 0.03917235, 0.068572708, 0.0350335
65]
 [2.594097966, 1.12190019, 0.17391593, 0.075661206, 0.203839255, 0.05553635
4]
 [5.58743157, 2.811020104, 0.297288077, 0.155598257, 0.343371876, 0.0848273
33]
 [12.687357132, 6.569018145, 0.722367394, 0.222826403, 0.85754295, 0.138173
58]
 [21.26981936, 10.87202168, 0.863249517, 0.381403743, 1.268057573, 0.195809
412]
 [64.11596754, 32.81239846, 2.22189177, 0.889011557, 3.920252231, 0.4298670
29]
 [230.927492573, 124.882455901, 5.251848781, 2.7871184, 14.912401246, 1.055
379152]
 [586.599723097, 639.725617118, 8.260552394, 5.771834001, 32.290931626, 1.8
93561637]
csa_g = map(csan) do n
    bfun, b_u0, b_p, brusselator_jac, brusselator_comp = makebrusselator!(PROBS, n)
    @time ts = map(ADJOINT_METHODS_G) do alg
        @info "Running $alg"
        solver = Rodas5(autodiff = AutoFiniteDiff())
        @time diffeq_sen_l2(bfun, b_u0, tspan, b_p, bt, solver; sensalg = alg, tols...)
        t = @elapsed diffeq_sen_l2(bfun, b_u0, tspan, b_p, bt, solver; sensalg = alg, tols...)
        return t
    end
    @show n, ts
    ts
end
0.001977 seconds (2.54 k allocations: 289.445 KiB)
  0.002873 seconds (6.17 k allocations: 488.875 KiB)
  0.002164 seconds (4.75 k allocations: 338.898 KiB)
  0.003212 seconds (8.92 k allocations: 575.453 KiB)
  0.493268 seconds (443.93 k allocations: 24.436 MiB, 93.99% compilation ti
me)
(n, ts) = (2, [0.001262843, 0.002455946, 0.001500031, 0.002805444])
  9.132297 seconds (7.52 M allocations: 372.286 MiB, 0.95% gc time, 99.92% 
compilation time)
  0.005720 seconds (7.19 k allocations: 682.188 KiB)
  9.312185 seconds (7.73 M allocations: 382.866 MiB, 0.97% gc time, 99.91% 
compilation time)
  0.007185 seconds (13.23 k allocations: 910.062 KiB)
 18.487747 seconds (15.31 M allocations: 759.952 MiB, 0.96% gc time, 99.68%
 compilation time)
(n, ts) = (3, [0.003861318, 0.005553688, 0.0053716, 0.007281889])
  9.178591 seconds (7.52 M allocations: 372.576 MiB, 1.03% gc time, 99.86% 
compilation time)
  0.010425 seconds (9.26 k allocations: 1.067 MiB)
  9.403251 seconds (7.74 M allocations: 383.305 MiB, 1.88% gc time, 99.80% 
compilation time)
  0.012790 seconds (18.10 k allocations: 1.426 MiB)
 18.660047 seconds (15.33 M allocations: 763.448 MiB, 1.45% gc time, 99.41%
 compilation time)
(n, ts) = (4, [0.010079942, 0.010244431, 0.015508861, 0.012650928])
  8.597912 seconds (5.83 M allocations: 279.212 MiB, 1.91% gc time, 99.60% 
compilation time)
  0.020216 seconds (11.74 k allocations: 1.742 MiB)
  8.632248 seconds (6.20 M allocations: 297.363 MiB, 0.93% gc time, 99.47% 
compilation time)
  0.023687 seconds (20.58 k allocations: 2.155 MiB)
 17.400487 seconds (12.10 M allocations: 588.371 MiB, 1.41% gc time, 98.56%
 compilation time)
(n, ts) = (5, [0.031942476, 0.020038735, 0.043054074, 0.023830813])
  8.539620 seconds (5.83 M allocations: 280.505 MiB, 1.06% gc time, 99.19% 
compilation time)
  0.035030 seconds (14.75 k allocations: 2.842 MiB)
  8.772557 seconds (6.20 M allocations: 298.834 MiB, 1.07% gc time, 98.84% 
compilation time)
  0.041140 seconds (27.09 k allocations: 3.487 MiB)
 17.638690 seconds (12.13 M allocations: 598.852 MiB, 1.04% gc time, 97.18%
 compilation time)
(n, ts) = (6, [0.067195425, 0.03478832, 0.09951842, 0.041750331])
  8.516429 seconds (5.64 M allocations: 273.274 MiB, 1.19% gc time, 97.71% 
compilation time)
  0.057910 seconds (18.22 k allocations: 4.432 MiB)
  8.739633 seconds (6.01 M allocations: 291.756 MiB, 0.90% gc time, 97.00% 
compilation time)
  0.065746 seconds (32.50 k allocations: 5.265 MiB)
 17.918657 seconds (11.78 M allocations: 595.184 MiB, 1.24% gc time, 93.75%
 compilation time)
(n, ts) = (7, [0.144045944, 0.100390684, 0.219824239, 0.06633376])
  0.311812 seconds (6.60 k allocations: 7.829 MiB)
  0.140833 seconds (22.29 k allocations: 6.766 MiB, 37.00% gc time)
  0.455795 seconds (19.32 k allocations: 8.627 MiB)
  0.155720 seconds (40.35 k allocations: 7.927 MiB, 35.65% gc time)
  2.025740 seconds (178.48 k allocations: 62.876 MiB, 5.31% gc time)
(n, ts) = (8, [0.311535643, 0.088565152, 0.452656584, 0.099789767])
  0.610036 seconds (7.51 k allocations: 12.114 MiB)
  0.146916 seconds (26.88 k allocations: 10.085 MiB)
  0.918568 seconds (22.78 k allocations: 13.252 MiB, 2.69% gc time)
  0.162601 seconds (47.70 k allocations: 11.604 MiB)
  3.719951 seconds (211.10 k allocations: 94.687 MiB, 2.23% gc time)
(n, ts) = (9, [0.64178441, 0.146224121, 0.888184607, 0.195922646])
  1.165480 seconds (8.77 k allocations: 17.895 MiB)
  0.219044 seconds (34.16 k allocations: 14.581 MiB)
  1.567767 seconds (23.96 k allocations: 19.216 MiB)
  0.235283 seconds (55.54 k allocations: 16.334 MiB)
  6.652948 seconds (246.20 k allocations: 136.629 MiB, 4.43% gc time)
(n, ts) = (10, [1.215053953, 0.384120528, 1.589022757, 0.266622931])
  4.153744 seconds (13.50 k allocations: 35.815 MiB, 0.64% gc time)
  0.500171 seconds (46.18 k allocations: 27.761 MiB, 2.48% gc time)
  5.344365 seconds (31.86 k allocations: 37.959 MiB, 0.24% gc time)
  0.713651 seconds (72.60 k allocations: 30.481 MiB, 27.79% gc time)
 21.490143 seconds (329.63 k allocations: 264.612 MiB, 2.46% gc time)
(n, ts) = (12, [4.176005334, 0.518152929, 5.372500771, 0.700359394])
 13.413497 seconds (16.55 k allocations: 84.682 MiB, 0.22% gc time)
  1.451386 seconds (73.53 k allocations: 63.578 MiB, 13.22% gc time)
 17.380219 seconds (40.63 k allocations: 88.781 MiB, 0.31% gc time)
  1.428808 seconds (104.84 k allocations: 68.045 MiB, 2.18% gc time)
 67.619343 seconds (472.47 k allocations: 610.750 MiB, 1.12% gc time)
(n, ts) = (15, [13.540071843, 1.305639031, 17.61057016, 1.477409435])
 27.731223 seconds (19.56 k allocations: 138.004 MiB, 0.48% gc time)
  2.769336 seconds (89.78 k allocations: 101.685 MiB, 8.97% gc time)
 52.268320 seconds (42.30 k allocations: 142.854 MiB, 0.21% gc time)
  2.693686 seconds (126.33 k allocations: 108.002 MiB, 9.11% gc time)
169.862872 seconds (557.30 k allocations: 981.669 MiB, 0.76% gc time)
(n, ts) = (17, [27.941985967, 2.392182254, 51.588983067, 2.465422016])
12-element Vector{Vector{Float64}}:
 [0.001262843, 0.002455946, 0.001500031, 0.002805444]
 [0.003861318, 0.005553688, 0.0053716, 0.007281889]
 [0.010079942, 0.010244431, 0.015508861, 0.012650928]
 [0.031942476, 0.020038735, 0.043054074, 0.023830813]
 [0.067195425, 0.03478832, 0.09951842, 0.041750331]
 [0.144045944, 0.100390684, 0.219824239, 0.06633376]
 [0.311535643, 0.088565152, 0.452656584, 0.099789767]
 [0.64178441, 0.146224121, 0.888184607, 0.195922646]
 [1.215053953, 0.384120528, 1.589022757, 0.266622931]
 [4.176005334, 0.518152929, 5.372500771, 0.700359394]
 [13.540071843, 1.305639031, 17.61057016, 1.477409435]
 [27.941985967, 2.392182254, 51.588983067, 2.465422016]
n_to_param(n) = 4n^2

lw = 2
ms = 0.5
plt1 = plot(title = "Sensitivity Scaling on Brusselator");
plot!(plt1, n_to_param.(forwarddiffn), forwarddiff, lab = "Forward-Mode DSAAD",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
#plot!(plt1, n_to_param.(reversediffn), reversediff, lab="Reverse-Mode DSAAD", lw=lw, marksize=ms, linestyle=:auto, marker=:auto);
csadata_iq = [[csa_iq[j][i] for j in eachindex(csa_iq)] for i in eachindex(csa_iq[1])]
csadata_g = [[csa_g[j][i] for j in eachindex(csa_g)] for i in eachindex(csa_g[1])]
plot!(plt1, n_to_param.(csan), csadata_iq[1], lab = "Interpolating CASA user-Jacobian",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(plt1, n_to_param.(csan), csadata_iq[2], lab = "Interpolating CASA AD-Jacobian",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(
    plt1, n_to_param.(csan), csadata_iq[3], lab = raw"Interpolating CASA AD-$v^{T}J$ seeding",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(plt1, n_to_param.(csan), csadata_iq[1 + 3], lab = "Quadrature CASA user-Jacobian",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(plt1, n_to_param.(csan), csadata_iq[2 + 3], lab = "Quadrature CASA AD-Jacobian",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(
    plt1, n_to_param.(csan), csadata_iq[3 + 3], lab = raw"Quadrature CASA AD-$v^{T}J$ seeding",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(plt1, n_to_param.(csan), csadata_g[1], lab = "Gauss CASA AD-Jacobian",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(
    plt1, n_to_param.(csan), csadata_g[2], lab = raw"Gauss CASA AD-$v^{T}J$ seeding",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(plt1, n_to_param.(csan), csadata_g[1 + 2], lab = "GaussKronrod CASA AD-Jacobian",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(
    plt1, n_to_param.(csan), csadata_g[2 + 2], lab = raw"GaussKronrod CASA AD-$v^{T}J$ seeding",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(plt1, n_to_param.(numdiffn), numdiff, lab = "Numerical Differentiation",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
xaxis!(plt1, "Number of Parameters", :log10);
yaxis!(plt1, "Runtime (s)", :log10);
plot!(plt1, legend = :outertopleft, size = (1200, 600))

VJP Choice Benchmarks

bt = 0:0.1:1
tspan = (0.0, 1.0)
csan = vcat(2:10, 12, 15, 17)
tols = (abstol = 1e-5, reltol = 1e-7)

_adjoint_methods = ntuple(4) do ii
    Alg = (InterpolatingAdjoint, QuadratureAdjoint, GaussAdjoint, GaussKronrodAdjoint)[ii]
    (
        advj1 = Alg(autodiff = true, autojacvec = EnzymeVJP()), # AD vJ (Enzyme)
        advj2 = Alg(autodiff = true, autojacvec = ReverseDiffVJP(false)), # AD vJ (ReverseDiff)
        advj3 = Alg(autodiff = true, autojacvec = ReverseDiffVJP(true)), # AD vJ (Compiled ReverseDiff)
        advj4 = Alg(autodiff = true, autojacvec = SciMLSensitivity.MooncakeVJP()) # AD vJ (Mooncake)
    )
end |> NamedTuple{(:interp, :quad, :gauss, :gausskronrod)}
adjoint_methods = mapreduce(collect, vcat, _adjoint_methods)
16-element Vector{SciMLBase.AbstractAdjointSensitivityAlgorithm{0, true, Va
l{:central}}}:
 SciMLSensitivity.InterpolatingAdjoint{0, true, Val{:central}, SciMLSensiti
vity.EnzymeVJP{EnzymeCore.ReverseMode{false, false, false, EnzymeCore.FFIAB
I, false, false}}}(SciMLSensitivity.EnzymeVJP{EnzymeCore.ReverseMode{false,
 false, false, EnzymeCore.FFIABI, false, false}}(0, EnzymeCore.ReverseMode{
false, false, false, EnzymeCore.FFIABI, false, false}()), false, false)
 SciMLSensitivity.InterpolatingAdjoint{0, true, Val{:central}, SciMLSensiti
vity.ReverseDiffVJP{false}}(SciMLSensitivity.ReverseDiffVJP{false}(), false
, false)
 SciMLSensitivity.InterpolatingAdjoint{0, true, Val{:central}, SciMLSensiti
vity.ReverseDiffVJP{true}}(SciMLSensitivity.ReverseDiffVJP{true}(), false, 
false)
 SciMLSensitivity.InterpolatingAdjoint{0, true, Val{:central}, SciMLSensiti
vity.MooncakeVJP}(SciMLSensitivity.MooncakeVJP(), false, false)
 SciMLSensitivity.QuadratureAdjoint{0, true, Val{:central}, SciMLSensitivit
y.EnzymeVJP{EnzymeCore.ReverseMode{false, false, false, EnzymeCore.FFIABI, 
false, false}}, Val{true}}(SciMLSensitivity.EnzymeVJP{EnzymeCore.ReverseMod
e{false, false, false, EnzymeCore.FFIABI, false, false}}(0, EnzymeCore.Reve
rseMode{false, false, false, EnzymeCore.FFIABI, false, false}()), 1.0e-6, 0
.001, Val{true}())
 SciMLSensitivity.QuadratureAdjoint{0, true, Val{:central}, SciMLSensitivit
y.ReverseDiffVJP{false}, Val{true}}(SciMLSensitivity.ReverseDiffVJP{false}(
), 1.0e-6, 0.001, Val{true}())
 SciMLSensitivity.QuadratureAdjoint{0, true, Val{:central}, SciMLSensitivit
y.ReverseDiffVJP{true}, Val{true}}(SciMLSensitivity.ReverseDiffVJP{true}(),
 1.0e-6, 0.001, Val{true}())
 SciMLSensitivity.QuadratureAdjoint{0, true, Val{:central}, SciMLSensitivit
y.MooncakeVJP, Val{true}}(SciMLSensitivity.MooncakeVJP(), 1.0e-6, 0.001, Va
l{true}())
 SciMLSensitivity.GaussAdjoint{0, true, Val{:central}, SciMLSensitivity.Enz
ymeVJP{EnzymeCore.ReverseMode{false, false, false, EnzymeCore.FFIABI, false
, false}}, Val{true}}(SciMLSensitivity.EnzymeVJP{EnzymeCore.ReverseMode{fal
se, false, false, EnzymeCore.FFIABI, false, false}}(0, EnzymeCore.ReverseMo
de{false, false, false, EnzymeCore.FFIABI, false, false}()), false, Val{tru
e}())
 SciMLSensitivity.GaussAdjoint{0, true, Val{:central}, SciMLSensitivity.Rev
erseDiffVJP{false}, Val{true}}(SciMLSensitivity.ReverseDiffVJP{false}(), fa
lse, Val{true}())
 SciMLSensitivity.GaussAdjoint{0, true, Val{:central}, SciMLSensitivity.Rev
erseDiffVJP{true}, Val{true}}(SciMLSensitivity.ReverseDiffVJP{true}(), fals
e, Val{true}())
 SciMLSensitivity.GaussAdjoint{0, true, Val{:central}, SciMLSensitivity.Moo
ncakeVJP, Val{true}}(SciMLSensitivity.MooncakeVJP(), false, Val{true}())
 SciMLSensitivity.GaussKronrodAdjoint{0, true, Val{:central}, SciMLSensitiv
ity.EnzymeVJP{EnzymeCore.ReverseMode{false, false, false, EnzymeCore.FFIABI
, false, false}}}(SciMLSensitivity.EnzymeVJP{EnzymeCore.ReverseMode{false, 
false, false, EnzymeCore.FFIABI, false, false}}(0, EnzymeCore.ReverseMode{f
alse, false, false, EnzymeCore.FFIABI, false, false}()), false)
 SciMLSensitivity.GaussKronrodAdjoint{0, true, Val{:central}, SciMLSensitiv
ity.ReverseDiffVJP{false}}(SciMLSensitivity.ReverseDiffVJP{false}(), false)
 SciMLSensitivity.GaussKronrodAdjoint{0, true, Val{:central}, SciMLSensitiv
ity.ReverseDiffVJP{true}}(SciMLSensitivity.ReverseDiffVJP{true}(), false)
 SciMLSensitivity.GaussKronrodAdjoint{0, true, Val{:central}, SciMLSensitiv
ity.MooncakeVJP}(SciMLSensitivity.MooncakeVJP(), false)

Warmup: compile all VJP backends before benchmarking.

let n = first(csan)
    bfun, b_u0, b_p, brusselator_jac, brusselator_comp = makebrusselator!(PROBS, n)
    solver = Rodas5(autodiff = AutoFiniteDiff())
    for alg in adjoint_methods
        f = SciMLSensitivity.alg_autodiff(alg) ? bfun :
            ODEFunction(bfun, jac = brusselator_jac)
        diffeq_sen_l2(f, b_u0, tspan, b_p, bt, solver; sensalg = alg, tols...)
    end
end
csavjp = map(csan) do n
    bfun, b_u0, b_p, brusselator_jac, brusselator_comp = makebrusselator!(PROBS, n)
    @time ts = map(adjoint_methods) do alg
        @info "Running $alg"
        f = SciMLSensitivity.alg_autodiff(alg) ? bfun :
            ODEFunction(bfun, jac = brusselator_jac)
        solver = Rodas5(autodiff = AutoFiniteDiff())
        @time diffeq_sen_l2(f, b_u0, tspan, b_p, bt, solver; sensalg = alg, tols...)
        t = @elapsed diffeq_sen_l2(f, b_u0, tspan, b_p, bt, solver; sensalg = alg, tols...)
        return t
    end
    @show n, ts
    ts
end
0.003498 seconds (7.34 k allocations: 569.445 KiB)
  0.102555 seconds (567.59 k allocations: 23.797 MiB)
  0.008709 seconds (3.63 k allocations: 295.617 KiB)
  1.244756 seconds (2.71 M allocations: 123.947 MiB, 4.23% gc time, 63.88% 
compilation time: <1% of which was recompilation)
  0.002558 seconds (7.90 k allocations: 499.922 KiB)
  0.103686 seconds (450.86 k allocations: 18.901 MiB)
  0.007000 seconds (4.40 k allocations: 307.938 KiB)
  0.002784 seconds (3.98 k allocations: 474.500 KiB)
  0.002768 seconds (6.17 k allocations: 488.875 KiB)
  0.071495 seconds (422.27 k allocations: 17.729 MiB)
  0.007003 seconds (4.74 k allocations: 343.312 KiB)
  0.002923 seconds (4.32 k allocations: 520.328 KiB)
  0.003209 seconds (8.92 k allocations: 575.453 KiB)
  0.189136 seconds (500.39 k allocations: 20.936 MiB, 52.22% gc time)
  0.008444 seconds (6.82 k allocations: 386.578 KiB)
  0.003366 seconds (6.67 k allocations: 604.844 KiB)
  3.767513 seconds (7.97 M allocations: 363.588 MiB, 4.02% gc time, 62.42% 
compilation time: <1% of which was recompilation)
(n, ts) = (2, [0.00357933, 0.094183289, 0.008210223, 0.003272981, 0.0022542
17, 0.09384601, 0.006543322, 0.002340117, 0.002430766, 0.066327519, 0.00644
4913, 0.002509646, 0.002863724, 0.103941123, 0.007919694, 0.002896453])
  0.010569 seconds (10.96 k allocations: 994.891 KiB)
  0.555127 seconds (2.11 M allocations: 92.709 MiB, 33.28% gc time)
  0.029560 seconds (5.10 k allocations: 556.945 KiB)
  0.009065 seconds (6.10 k allocations: 993.852 KiB)
  0.005754 seconds (9.10 k allocations: 695.078 KiB)
  0.299815 seconds (1.29 M allocations: 56.839 MiB, 25.79% gc time)
  0.018645 seconds (7.19 k allocations: 556.766 KiB)
  0.005163 seconds (4.45 k allocations: 671.109 KiB)
  0.005867 seconds (7.19 k allocations: 682.156 KiB)
  0.184316 seconds (1.13 M allocations: 49.919 MiB)
  0.018047 seconds (7.72 k allocations: 615.141 KiB)
  0.005472 seconds (4.83 k allocations: 716.219 KiB)
  0.007803 seconds (13.23 k allocations: 910.031 KiB)
  0.298279 seconds (1.62 M allocations: 71.162 MiB)
  0.024628 seconds (11.85 k allocations: 717.047 KiB)
  0.007332 seconds (9.72 k allocations: 938.125 KiB)
  2.856632 seconds (12.51 M allocations: 561.239 MiB, 11.61% gc time)
(n, ts) = (3, [0.01030296, 0.379802027, 0.029382481, 0.008891039, 0.0055839
58, 0.250016803, 0.018248305, 0.005019741, 0.005721407, 0.238887888, 0.0176
22169, 0.005224779, 0.007934165, 0.331477825, 0.02438136, 0.007273548])
  0.021788 seconds (16.77 k allocations: 1.758 MiB)
  1.192067 seconds (6.26 M allocations: 263.560 MiB, 5.35% gc time)
  0.091971 seconds (7.08 k allocations: 1.045 MiB)
  0.023078 seconds (8.47 k allocations: 1.764 MiB)
  0.009935 seconds (11.39 k allocations: 1005.797 KiB)
  0.627437 seconds (3.28 M allocations: 138.065 MiB, 5.31% gc time)
  0.044597 seconds (11.13 k allocations: 893.688 KiB)
  0.011113 seconds (5.44 k allocations: 1.003 MiB)
  0.010339 seconds (9.26 k allocations: 1.067 MiB)
  0.585433 seconds (2.90 M allocations: 122.115 MiB, 10.85% gc time)
  0.040941 seconds (11.68 k allocations: 1.042 MiB)
  0.011588 seconds (5.75 k allocations: 1.145 MiB)
  0.012650 seconds (18.10 k allocations: 1.426 MiB)
  0.809978 seconds (4.19 M allocations: 176.445 MiB, 4.12% gc time)
  0.057370 seconds (17.60 k allocations: 1.214 MiB)
  0.015760 seconds (12.84 k allocations: 1.495 MiB)
  7.227587 seconds (33.54 M allocations: 1.399 GiB, 5.90% gc time)
(n, ts) = (4, [0.021462316, 1.278853252, 0.089949413, 0.022928978, 0.009916
713, 0.614839037, 0.044023276, 0.011074977, 0.010255791, 0.554122346, 0.040
514827, 0.011186045, 0.012591048, 0.844452788, 0.057359251, 0.015506671])
  0.045835 seconds (23.80 k allocations: 3.015 MiB)
  2.780436 seconds (14.56 M allocations: 645.061 MiB, 6.61% gc time)
  0.197901 seconds (9.63 k allocations: 1.973 MiB)
  0.052093 seconds (11.29 k allocations: 3.000 MiB)
  0.018943 seconds (14.40 k allocations: 1.447 MiB)
  1.412798 seconds (7.16 M allocations: 317.337 MiB, 8.79% gc time)
  0.097020 seconds (16.21 k allocations: 1.387 MiB)
  0.022401 seconds (6.66 k allocations: 1.460 MiB)
  0.019573 seconds (11.74 k allocations: 1.742 MiB)
  1.226405 seconds (6.20 M allocations: 275.096 MiB, 8.09% gc time)
  0.085565 seconds (16.75 k allocations: 1.802 MiB)
  0.022608 seconds (6.74 k allocations: 1.815 MiB)
  0.022685 seconds (20.58 k allocations: 2.155 MiB)
  1.587518 seconds (8.21 M allocations: 364.006 MiB, 8.05% gc time)
  0.124561 seconds (22.67 k allocations: 2.027 MiB)
  0.028411 seconds (13.83 k allocations: 2.218 MiB)
 15.427244 seconds (72.62 M allocations: 3.177 GiB, 6.21% gc time)
(n, ts) = (5, [0.045540548, 2.782373351, 0.200580418, 0.051970631, 0.018867
821, 1.359150809, 0.096894564, 0.022643739, 0.019361058, 1.17461156, 0.0852
44371, 0.022409221, 0.02249108, 1.601062699, 0.127334138, 0.028251527])
  0.087212 seconds (32.42 k allocations: 5.000 MiB)
  5.558120 seconds (29.36 M allocations: 1.230 GiB, 5.91% gc time)
  0.445001 seconds (12.76 k allocations: 3.543 MiB)
  0.102894 seconds (14.84 k allocations: 5.060 MiB)
  0.033587 seconds (17.82 k allocations: 2.003 MiB)
  2.607524 seconds (13.65 M allocations: 585.184 MiB, 6.07% gc time)
  0.182571 seconds (22.39 k allocations: 1.996 MiB)
  0.041203 seconds (8.16 k allocations: 2.207 MiB)
  0.035265 seconds (14.75 k allocations: 2.842 MiB)
  2.240003 seconds (11.83 M allocations: 508.288 MiB, 5.77% gc time)
  0.199605 seconds (22.97 k allocations: 2.982 MiB)
  0.041603 seconds (8.11 k allocations: 3.107 MiB)
  0.040392 seconds (27.09 k allocations: 3.487 MiB)
  2.990970 seconds (15.96 M allocations: 685.151 MiB, 5.59% gc time)
  0.220237 seconds (31.13 k allocations: 3.359 MiB)
  0.053225 seconds (17.94 k allocations: 3.739 MiB)
 29.804666 seconds (142.04 M allocations: 6.013 GiB, 5.46% gc time)
(n, ts) = (6, [0.087207949, 5.607523723, 0.447168389, 0.10261769, 0.0331982
29, 2.551256747, 0.18263494, 0.041244943, 0.034492921, 2.230411441, 0.18640
4759, 0.042415006, 0.040155199, 3.037685831, 0.221543037, 0.052630757])
  0.170501 seconds (44.24 k allocations: 7.964 MiB)
 10.778985 seconds (55.63 M allocations: 2.474 GiB, 6.31% gc time)
  1.051191 seconds (16.48 k allocations: 5.938 MiB)
  0.214245 seconds (19.60 k allocations: 7.989 MiB)
  0.053761 seconds (21.85 k allocations: 2.693 MiB)
  4.624166 seconds (23.94 M allocations: 1.064 GiB, 6.20% gc time)
  0.321290 seconds (29.70 k allocations: 2.794 MiB)
  0.071469 seconds (9.77 k allocations: 2.887 MiB)
  0.056666 seconds (18.22 k allocations: 4.432 MiB)
  4.050053 seconds (20.72 M allocations: 945.120 MiB, 6.46% gc time)
  0.291727 seconds (30.28 k allocations: 4.716 MiB)
  0.072929 seconds (9.50 k allocations: 4.688 MiB)
  0.064064 seconds (32.50 k allocations: 5.265 MiB)
  5.405602 seconds (27.37 M allocations: 1.218 GiB, 6.01% gc time)
  0.488717 seconds (39.80 k allocations: 5.243 MiB)
  0.090315 seconds (20.70 k allocations: 5.495 MiB)
 55.441126 seconds (255.93 M allocations: 11.478 GiB, 5.65% gc time)
(n, ts) = (7, [0.173051145, 10.801715073, 1.055847101, 0.212338219, 0.05302
9055, 4.600069078, 0.323776779, 0.07120568, 0.056372875, 3.983449361, 0.293
602012, 0.072166926, 0.063923742, 5.246410839, 0.508434177, 0.089762314])
  0.284392 seconds (56.40 k allocations: 12.095 MiB)
 17.969693 seconds (93.83 M allocations: 4.048 GiB, 6.12% gc time)
  1.449346 seconds (20.71 k allocations: 9.473 MiB)
  0.366354 seconds (24.48 k allocations: 12.083 MiB)
  0.084212 seconds (27.71 k allocations: 3.636 MiB)
  7.998531 seconds (40.66 M allocations: 1.752 GiB, 10.71% gc time)
  0.716441 seconds (38.62 k allocations: 3.761 MiB)
  0.119827 seconds (12.41 k allocations: 3.817 MiB)
  0.087573 seconds (22.29 k allocations: 6.766 MiB)
  6.562415 seconds (34.12 M allocations: 1.474 GiB, 5.94% gc time)
  0.478873 seconds (38.73 k allocations: 7.160 MiB)
  0.118159 seconds (11.14 k allocations: 7.013 MiB)
  0.098195 seconds (40.35 k allocations: 7.927 MiB)
  8.714709 seconds (45.24 M allocations: 1.953 GiB, 6.06% gc time)
  0.639200 seconds (51.03 k allocations: 7.950 MiB)
  0.145010 seconds (24.90 k allocations: 8.106 MiB)
 91.318967 seconds (428.43 M allocations: 18.630 GiB, 5.91% gc time)
(n, ts) = (8, [0.283189062, 18.005933975, 1.358903657, 0.362674095, 0.08355
4489, 7.692824616, 0.61045138, 0.119612362, 0.087259919, 6.595265525, 0.476
320281, 0.117615683, 0.097640789, 8.786335571, 0.630147916, 0.145100126])
  0.683006 seconds (107.37 k allocations: 20.228 MiB)
 30.467188 seconds (160.08 M allocations: 6.759 GiB, 5.95% gc time)
  2.321723 seconds (25.49 k allocations: 14.548 MiB)
  0.630827 seconds (31.97 k allocations: 18.106 MiB)
  0.137693 seconds (33.16 k allocations: 4.793 MiB)
 12.208722 seconds (63.44 M allocations: 2.676 GiB, 8.78% gc time)
  0.861986 seconds (48.17 k allocations: 4.988 MiB)
  0.189499 seconds (14.59 k allocations: 4.984 MiB)
  0.144530 seconds (26.88 k allocations: 10.085 MiB)
 10.231270 seconds (53.29 M allocations: 2.254 GiB, 5.92% gc time)
  0.753061 seconds (48.28 k allocations: 10.603 MiB)
  0.189602 seconds (12.98 k allocations: 10.345 MiB)
  0.160288 seconds (47.70 k allocations: 11.604 MiB)
 13.058479 seconds (68.74 M allocations: 2.905 GiB, 5.75% gc time)
  1.232862 seconds (61.73 k allocations: 11.646 MiB)
  0.230444 seconds (29.22 k allocations: 11.814 MiB)
146.952136 seconds (692.06 M allocations: 29.450 GiB, 5.63% gc time)
(n, ts) = (9, [0.717109066, 30.624253989, 2.359922965, 0.633421316, 0.13889
2631, 11.906786011, 0.864173848, 0.189606119, 0.144029821, 10.228031614, 0.
756191979, 0.189723228, 0.193370688, 13.122467513, 1.118886682, 0.230099266
])
  0.822824 seconds (88.59 k allocations: 25.620 MiB, 4.68% gc time)
 46.350225 seconds (234.14 M allocations: 9.730 GiB, 6.72% gc time)
  3.584621 seconds (30.84 k allocations: 21.377 MiB, 0.95% gc time)
  1.028322 seconds (37.36 k allocations: 25.523 MiB)
  0.193850 seconds (39.23 k allocations: 6.158 MiB)
 18.134978 seconds (94.85 M allocations: 3.937 GiB, 5.96% gc time)
  1.275175 seconds (58.85 k allocations: 6.430 MiB)
  0.284454 seconds (17.02 k allocations: 6.333 MiB)
  0.216525 seconds (34.16 k allocations: 14.581 MiB)
 16.495200 seconds (85.57 M allocations: 3.561 GiB, 5.62% gc time)
  1.438721 seconds (58.97 k allocations: 15.100 MiB)
  0.296940 seconds (15.90 k allocations: 14.821 MiB)
  0.234178 seconds (55.54 k allocations: 16.333 MiB)
 20.399234 seconds (104.93 M allocations: 4.364 GiB, 5.81% gc time)
  1.462550 seconds (72.68 k allocations: 16.351 MiB)
  0.347343 seconds (32.59 k allocations: 16.521 MiB)
224.981678 seconds (1.04 G allocations: 43.549 GiB, 5.58% gc time)
(n, ts) = (10, [0.781561902, 45.772287025, 3.49693886, 1.074921273, 0.19360
6375, 18.216835184, 1.282054869, 0.407848052, 0.215832087, 16.642934089, 1.
648352211, 0.32761436, 0.233865626, 20.234351809, 1.464208373, 0.38271255])
  2.064751 seconds (130.62 k allocations: 48.696 MiB)
 95.946192 seconds (499.22 M allocations: 21.553 GiB, 6.41% gc time)
  8.372802 seconds (43.90 k allocations: 42.393 MiB, 0.37% gc time)
  3.185799 seconds (54.66 k allocations: 48.844 MiB, 5.91% gc time)
  0.485238 seconds (53.61 k allocations: 9.886 MiB)
 36.251090 seconds (191.77 M allocations: 8.270 GiB, 6.48% gc time)
  2.637927 seconds (83.88 k allocations: 10.438 MiB)
  0.636772 seconds (23.09 k allocations: 10.766 MiB, 5.56% gc time)
  0.473279 seconds (46.18 k allocations: 27.760 MiB)
 32.531326 seconds (167.31 M allocations: 7.234 GiB, 7.07% gc time)
  2.515598 seconds (84.47 k allocations: 28.740 MiB, 1.42% gc time)
  0.749088 seconds (21.31 k allocations: 28.713 MiB, 10.58% gc time)
  0.527538 seconds (72.60 k allocations: 30.481 MiB)
 38.970308 seconds (204.15 M allocations: 8.823 GiB, 6.23% gc time)
  3.282901 seconds (102.30 k allocations: 30.938 MiB)
  0.754170 seconds (43.04 k allocations: 31.504 MiB)
458.154045 seconds (2.13 G allocations: 92.446 GiB, 6.02% gc time)
(n, ts) = (12, [2.126064473, 95.772762116, 8.394698919, 2.576371745, 0.4618
78529, 35.99736931, 2.654072956, 0.593293707, 0.55562367, 32.343255805, 2.6
51178352, 0.8568525, 0.532520398, 39.233807236, 3.173992531, 0.80275542])
  5.083548 seconds (207.92 k allocations: 109.801 MiB, 1.54% gc time)
243.174347 seconds (1.25 G allocations: 52.110 GiB, 6.54% gc time)
 23.247176 seconds (66.69 k allocations: 99.539 MiB, 0.15% gc time)
  7.965201 seconds (104.65 k allocations: 112.650 MiB, 3.26% gc time)
  1.066517 seconds (79.53 k allocations: 18.622 MiB)
 86.794626 seconds (458.25 M allocations: 19.037 GiB, 6.43% gc time)
  6.594752 seconds (129.40 k allocations: 19.505 MiB)
  1.497955 seconds (33.47 k allocations: 19.487 MiB)
  1.318228 seconds (73.53 k allocations: 63.577 MiB, 2.40% gc time)
 80.848889 seconds (414.43 M allocations: 17.262 GiB, 6.69% gc time)
  7.599819 seconds (130.02 k allocations: 64.799 MiB, 0.47% gc time)
  1.830667 seconds (31.34 k allocations: 64.371 MiB, 4.05% gc time)
  1.408392 seconds (104.84 k allocations: 68.044 MiB, 2.67% gc time)
 93.825254 seconds (484.08 M allocations: 20.156 GiB, 6.32% gc time)
  9.372143 seconds (151.47 k allocations: 68.739 MiB, 0.36% gc time)
  1.988884 seconds (58.26 k allocations: 69.166 MiB, 2.31% gc time)
1137.794919 seconds (5.22 G allocations: 218.650 GiB, 5.83% gc time)
(n, ts) = (15, [4.405508722, 240.768219103, 21.092031701, 7.979428498, 1.09
4051443, 86.384456247, 6.628444017, 1.534533918, 1.435663909, 80.244023123,
 7.062599139, 1.932887896, 1.386303152, 92.331647464, 7.764801673, 2.087719
974])
  8.602305 seconds (248.09 k allocations: 173.692 MiB, 1.99% gc time)
376.438230 seconds (1.92 G allocations: 83.774 GiB, 6.53% gc time)
 34.622631 seconds (84.69 k allocations: 161.743 MiB, 0.18% gc time)
 10.703464 seconds (101.67 k allocations: 173.513 MiB, 1.32% gc time)
  1.871627 seconds (100.02 k allocations: 27.201 MiB, 1.80% gc time)
142.471792 seconds (749.62 M allocations: 32.584 GiB, 7.05% gc time)
 11.974254 seconds (165.37 k allocations: 28.551 MiB, 0.29% gc time)
  2.809907 seconds (41.66 k allocations: 28.016 MiB, 2.32% gc time)
  2.413948 seconds (89.78 k allocations: 101.685 MiB, 3.29% gc time)
130.423152 seconds (678.24 M allocations: 29.556 GiB, 6.78% gc time)
 10.619360 seconds (166.01 k allocations: 103.646 MiB, 0.36% gc time)
  3.326598 seconds (38.78 k allocations: 102.582 MiB, 6.32% gc time)
  2.428297 seconds (126.33 k allocations: 108.000 MiB, 2.83% gc time)
149.670914 seconds (772.63 M allocations: 33.660 GiB, 6.92% gc time)
 15.474737 seconds (188.60 k allocations: 108.863 MiB, 0.26% gc time)
  3.361333 seconds (66.17 k allocations: 108.531 MiB, 1.31% gc time)
1806.170842 seconds (8.25 G allocations: 361.545 GiB, 6.05% gc time)
(n, ts) = (17, [8.684872401, 370.445731333, 35.454873012, 10.900776718, 1.8
26238911, 140.608069171, 13.46019005, 2.655477592, 2.796960139, 130.2751240
64, 10.732904286, 3.118774368, 2.470396791, 149.926873048, 12.213299233, 3.
343376068])
12-element Vector{Vector{Float64}}:
 [0.00357933, 0.094183289, 0.008210223, 0.003272981, 0.002254217, 0.0938460
1, 0.006543322, 0.002340117, 0.002430766, 0.066327519, 0.006444913, 0.00250
9646, 0.002863724, 0.103941123, 0.007919694, 0.002896453]
 [0.01030296, 0.379802027, 0.029382481, 0.008891039, 0.005583958, 0.2500168
03, 0.018248305, 0.005019741, 0.005721407, 0.238887888, 0.017622169, 0.0052
24779, 0.007934165, 0.331477825, 0.02438136, 0.007273548]
 [0.021462316, 1.278853252, 0.089949413, 0.022928978, 0.009916713, 0.614839
037, 0.044023276, 0.011074977, 0.010255791, 0.554122346, 0.040514827, 0.011
186045, 0.012591048, 0.844452788, 0.057359251, 0.015506671]
 [0.045540548, 2.782373351, 0.200580418, 0.051970631, 0.018867821, 1.359150
809, 0.096894564, 0.022643739, 0.019361058, 1.17461156, 0.085244371, 0.0224
09221, 0.02249108, 1.601062699, 0.127334138, 0.028251527]
 [0.087207949, 5.607523723, 0.447168389, 0.10261769, 0.033198229, 2.5512567
47, 0.18263494, 0.041244943, 0.034492921, 2.230411441, 0.186404759, 0.04241
5006, 0.040155199, 3.037685831, 0.221543037, 0.052630757]
 [0.173051145, 10.801715073, 1.055847101, 0.212338219, 0.053029055, 4.60006
9078, 0.323776779, 0.07120568, 0.056372875, 3.983449361, 0.293602012, 0.072
166926, 0.063923742, 5.246410839, 0.508434177, 0.089762314]
 [0.283189062, 18.005933975, 1.358903657, 0.362674095, 0.083554489, 7.69282
4616, 0.61045138, 0.119612362, 0.087259919, 6.595265525, 0.476320281, 0.117
615683, 0.097640789, 8.786335571, 0.630147916, 0.145100126]
 [0.717109066, 30.624253989, 2.359922965, 0.633421316, 0.138892631, 11.9067
86011, 0.864173848, 0.189606119, 0.144029821, 10.228031614, 0.756191979, 0.
189723228, 0.193370688, 13.122467513, 1.118886682, 0.230099266]
 [0.781561902, 45.772287025, 3.49693886, 1.074921273, 0.193606375, 18.21683
5184, 1.282054869, 0.407848052, 0.215832087, 16.642934089, 1.648352211, 0.3
2761436, 0.233865626, 20.234351809, 1.464208373, 0.38271255]
 [2.126064473, 95.772762116, 8.394698919, 2.576371745, 0.461878529, 35.9973
6931, 2.654072956, 0.593293707, 0.55562367, 32.343255805, 2.651178352, 0.85
68525, 0.532520398, 39.233807236, 3.173992531, 0.80275542]
 [4.405508722, 240.768219103, 21.092031701, 7.979428498, 1.094051443, 86.38
4456247, 6.628444017, 1.534533918, 1.435663909, 80.244023123, 7.062599139, 
1.932887896, 1.386303152, 92.331647464, 7.764801673, 2.087719974]
 [8.684872401, 370.445731333, 35.454873012, 10.900776718, 1.826238911, 140.
608069171, 13.46019005, 2.655477592, 2.796960139, 130.275124064, 10.7329042
86, 3.118774368, 2.470396791, 149.926873048, 12.213299233, 3.343376068]
csacompare = [[csavjp[j][i] for j in eachindex(csavjp)] for i in eachindex(csavjp[1])]

plt_interp = plot(title = "Brusselator interpolating adjoint VJP scaling");
plot!(plt_interp, n_to_param.(csan), csadata_iq[2], lab = "AD-Jacobian",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(plt_interp, n_to_param.(csan), csacompare[1], lab = raw"EnzymeVJP",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(plt_interp, n_to_param.(csan), csacompare[2], lab = raw"ReverseDiffVJP",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(plt_interp, n_to_param.(csan), csacompare[3], lab = raw"Compiled ReverseDiffVJP",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(plt_interp, n_to_param.(csan), csacompare[4], lab = raw"MooncakeVJP",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
xaxis!(plt_interp, "Number of Parameters", :log10);
yaxis!(plt_interp, "Runtime (s)", :log10);
plot!(plt_interp, legend = :outertopleft, size = (1200, 600))

plt2 = plot(title = "Brusselator quadrature adjoint VJP scaling");
plot!(plt2, n_to_param.(csan), csadata_iq[2 + 3], lab = "AD-Jacobian",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(plt2, n_to_param.(csan), csacompare[1 + 4], lab = raw"EnzymeVJP",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(plt2, n_to_param.(csan), csacompare[2 + 4], lab = raw"ReverseDiffVJP",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(plt2, n_to_param.(csan), csacompare[3 + 4], lab = raw"Compiled ReverseDiffVJP",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(plt2, n_to_param.(csan), csacompare[4 + 4], lab = raw"MooncakeVJP",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
xaxis!(plt2, "Number of Parameters", :log10);
yaxis!(plt2, "Runtime (s)", :log10);
plot!(plt2, legend = :outertopleft, size = (1200, 600))

plt_gauss = plot(title = "Brusselator Gauss adjoint VJP scaling");
plot!(plt_gauss, n_to_param.(csan), csadata_g[1], lab = "AD-Jacobian",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(plt_gauss, n_to_param.(csan), csacompare[1 + 8], lab = raw"EnzymeVJP",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(plt_gauss, n_to_param.(csan), csacompare[2 + 8], lab = raw"ReverseDiffVJP",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(plt_gauss, n_to_param.(csan), csacompare[3 + 8], lab = raw"Compiled ReverseDiffVJP",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(plt_gauss, n_to_param.(csan), csacompare[4 + 8], lab = raw"MooncakeVJP",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
xaxis!(plt_gauss, "Number of Parameters", :log10);
yaxis!(plt_gauss, "Runtime (s)", :log10);
plot!(plt_gauss, legend = :outertopleft, size = (1200, 600))

plt_gk = plot(title = "Brusselator GaussKronrod adjoint VJP scaling");
plot!(plt_gk, n_to_param.(csan), csadata_g[1 + 2], lab = "AD-Jacobian",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(plt_gk, n_to_param.(csan), csacompare[1 + 12], lab = raw"EnzymeVJP",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(plt_gk, n_to_param.(csan), csacompare[2 + 12], lab = raw"ReverseDiffVJP",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(plt_gk, n_to_param.(csan), csacompare[3 + 12], lab = raw"Compiled ReverseDiffVJP",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(plt_gk, n_to_param.(csan), csacompare[4 + 12], lab = raw"MooncakeVJP",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
xaxis!(plt_gk, "Number of Parameters", :log10);
yaxis!(plt_gk, "Runtime (s)", :log10);
plot!(plt_gk, legend = :outertopleft, size = (1200, 600))

SUNDIALS CVODES C Adjoint Benchmarks

SundialsAdjoint uses the SUNDIALS CVODES C adjoint interface (CVodeAdjInit/CVodeF/CVodeB): the forward pass is re-integrated by CVODES with checkpointing, and both the backward (adjoint) pass and the parameter-gradient quadrature are integrated by CVODES itself. The vector-Jacobian products inside the backward pass use the same autojacvec machinery as the native Julia adjoints, so comparing against GaussAdjoint with the same vjp choice measures the difference of the surrounding ODE-integration and checkpointing machinery (CVODES in C vs the native Julia adjoints). Note that SundialsAdjoint requires CVODE_BDF/CVODE_Adams as the solver, while the native adjoint rows above use Rodas5, so the forward solver differs as well.

using Sundials

sundials_configs = [
    ("CVODES Dense EnzymeVJP", CVODE_BDF(),
        SundialsAdjoint(autojacvec = EnzymeVJP())),
    ("CVODES Dense Compiled ReverseDiffVJP", CVODE_BDF(),
        SundialsAdjoint(autojacvec = ReverseDiffVJP(true))),
    ("CVODES GMRES EnzymeVJP", CVODE_BDF(linear_solver = :GMRES),
        SundialsAdjoint(autojacvec = EnzymeVJP())),
    ("CVODES GMRES Compiled ReverseDiffVJP", CVODE_BDF(linear_solver = :GMRES),
        SundialsAdjoint(autojacvec = ReverseDiffVJP(true))),
]
4-element Vector{Tuple{String, Sundials.CVODE_BDF{:Newton, LinearSolver, No
thing, Nothing} where LinearSolver, SciMLSensitivity.SundialsAdjoint{0, tru
e, Val{:central}}}}:
 ("CVODES Dense EnzymeVJP", Sundials.CVODE_BDF{:Newton, :Dense, Nothing, No
thing}(0, 0, 0, false, 10, 5, 7, 3, 10, nothing, nothing, 0), SciMLSensitiv
ity.SundialsAdjoint{0, true, Val{:central}, SciMLSensitivity.EnzymeVJP{Enzy
meCore.ReverseMode{false, false, false, EnzymeCore.FFIABI, false, false}}}(
SciMLSensitivity.EnzymeVJP{EnzymeCore.ReverseMode{false, false, false, Enzy
meCore.FFIABI, false, false}}(0, EnzymeCore.ReverseMode{false, false, false
, EnzymeCore.FFIABI, false, false}()), 150, :hermite, true))
 ("CVODES Dense Compiled ReverseDiffVJP", Sundials.CVODE_BDF{:Newton, :Dens
e, Nothing, Nothing}(0, 0, 0, false, 10, 5, 7, 3, 10, nothing, nothing, 0),
 SciMLSensitivity.SundialsAdjoint{0, true, Val{:central}, SciMLSensitivity.
ReverseDiffVJP{true}}(SciMLSensitivity.ReverseDiffVJP{true}(), 150, :hermit
e, true))
 ("CVODES GMRES EnzymeVJP", Sundials.CVODE_BDF{:Newton, :GMRES, Nothing, No
thing}(0, 0, 0, false, 10, 5, 7, 3, 10, nothing, nothing, 0), SciMLSensitiv
ity.SundialsAdjoint{0, true, Val{:central}, SciMLSensitivity.EnzymeVJP{Enzy
meCore.ReverseMode{false, false, false, EnzymeCore.FFIABI, false, false}}}(
SciMLSensitivity.EnzymeVJP{EnzymeCore.ReverseMode{false, false, false, Enzy
meCore.FFIABI, false, false}}(0, EnzymeCore.ReverseMode{false, false, false
, EnzymeCore.FFIABI, false, false}()), 150, :hermite, true))
 ("CVODES GMRES Compiled ReverseDiffVJP", Sundials.CVODE_BDF{:Newton, :GMRE
S, Nothing, Nothing}(0, 0, 0, false, 10, 5, 7, 3, 10, nothing, nothing, 0),
 SciMLSensitivity.SundialsAdjoint{0, true, Val{:central}, SciMLSensitivity.
ReverseDiffVJP{true}}(SciMLSensitivity.ReverseDiffVJP{true}(), 150, :hermit
e, true))

Check that the CVODES C adjoint returns the same gradient as the native GaussAdjoint before timing it:

let n = first(csan)
    bfun, b_u0, b_p, brusselator_jac, brusselator_comp = makebrusselator!(PROBS, n)
    solver = Rodas5(autodiff = AutoFiniteDiff())
    gauss = GaussAdjoint(autodiff = true, autojacvec = EnzymeVJP())
    du0_ref, dp_ref = diffeq_sen_l2(
        bfun, b_u0, tspan, b_p, bt, solver; sensalg = gauss, tols...)
    for (name, alg, sensealg) in sundials_configs
        du0, dp = diffeq_sen_l2(
            bfun, b_u0, tspan, b_p, bt, alg; sensalg = sensealg, tols...)
        err_du0 = norm(du0 - du0_ref) / norm(du0_ref)
        err_dp = norm(vec(dp) - vec(dp_ref)) / norm(vec(dp_ref))
        @show name, err_du0, err_dp
    end
end
(name, err_du0, err_dp) = ("CVODES Dense EnzymeVJP", 1.576392245414933e-6, 
1.1362910965483166e-5)
(name, err_du0, err_dp) = ("CVODES Dense Compiled ReverseDiffVJP", 1.576392
118216258e-6, 1.1362910575670025e-5)
(name, err_du0, err_dp) = ("CVODES GMRES EnzymeVJP", 1.4691746698600681e-6,
 1.08063323821113e-5)
(name, err_du0, err_dp) = ("CVODES GMRES Compiled ReverseDiffVJP", 1.469174
542203823e-6, 1.0806332171648812e-5)
csa_sundials = map(csan) do n
    bfun, b_u0, b_p, brusselator_jac, brusselator_comp = makebrusselator!(PROBS, n)
    @time ts = map(sundials_configs) do (name, alg, sensealg)
        @info "Running $name"
        @time diffeq_sen_l2(bfun, b_u0, tspan, b_p, bt, alg; sensalg = sensealg, tols...)
        t = @elapsed diffeq_sen_l2(
            bfun, b_u0, tspan, b_p, bt, alg; sensalg = sensealg, tols...)
        return t
    end
    @show n, ts
    ts
end
0.002473 seconds (7.28 k allocations: 319.281 KiB)
  0.005579 seconds (6.96 k allocations: 249.562 KiB)
  0.002838 seconds (9.43 k allocations: 405.531 KiB)
  0.006594 seconds (8.32 k allocations: 286.875 KiB)
  0.363509 seconds (511.97 k allocations: 25.129 MiB, 90.03% compilation ti
me)
(n, ts) = (2, [0.002086178, 0.005062802, 0.002773134, 0.006453053])
  0.005205 seconds (11.21 k allocations: 475.281 KiB)
  0.015527 seconds (12.46 k allocations: 461.000 KiB)
  0.006457 seconds (15.21 k allocations: 618.844 KiB)
  0.018666 seconds (15.56 k allocations: 545.438 KiB)
  0.091794 seconds (109.43 k allocations: 4.145 MiB)
(n, ts) = (3, [0.005060671, 0.015249293, 0.006486894, 0.018511246])
  0.008035 seconds (13.55 k allocations: 577.594 KiB)
  0.032119 seconds (17.89 k allocations: 684.312 KiB)
  0.008424 seconds (16.26 k allocations: 667.219 KiB)
  0.034156 seconds (20.04 k allocations: 741.656 KiB)
  0.166780 seconds (136.03 k allocations: 5.259 MiB)
(n, ts) = (4, [0.008281873, 0.032160047, 0.008517012, 0.034684933])
  0.013110 seconds (15.82 k allocations: 673.938 KiB)
  0.057480 seconds (24.53 k allocations: 983.406 KiB)
  0.011665 seconds (18.01 k allocations: 728.906 KiB)
  0.054962 seconds (26.79 k allocations: 1.018 MiB)
  0.269789 seconds (170.84 k allocations: 6.739 MiB)
(n, ts) = (5, [0.012997766, 0.054286512, 0.011679164, 0.053198159])
  0.022283 seconds (19.61 k allocations: 842.422 KiB)
  0.102413 seconds (33.08 k allocations: 1.304 MiB)
  0.017324 seconds (21.76 k allocations: 869.609 KiB)
  0.095864 seconds (35.90 k allocations: 1.373 MiB)
  0.483422 seconds (221.26 k allocations: 8.742 MiB)
(n, ts) = (6, [0.022679452, 0.107991858, 0.01750931, 0.096944701])
  0.033726 seconds (22.45 k allocations: 964.359 KiB)
  0.171887 seconds (42.32 k allocations: 1.732 MiB)
  0.021004 seconds (22.60 k allocations: 906.984 KiB)
  0.126830 seconds (44.07 k allocations: 1.771 MiB)
  0.695932 seconds (263.45 k allocations: 10.703 MiB)
(n, ts) = (7, [0.033919088, 0.153759, 0.02117437, 0.133210426])
  0.056150 seconds (28.34 k allocations: 1.193 MiB)
  0.293944 seconds (54.50 k allocations: 2.215 MiB)
  0.033031 seconds (35.38 k allocations: 1.257 MiB)
  0.184771 seconds (64.20 k allocations: 2.451 MiB)
  1.124065 seconds (365.40 k allocations: 14.276 MiB)
(n, ts) = (8, [0.05648885, 0.277458679, 0.032985403, 0.188745231])
  0.092530 seconds (34.60 k allocations: 1.442 MiB)
  0.394098 seconds (68.48 k allocations: 2.765 MiB)
  0.095622 seconds (36.72 k allocations: 1.328 MiB, 43.24% gc time, 9.77% c
ompilation time)
  0.245516 seconds (74.45 k allocations: 2.900 MiB)
  1.600173 seconds (428.99 k allocations: 16.912 MiB, 2.58% gc time, 0.58% 
compilation time)
(n, ts) = (9, [0.092093779, 0.388363181, 0.044032041, 0.247179641])
  0.141899 seconds (39.46 k allocations: 1.668 MiB)
  0.627596 seconds (81.90 k allocations: 3.322 MiB)
  0.053260 seconds (32.61 k allocations: 1.296 MiB)
  0.354271 seconds (79.64 k allocations: 3.245 MiB)
  2.370902 seconds (467.79 k allocations: 19.104 MiB)
(n, ts) = (10, [0.141058731, 0.60759204, 0.053315778, 0.391273765])
  0.330021 seconds (53.48 k allocations: 2.260 MiB)
  1.241289 seconds (115.78 k allocations: 4.784 MiB)
  0.107943 seconds (50.52 k allocations: 1.957 MiB)
  0.819694 seconds (119.00 k allocations: 4.847 MiB)
  5.061984 seconds (678.09 k allocations: 27.740 MiB)
(n, ts) = (12, [0.328350491, 1.26811746, 0.107561791, 0.858010232])
  1.040991 seconds (78.77 k allocations: 3.339 MiB)
  3.822698 seconds (177.69 k allocations: 7.285 MiB)
  0.190030 seconds (61.42 k allocations: 2.427 MiB)
  1.610268 seconds (172.19 k allocations: 7.103 MiB)
 13.307436 seconds (980.72 k allocations: 40.351 MiB)
(n, ts) = (15, [1.03666319, 3.914331815, 0.189125469, 1.501709407])
  1.882694 seconds (98.57 k allocations: 4.185 MiB)
  5.125975 seconds (226.44 k allocations: 9.461 MiB)
  0.249960 seconds (61.78 k allocations: 2.495 MiB)
  1.829694 seconds (207.07 k allocations: 8.895 MiB)
 18.292115 seconds (1.19 M allocations: 50.113 MiB, 0.49% gc time)
(n, ts) = (17, [1.88243637, 5.125844432, 0.248855081, 1.945141334])
12-element Vector{Vector{Float64}}:
 [0.002086178, 0.005062802, 0.002773134, 0.006453053]
 [0.005060671, 0.015249293, 0.006486894, 0.018511246]
 [0.008281873, 0.032160047, 0.008517012, 0.034684933]
 [0.012997766, 0.054286512, 0.011679164, 0.053198159]
 [0.022679452, 0.107991858, 0.01750931, 0.096944701]
 [0.033919088, 0.153759, 0.02117437, 0.133210426]
 [0.05648885, 0.277458679, 0.032985403, 0.188745231]
 [0.092093779, 0.388363181, 0.044032041, 0.247179641]
 [0.141058731, 0.60759204, 0.053315778, 0.391273765]
 [0.328350491, 1.26811746, 0.107561791, 0.858010232]
 [1.03666319, 3.914331815, 0.189125469, 1.501709407]
 [1.88243637, 5.125844432, 0.248855081, 1.945141334]
csadata_sundials = [[csa_sundials[j][i] for j in eachindex(csa_sundials)]
                    for i in eachindex(csa_sundials[1])]

plt_sundials = plot(title = "Brusselator CVODES C adjoint vs native adjoints");
plot!(plt_sundials, n_to_param.(csan), csadata_g[2],
    lab = raw"GaussAdjoint EnzymeVJP (Rodas5)",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(plt_sundials, n_to_param.(csan), csacompare[3 + 8],
    lab = raw"GaussAdjoint Compiled ReverseDiffVJP (Rodas5)",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
plot!(plt_sundials, n_to_param.(csan), csacompare[1],
    lab = raw"InterpolatingAdjoint EnzymeVJP (Rodas5)",
    lw = lw, marksize = ms, linestyle = :auto, marker = :auto);
for (i, (name, _, _)) in enumerate(sundials_configs)
    plot!(plt_sundials, n_to_param.(csan), csadata_sundials[i], lab = name,
        lw = lw, marksize = ms, linestyle = :auto, marker = :auto)
end
xaxis!(plt_sundials, "Number of Parameters", :log10);
yaxis!(plt_sundials, "Runtime (s)", :log10);
plot!(plt_sundials, legend = :outertopleft, size = (1200, 600))

Same-solver adjoint comparison

The plot above conflates the forward solver choice with the adjoint method. To isolate the adjoint machinery itself, compare the CVODES C adjoint against the native GaussAdjoint with the same ODE solver: CVODE_BDF with dense and GMRES linear solvers, using EnzymeVJP for both adjoints. FBDF (the native Julia fixed-leading-coefficient BDF analogue of CVODE) with GaussAdjoint is included as the all-Julia BDF counterpart, with dense LU and matrix-free KrylovJL_GMRES linear solvers.

GaussAdjoint with a Sundials solver currently requires a dense forward solution (with saveat the checkpointing path errors on the interpolation type), so the GaussAdjoint rows use dense = true forward solves while the SundialsAdjoint rows use the saveat forward solve plus the internal checkpointed CVodeF re-integration; each method is measured in its natural mode of use.

using LinearSolve

function diffeq_sen_l2_dense(df, u0, tspan, p, t, alg;
        abstol = 1e-5, reltol = 1e-7, sensalg, kwargs...)
    prob = ODEProblem{true, SciMLBase.FullSpecialize}(df, u0, tspan, p)
    sol = solve(prob, alg, abstol = abstol, reltol = reltol, dense = true; kwargs...)
    dg(out, u, p, t, i) = (out.=u .- 1.0)
    adjoint_sensitivities(sol, alg; t, abstol = abstol, dgdu_discrete = dg,
        reltol = reltol, sensealg = sensalg)
end

gauss_enz = GaussAdjoint(autodiff = true, autojacvec = EnzymeVJP())
sundials_enz = SundialsAdjoint(autojacvec = EnzymeVJP())

same_solver_configs = [
    ("CVODE_BDF Dense + SundialsAdjoint", CVODE_BDF(), sundials_enz, false),
    ("CVODE_BDF Dense + GaussAdjoint", CVODE_BDF(), gauss_enz, true),
    ("CVODE_BDF GMRES + SundialsAdjoint",
        CVODE_BDF(linear_solver = :GMRES), sundials_enz, false),
    ("CVODE_BDF GMRES + GaussAdjoint",
        CVODE_BDF(linear_solver = :GMRES), gauss_enz, true),
    ("FBDF + GaussAdjoint", FBDF(autodiff = AutoFiniteDiff()), gauss_enz, true),
    ("FBDF GMRES + GaussAdjoint",
        FBDF(linsolve = KrylovJL_GMRES(), autodiff = AutoFiniteDiff()), gauss_enz, true),
]

function run_same_solver(bfun, b_u0, b_p, alg, sensealg, dense_forward)
    return dense_forward ?
        diffeq_sen_l2_dense(bfun, b_u0, tspan, b_p, bt, alg; sensalg = sensealg, tols...) :
        diffeq_sen_l2(bfun, b_u0, tspan, b_p, bt, alg; sensalg = sensealg, tols...)
end
run_same_solver (generic function with 1 method)
let n = first(csan)
    bfun, b_u0, b_p, brusselator_jac, brusselator_comp = makebrusselator!(PROBS, n)
    solver = Rodas5(autodiff = AutoFiniteDiff())
    du0_ref, dp_ref = diffeq_sen_l2(
        bfun, b_u0, tspan, b_p, bt, solver; sensalg = gauss_enz, tols...)
    for (name, alg, sensealg, dense_forward) in same_solver_configs
        du0, dp = run_same_solver(bfun, b_u0, b_p, alg, sensealg, dense_forward)
        err_du0 = norm(du0 - du0_ref) / norm(du0_ref)
        err_dp = norm(vec(dp) - vec(dp_ref)) / norm(vec(dp_ref))
        @show name, err_du0, err_dp
    end
end
(name, err_du0, err_dp) = ("CVODE_BDF Dense + SundialsAdjoint", 1.576392245
414933e-6, 1.1362910965483166e-5)
(name, err_du0, err_dp) = ("CVODE_BDF Dense + GaussAdjoint", 1.572195778601
5652e-6, 1.8830159827863253e-5)
(name, err_du0, err_dp) = ("CVODE_BDF GMRES + SundialsAdjoint", 1.469174669
8600681e-6, 1.08063323821113e-5)
(name, err_du0, err_dp) = ("CVODE_BDF GMRES + GaussAdjoint", 1.578966610320
9478e-6, 1.8891961107264642e-5)
(name, err_du0, err_dp) = ("FBDF + GaussAdjoint", 4.1711912972772977e-7, 2.
3283254883689787e-5)
(name, err_du0, err_dp) = ("FBDF GMRES + GaussAdjoint", 9.621092548258918e-
7, 2.387569401387661e-5)
csa_same_solver = map(csan) do n
    bfun, b_u0, b_p, brusselator_jac, brusselator_comp = makebrusselator!(PROBS, n)
    @time ts = map(same_solver_configs) do (name, alg, sensealg, dense_forward)
        @info "Running $name"
        @time run_same_solver(bfun, b_u0, b_p, alg, sensealg, dense_forward)
        t = @elapsed run_same_solver(bfun, b_u0, b_p, alg, sensealg, dense_forward)
        return t
    end
    @show n, ts
    ts
end
0.002188 seconds (7.28 k allocations: 319.281 KiB)
  0.002764 seconds (7.56 k allocations: 426.625 KiB)
  0.002491 seconds (9.43 k allocations: 405.531 KiB)
  0.003230 seconds (9.22 k allocations: 495.453 KiB)
  0.002418 seconds (5.50 k allocations: 405.812 KiB)
  0.009391 seconds (36.68 k allocations: 1.940 MiB)
  0.425004 seconds (825.80 k allocations: 41.259 MiB, 89.10% compilation ti
me)
(n, ts) = (2, [0.00190136, 0.002473006, 0.002361067, 0.002976003, 0.0020517
98, 0.00885639])
  0.004919 seconds (11.21 k allocations: 475.281 KiB)
  0.005975 seconds (11.07 k allocations: 617.133 KiB)
  0.005780 seconds (15.21 k allocations: 618.844 KiB)
  0.006905 seconds (14.03 k allocations: 732.977 KiB)
  0.004777 seconds (7.43 k allocations: 609.484 KiB)
  0.023603 seconds (64.76 k allocations: 3.796 MiB)
  0.103302 seconds (248.24 k allocations: 13.620 MiB)
(n, ts) = (3, [0.004771473, 0.006042376, 0.005618909, 0.006840671, 0.004423
205, 0.022798981])
  0.007791 seconds (13.55 k allocations: 577.594 KiB)
  0.008198 seconds (12.57 k allocations: 714.039 KiB)
  0.007420 seconds (16.26 k allocations: 667.219 KiB)
  0.008412 seconds (14.58 k allocations: 793.398 KiB)
  0.006378 seconds (7.54 k allocations: 737.359 KiB)
  0.032961 seconds (80.79 k allocations: 5.176 MiB)
  0.141721 seconds (291.40 k allocations: 17.231 MiB)
(n, ts) = (4, [0.007413577, 0.008128144, 0.007284139, 0.008363003, 0.006027
336, 0.032710475])
  0.012770 seconds (15.82 k allocations: 673.938 KiB)
  0.013563 seconds (14.85 k allocations: 873.758 KiB)
  0.010791 seconds (18.01 k allocations: 728.906 KiB)
  0.011868 seconds (15.90 k allocations: 900.789 KiB)
  0.009919 seconds (7.63 k allocations: 942.953 KiB)
  0.050161 seconds (94.58 k allocations: 6.476 MiB)
  0.217045 seconds (334.40 k allocations: 21.064 MiB)
(n, ts) = (5, [0.012403469, 0.013349415, 0.010666249, 0.011785663, 0.009517
296, 0.049586869])
  0.022508 seconds (19.61 k allocations: 842.422 KiB)
  0.022030 seconds (17.33 k allocations: 1.009 MiB)
  0.017098 seconds (21.76 k allocations: 869.609 KiB)
  0.017359 seconds (18.18 k allocations: 1.051 MiB)
  0.015266 seconds (7.83 k allocations: 1.236 MiB)
  0.137513 seconds (114.08 k allocations: 7.992 MiB, 37.42% gc time)
  0.404357 seconds (398.44 k allocations: 25.985 MiB, 12.72% gc time)
(n, ts) = (6, [0.021957796, 0.021963546, 0.016848255, 0.01761928, 0.0146383
37, 0.078780544])
  0.033661 seconds (22.45 k allocations: 964.359 KiB)
  0.033545 seconds (20.42 k allocations: 1.200 MiB)
  0.020513 seconds (22.60 k allocations: 906.984 KiB)
  0.021613 seconds (19.18 k allocations: 1.179 MiB)
  0.025723 seconds (8.25 k allocations: 1.701 MiB)
  0.103495 seconds (134.25 k allocations: 9.484 MiB)
  0.476693 seconds (455.14 k allocations: 30.847 MiB)
(n, ts) = (7, [0.033321011, 0.033652769, 0.020411305, 0.021500198, 0.025191
728, 0.103290205])
  0.056500 seconds (28.34 k allocations: 1.193 MiB)
  0.052646 seconds (24.20 k allocations: 1.461 MiB)
  0.033317 seconds (35.38 k allocations: 1.257 MiB)
  0.028960 seconds (23.19 k allocations: 1.527 MiB)
  0.034952 seconds (8.07 k allocations: 2.166 MiB)
  0.140714 seconds (161.03 k allocations: 12.050 MiB)
  0.692790 seconds (561.24 k allocations: 39.373 MiB)
(n, ts) = (8, [0.056044153, 0.052248944, 0.033190482, 0.029023986, 0.034630
534, 0.139681119])
  0.091963 seconds (34.60 k allocations: 1.442 MiB)
  0.085079 seconds (29.31 k allocations: 1.761 MiB)
  0.044822 seconds (36.66 k allocations: 1.325 MiB)
  0.039760 seconds (24.50 k allocations: 1.751 MiB)
  0.065133 seconds (8.32 k allocations: 2.881 MiB)
  0.255534 seconds (182.26 k allocations: 14.536 MiB, 24.70% gc time)
  1.098378 seconds (632.14 k allocations: 47.458 MiB, 5.75% gc time)
(n, ts) = (9, [0.093658169, 0.085309577, 0.044281349, 0.039251508, 0.065023
212, 0.187508358])
  0.141230 seconds (39.46 k allocations: 1.668 MiB)
  0.126528 seconds (32.80 k allocations: 1.974 MiB)
  0.052832 seconds (32.61 k allocations: 1.296 MiB)
  0.046715 seconds (24.70 k allocations: 1.831 MiB)
  0.086122 seconds (8.44 k allocations: 3.785 MiB)
  0.244318 seconds (214.93 k allocations: 17.122 MiB)
  1.400269 seconds (706.70 k allocations: 55.415 MiB)
(n, ts) = (10, [0.14109343, 0.125941926, 0.053986025, 0.047086284, 0.085847
463, 0.247554958])
  0.329792 seconds (53.48 k allocations: 2.260 MiB)
  0.286264 seconds (43.99 k allocations: 2.705 MiB)
  0.108849 seconds (50.52 k allocations: 1.957 MiB)
  0.106015 seconds (38.95 k allocations: 2.892 MiB)
  0.249819 seconds (9.66 k allocations: 6.391 MiB)
  0.375849 seconds (265.69 k allocations: 22.871 MiB)
  2.984550 seconds (925.41 k allocations: 78.219 MiB, 2.34% gc time)
(n, ts) = (12, [0.328819947, 0.288407606, 0.108140057, 0.138838114, 0.25096
1478, 0.41138527])
  1.040379 seconds (78.77 k allocations: 3.339 MiB)
  0.892269 seconds (63.86 k allocations: 4.027 MiB)
  0.191343 seconds (61.42 k allocations: 2.427 MiB)
  0.183842 seconds (48.43 k allocations: 4.218 MiB)
  0.740938 seconds (10.67 k allocations: 13.194 MiB)
  0.769449 seconds (385.68 k allocations: 36.273 MiB, 13.26% gc time)
  7.555266 seconds (1.30 M allocations: 127.022 MiB, 1.86% gc time)
(n, ts) = (15, [1.036677607, 0.893899475, 0.191667454, 0.183949688, 0.76903
1333, 0.659933841])
  1.889528 seconds (98.57 k allocations: 4.185 MiB)
  1.640299 seconds (79.31 k allocations: 4.982 MiB)
  0.246868 seconds (61.78 k allocations: 2.495 MiB)
  0.450026 seconds (48.09 k allocations: 4.930 MiB, 47.34% gc time)
  1.459557 seconds (11.41 k allocations: 20.257 MiB)
  0.948305 seconds (446.16 k allocations: 46.481 MiB)
 13.168929 seconds (1.49 M allocations: 166.726 MiB, 1.88% gc time)
(n, ts) = (17, [1.891517543, 1.641265601, 0.246893831, 0.2313017, 1.5473722
03, 0.973605524])
12-element Vector{Vector{Float64}}:
 [0.00190136, 0.002473006, 0.002361067, 0.002976003, 0.002051798, 0.0088563
9]
 [0.004771473, 0.006042376, 0.005618909, 0.006840671, 0.004423205, 0.022798
981]
 [0.007413577, 0.008128144, 0.007284139, 0.008363003, 0.006027336, 0.032710
475]
 [0.012403469, 0.013349415, 0.010666249, 0.011785663, 0.009517296, 0.049586
869]
 [0.021957796, 0.021963546, 0.016848255, 0.01761928, 0.014638337, 0.0787805
44]
 [0.033321011, 0.033652769, 0.020411305, 0.021500198, 0.025191728, 0.103290
205]
 [0.056044153, 0.052248944, 0.033190482, 0.029023986, 0.034630534, 0.139681
119]
 [0.093658169, 0.085309577, 0.044281349, 0.039251508, 0.065023212, 0.187508
358]
 [0.14109343, 0.125941926, 0.053986025, 0.047086284, 0.085847463, 0.2475549
58]
 [0.328819947, 0.288407606, 0.108140057, 0.138838114, 0.250961478, 0.411385
27]
 [1.036677607, 0.893899475, 0.191667454, 0.183949688, 0.769031333, 0.659933
841]
 [1.891517543, 1.641265601, 0.246893831, 0.2313017, 1.547372203, 0.97360552
4]
csadata_same_solver = [[csa_same_solver[j][i] for j in eachindex(csa_same_solver)]
                       for i in eachindex(csa_same_solver[1])]

plt_same = plot(title = "Brusselator same-solver adjoint comparison (EnzymeVJP)");
for (i, (name, _, _, _)) in enumerate(same_solver_configs)
    plot!(plt_same, n_to_param.(csan), csadata_same_solver[i], lab = name,
        lw = lw, marksize = ms, linestyle = :auto, marker = :auto)
end
xaxis!(plt_same, "Number of Parameters", :log10);
yaxis!(plt_same, "Runtime (s)", :log10);
plot!(plt_same, legend = :outertopleft, size = (1200, 600))

Gradient accuracy scaling

The Brusselator's stiffness grows with the grid refinement (the diffusion rate scales as 1/dx^2 ∝ N^2), so sweeping the size also sweeps the stiffness. This measures the gradient accuracy of each adjoint configuration at the benchmark tolerances against a tight (abstol = reltol = 1e-12) reference. The reference is cross-validated at every size by computing it with two completely independent implementations (the CVODES C adjoint and the native GaussAdjoint); the printed agreement is their relative difference, which was at most ~3e-10 at every size in the latest run.

accuracy_configs = vcat(
    same_solver_configs,
    [("Rodas5 + GaussAdjoint", Rodas5(autodiff = AutoFiniteDiff()), gauss_enz, false)],
)

relerr(a, b) = norm(vec(a) .- vec(b)) / norm(vec(b))

csa_accuracy = map(csan) do n
    bfun, b_u0, b_p, brusselator_jac, brusselator_comp = makebrusselator!(PROBS, n)
    prob = ODEProblem{true, SciMLBase.FullSpecialize}(bfun, b_u0, tspan, b_p)
    dg(out, u, p, t, i) = (out .= u .- 1.0)

    ref_sol = solve(prob, CVODE_BDF(linear_solver = :GMRES),
        abstol = 1e-12, reltol = 1e-12, dense = true)
    du0_ref, dp_ref = adjoint_sensitivities(
        ref_sol, CVODE_BDF(linear_solver = :GMRES); t = collect(bt),
        dgdu_discrete = dg, abstol = 1e-12, reltol = 1e-12, sensealg = gauss_enz)

    ref_sol2 = solve(prob, CVODE_BDF(linear_solver = :GMRES),
        abstol = 1e-12, reltol = 1e-12, saveat = bt)
    du0_ref2, dp_ref2 = adjoint_sensitivities(
        ref_sol2, CVODE_BDF(linear_solver = :GMRES); t = collect(bt),
        dgdu_discrete = dg, abstol = 1e-12, reltol = 1e-12,
        sensealg = SundialsAdjoint(autojacvec = EnzymeVJP()))
    agreement = max(relerr(dp_ref2, dp_ref), relerr(du0_ref2, du0_ref))

    errs = map(accuracy_configs) do (name, alg, sensealg, dense_forward)
        du0, dp = run_same_solver(bfun, b_u0, b_p, alg, sensealg, dense_forward)
        relerr(dp, dp_ref)
    end
    @show n, agreement, errs
    errs
end
(n, agreement, errs) = (2, 2.1963974702304678e-11, [1.1360036542730342e-5, 
1.8833104332738728e-5, 1.0803476551076878e-5, 1.8894905218340053e-5, 2.3286
287656629697e-5, 2.3878725506819976e-5, 3.0394201713227987e-9])
(n, agreement, errs) = (3, 3.010588434012082e-10, [0.00028737837337994515, 
2.80086857116937e-5, 0.00023912428713242524, 3.714568563653246e-5, 0.001505
4829082979052, 0.0014766321747167997, 4.117603080409098e-5])
(n, agreement, errs) = (4, 1.947292936038399e-10, [0.0001687926610456094, 1
.874971527376774e-5, 0.00015883934120616695, 3.691915020869242e-5, 0.001018
0981813033077, 0.0010587913364942991, 1.581329703357047e-5])
(n, agreement, errs) = (5, 1.1458770290019084e-10, [0.00012861447619119944,
 2.1358797849169975e-5, 0.000139683208049175, 1.55745766292639e-5, 0.000778
2160218531484, 0.0007979527317029229, 1.132856087354418e-5])
(n, agreement, errs) = (6, 1.1084599668766109e-10, [0.00010956844524867612,
 2.0269472136113163e-5, 0.00010943230107359006, 2.3143298643346274e-5, 0.00
06166101072026013, 0.0006496243246276735, 8.017922704767997e-6])
(n, agreement, errs) = (7, 9.654055320116652e-11, [9.417453384377807e-5, 1.
6873080971675968e-5, 9.234561714067285e-5, 1.3570524356883525e-5, 0.0005370
742416252596, 0.0005447209207241839, 5.06464820260714e-6])
(n, agreement, errs) = (8, 9.794632556352656e-11, [8.831634072154492e-5, 1.
6851933957859684e-5, 8.557186511651593e-5, 9.322418419792044e-6, 0.00055794
43464255073, 0.0005556296118236512, 3.5306893395213555e-6])
(n, agreement, errs) = (9, 8.892003424263953e-11, [8.565053270513868e-5, 1.
462096461377351e-5, 8.083127239464406e-5, 4.795229093527008e-6, 0.000464146
1425833273, 0.0004611103756681643, 2.9710570433722095e-6])
(n, agreement, errs) = (10, 9.066173619606881e-11, [7.682476620291043e-5, 1
.8437614241217117e-5, 8.89644356071134e-5, 3.104678797113946e-5, 0.00047807
01676914919, 0.0004783628725488691, 2.794065437343273e-6])
(n, agreement, errs) = (12, 7.67140688387645e-11, [7.411141828954412e-5, 2.
2215058788081397e-5, 6.526502659179764e-5, 1.2430558885657634e-5, 0.0004238
766249800418, 0.00042744908731910143, 2.730558096856138e-6])
(n, agreement, errs) = (15, 6.721550805195372e-11, [7.380143188028393e-5, 9
.487618540479155e-6, 9.757134118760027e-5, 4.297504300599217e-5, 0.00044138
559604753525, 0.0004439295999529228, 2.9139178684502948e-6])
(n, agreement, errs) = (17, 7.560204822537469e-11, [6.830653284894397e-5, 1
.0787115927338272e-5, 0.0001115941045792728, 7.272612520132917e-5, 0.000414
3204376360607, 0.00041593354967809336, 3.89049367453392e-6])
12-element Vector{Vector{Float64}}:
 [1.1360036542730342e-5, 1.8833104332738728e-5, 1.0803476551076878e-5, 1.88
94905218340053e-5, 2.3286287656629697e-5, 2.3878725506819976e-5, 3.03942017
13227987e-9]
 [0.00028737837337994515, 2.80086857116937e-5, 0.00023912428713242524, 3.71
4568563653246e-5, 0.0015054829082979052, 0.0014766321747167997, 4.117603080
409098e-5]
 [0.0001687926610456094, 1.874971527376774e-5, 0.00015883934120616695, 3.69
1915020869242e-5, 0.0010180981813033077, 0.0010587913364942991, 1.581329703
357047e-5]
 [0.00012861447619119944, 2.1358797849169975e-5, 0.000139683208049175, 1.55
745766292639e-5, 0.0007782160218531484, 0.0007979527317029229, 1.1328560873
54418e-5]
 [0.00010956844524867612, 2.0269472136113163e-5, 0.00010943230107359006, 2.
3143298643346274e-5, 0.0006166101072026013, 0.0006496243246276735, 8.017922
704767997e-6]
 [9.417453384377807e-5, 1.6873080971675968e-5, 9.234561714067285e-5, 1.3570
524356883525e-5, 0.0005370742416252596, 0.0005447209207241839, 5.0646482026
0714e-6]
 [8.831634072154492e-5, 1.6851933957859684e-5, 8.557186511651593e-5, 9.3224
18419792044e-6, 0.0005579443464255073, 0.0005556296118236512, 3.53068933952
13555e-6]
 [8.565053270513868e-5, 1.462096461377351e-5, 8.083127239464406e-5, 4.79522
9093527008e-6, 0.0004641461425833273, 0.0004611103756681643, 2.971057043372
2095e-6]
 [7.682476620291043e-5, 1.8437614241217117e-5, 8.89644356071134e-5, 3.10467
8797113946e-5, 0.0004780701676914919, 0.0004783628725488691, 2.794065437343
273e-6]
 [7.411141828954412e-5, 2.2215058788081397e-5, 6.526502659179764e-5, 1.2430
558885657634e-5, 0.0004238766249800418, 0.00042744908731910143, 2.730558096
856138e-6]
 [7.380143188028393e-5, 9.487618540479155e-6, 9.757134118760027e-5, 4.29750
4300599217e-5, 0.00044138559604753525, 0.0004439295999529228, 2.91391786845
02948e-6]
 [6.830653284894397e-5, 1.0787115927338272e-5, 0.0001115941045792728, 7.272
612520132917e-5, 0.0004143204376360607, 0.00041593354967809336, 3.890493674
53392e-6]
csadata_accuracy = [[csa_accuracy[j][i] for j in eachindex(csa_accuracy)]
                    for i in eachindex(csa_accuracy[1])]

plt_acc = plot(title = "Brusselator adjoint gradient accuracy scaling");
for (i, (name, _, _, _)) in enumerate(accuracy_configs)
    plot!(plt_acc, n_to_param.(csan), csadata_accuracy[i], lab = name,
        lw = lw, marksize = ms, linestyle = :auto, marker = :auto)
end
xaxis!(plt_acc, "Number of Parameters", :log10);
yaxis!(plt_acc, "Relative L2 error of dG/dp", :log10);
plot!(plt_acc, legend = :outertopleft, size = (1200, 600))

Peak Memory Benchmarks

Measures the memory consumed by each sensitivity computation. Each configuration runs in a separate subprocess. We measure current RSS (from /proc/self/statm) before and after the computation to isolate the memory used by the sensitivity solve from the large fixed cost of package loading.

const CHILD_PREAMBLE = raw"""
using OrdinaryDiffEq, OrdinaryDiffEqRosenbrock, ReverseDiff, ForwardDiff, FiniteDiff,
      SciMLSensitivity
using LinearAlgebra, Mooncake

function get_rss_mib()
    statm = read("/proc/self/statm", String)
    resident_pages = parse(Int, split(statm)[2])
    return resident_pages * 4096 / (1024^2)
end

function makebrusselator(N = 8)
    xyd_brusselator = range(0, stop = 1, length = N)
    function limit(a, N)
        if a == N+1
            return 1
        elseif a == 0
            return N
        else
            return a
        end
    end
    brusselator_f(x, y, t) = ifelse(
        (((x-0.3)^2 + (y-0.6)^2) <= 0.1^2) &&
        (t >= 1.1), 5.0, 0.0)
    brusselator_2d_loop = let N=N, xyd=xyd_brusselator, dx=step(xyd_brusselator)
        function brusselator_2d_loop(du, u, p, t)
            @inbounds begin
                ii1 = N^2
                ii2 = ii1+N^2
                ii3 = ii2+2(N^2)
                A = @view p[1:ii1]
                B = @view p[(ii1 + 1):ii2]
                α = @view p[(ii2 + 1):ii3]
                II = LinearIndices((N, N, 2))
                for I in CartesianIndices((N, N))
                    x = xyd[I[1]]
                    y = xyd[I[2]]
                    i = I[1]
                    j = I[2]
                    ip1 = limit(i+1, N);
                    im1 = limit(i-1, N)
                    jp1 = limit(j+1, N);
                    jm1 = limit(j-1, N)
                    du[II[i, j, 1]] = α[II[
                                          i, j, 1]]*(u[II[im1, j, 1]] + u[II[ip1, j, 1]] +
                                                     u[II[i, jp1, 1]] + u[II[i, jm1, 1]] -
                                                     4u[II[i, j, 1]])/dx^2 +
                                      B[II[i, j, 1]] + u[II[i, j, 1]]^2*u[II[i, j, 2]] -
                                      (A[II[i, j, 1]] + 1)*u[II[i, j, 1]] +
                                      brusselator_f(x, y, t)
                end
                for I in CartesianIndices((N, N))
                    i = I[1]
                    j = I[2]
                    ip1 = limit(i+1, N)
                    im1 = limit(i-1, N)
                    jp1 = limit(j+1, N)
                    jm1 = limit(j-1, N)
                    du[II[i, j, 2]] = α[II[
                        i, j, 2]]*(u[II[im1, j, 2]] + u[II[ip1, j, 2]] + u[II[i, jp1, 2]] +
                                   u[II[i, jm1, 2]] - 4u[II[i, j, 2]])/dx^2 +
                                      A[II[i, j, 1]]*u[II[i, j, 1]] -
                                      u[II[i, j, 1]]^2*u[II[i, j, 2]]
                end
                return nothing
            end
        end
    end
    function init_brusselator_2d(xyd)
        N = length(xyd)
        u = zeros(N, N, 2)
        for I in CartesianIndices((N, N))
            x = xyd[I[1]]
            y = xyd[I[2]]
            u[I, 1] = 22*(y*(1-y))^(3/2)
            u[I, 2] = 27*(x*(1-x))^(3/2)
        end
        vec(u)
    end
    dx = step(xyd_brusselator)
    e1 = ones(N-1)
    off = N-1
    e4 = ones(N-off)
    T = diagm(0=>-2ones(N), -1=>e1, 1=>e1, off=>e4, -off=>e4) ./ dx^2
    Ie = Matrix{Float64}(I, N, N)
    Op = kron(Ie, T) + kron(T, Ie)
    brusselator_jac = let N=N
        (J, a, p, t) -> begin
            ii1 = N^2
            ii2 = ii1+N^2
            ii3 = ii2+2(N^2)
            A = @view p[1:ii1]
            B = @view p[(ii1 + 1):ii2]
            α = @view p[(ii2 + 1):ii3]
            u = @view a[1:(end ÷ 2)]
            v = @view a[(end ÷ 2 + 1):end]
            N2 = length(a)÷2
            α1 = @view α[1:(end ÷ 2)]
            α2 = @view α[(end ÷ 2 + 1):end]
            fill!(J, 0)
            J[1:N2, 1:N2] .= α1 .* Op
            J[(N2 + 1):end, (N2 + 1):end] .= α2 .* Op
            J1 = @view J[1:N2, 1:N2]
            J2 = @view J[(N2 + 1):end, 1:N2]
            J3 = @view J[1:N2, (N2 + 1):end]
            J4 = @view J[(N2 + 1):end, (N2 + 1):end]
            J1[diagind(J1)] .+= @. 2u*v-(A+1)
            J2[diagind(J2)] .= @. A-2u*v
            J3[diagind(J3)] .= @. u^2
            J4[diagind(J4)] .+= @. -u^2
            nothing
        end
    end
    u0 = init_brusselator_2d(xyd_brusselator)
    p = [fill(3.4, N^2); fill(1.0, N^2); fill(10.0, 2*N^2)]
    brusselator_2d_loop, u0, p, brusselator_jac
end

Base.vec(v::Adjoint{<:Real, <:AbstractVector}) = vec(v')

bt = 0:0.1:1
tspan = (0.0, 1.0)
tols = (abstol = 1e-5, reltol = 1e-7)

function auto_sen_l2(
        f, u0, tspan, p, t, alg = Tsit5(); diffalg = ReverseDiff.gradient, kwargs...)
    test_f(p) = begin
        prob = ODEProblem{true, SciMLBase.FullSpecialize}(f, convert.(eltype(p), u0), tspan, p)
        sol = solve(prob, alg, saveat = t; kwargs...)
        sum(sol.u) do x
            sum(z->(1-z)^2/2, x)
        end
    end
    diffalg(test_f, p)
end

@inline function diffeq_sen_l2(df, u0, tspan, p, t, alg = Tsit5();
        abstol = 1e-5, reltol = 1e-7, iabstol = abstol, ireltol = reltol,
        sensalg = SensitivityAlg(), kwargs...)
    prob = ODEProblem{true, SciMLBase.FullSpecialize}(df, u0, tspan, p)
    saveat = tspan[1] != t[1] && tspan[end] != t[end] ? vcat(tspan[1], t, tspan[end]) : t
    sol = solve(prob, alg, abstol = abstol, reltol = reltol, saveat = saveat; kwargs...)
    dg(out, u, p, t, i) = (out.=u .- 1.0)
    adjoint_sensitivities(sol, alg; t, abstol = abstol, dgdu_discrete = dg,
        reltol = reltol, sensealg = sensalg)
end
"""

const PROJECT_DIR = @__DIR__

function run_memory_benchmark(n::Int, method_setup::String)
    child_script = CHILD_PREAMBLE * """

    n = $(n)
    bfun, b_u0, b_p, brusselator_jac = makebrusselator(n)

    GC.gc(); GC.gc()
    rss_before = get_rss_mib()

    """ * method_setup * """

    GC.gc(); GC.gc()
    rss_after = get_rss_mib()

    println("BRUSSMEM_TIMING:", t)
    println("BRUSSMEM_RSS_BEFORE:", rss_before)
    println("BRUSSMEM_RSS_AFTER:", rss_after)
    """

    try
        output = read(
            `$(Base.julia_cmd()) --project=$(PROJECT_DIR) -e $(child_script)`, String)
        time_m = match(r"BRUSSMEM_TIMING:([\d.eE+-]+)", output)
        before_m = match(r"BRUSSMEM_RSS_BEFORE:([\d.eE+-]+)", output)
        after_m = match(r"BRUSSMEM_RSS_AFTER:([\d.eE+-]+)", output)
        if time_m === nothing || before_m === nothing || after_m === nothing
            @warn "Failed to parse subprocess output" n output
            return (; rss_before = NaN, rss_after = NaN, delta_mib = NaN, timing = NaN)
        end
        timing = parse(Float64, time_m.captures[1])
        rss_before = parse(Float64, before_m.captures[1])
        rss_after = parse(Float64, after_m.captures[1])
        delta_mib = rss_after - rss_before
        return (; rss_before, rss_after, delta_mib, timing)
    catch e
        @warn "Subprocess failed" n exception = (e, catch_backtrace())
        return (; rss_before = NaN, rss_after = NaN, delta_mib = NaN, timing = NaN)
    end
end

mem_sizes = [2, 4, 6, 8, 10, 12]
6-element Vector{Int64}:
  2
  4
  6
  8
 10
 12
forwarddiff_mem = map(mem_sizes) do n
    result = run_memory_benchmark(n, """
    auto_sen_l2(bfun, b_u0, tspan, b_p, bt, Rodas5();
        diffalg = ForwardDiff.gradient, tols...)
    t = @elapsed auto_sen_l2(bfun, b_u0, tspan, b_p, bt, Rodas5();
        diffalg = ForwardDiff.gradient, tols...)
    """)
    @show n, result
    result
end
(n, result) = (2, (rss_before = 1057.6015625, rss_after = 1211.625, delta_m
ib = 154.0234375, timing = 0.001595871))
(n, result) = (4, (rss_before = 1047.94140625, rss_after = 1251.41015625, d
elta_mib = 203.46875, timing = 0.040372941))
(n, result) = (6, (rss_before = 1076.875, rss_after = 1230.28125, delta_mib
 = 153.40625, timing = 0.449492779))
(n, result) = (8, (rss_before = 1061.01953125, rss_after = 1230.3125, delta
_mib = 169.29296875, timing = 2.332212718))
(n, result) = (10, (rss_before = 1055.66796875, rss_after = 1236.26953125, 
delta_mib = 180.6015625, timing = 9.810840266))
(n, result) = (12, (rss_before = 1048.5390625, rss_after = 1261.6484375, de
lta_mib = 213.109375, timing = 31.595231227))
6-element Vector{@NamedTuple{rss_before::Float64, rss_after::Float64, delta
_mib::Float64, timing::Float64}}:
 (rss_before = 1057.6015625, rss_after = 1211.625, delta_mib = 154.0234375,
 timing = 0.001595871)
 (rss_before = 1047.94140625, rss_after = 1251.41015625, delta_mib = 203.46
875, timing = 0.040372941)
 (rss_before = 1076.875, rss_after = 1230.28125, delta_mib = 153.40625, tim
ing = 0.449492779)
 (rss_before = 1061.01953125, rss_after = 1230.3125, delta_mib = 169.292968
75, timing = 2.332212718)
 (rss_before = 1055.66796875, rss_after = 1236.26953125, delta_mib = 180.60
15625, timing = 9.810840266)
 (rss_before = 1048.5390625, rss_after = 1261.6484375, delta_mib = 213.1093
75, timing = 31.595231227)
numdiff_mem = map(mem_sizes) do n
    result = run_memory_benchmark(n, """
    auto_sen_l2(bfun, b_u0, tspan, b_p, bt, Rodas5();
        diffalg = FiniteDiff.finite_difference_gradient, tols...)
    t = @elapsed auto_sen_l2(bfun, b_u0, tspan, b_p, bt, Rodas5();
        diffalg = FiniteDiff.finite_difference_gradient, tols...)
    """)
    @show n, result
    result
end
(n, result) = (2, (rss_before = 1054.09375, rss_after = 1168.4609375, delta
_mib = 114.3671875, timing = 0.004104796))
(n, result) = (4, (rss_before = 1050.1953125, rss_after = 1172.1640625, del
ta_mib = 121.96875, timing = 0.125163448))
(n, result) = (6, (rss_before = 1054.73046875, rss_after = 1185.44921875, d
elta_mib = 130.71875, timing = 0.827041386))
(n, result) = (8, (rss_before = 1058.1796875, rss_after = 1184.7421875, del
ta_mib = 126.5625, timing = 4.384329712))
(n, result) = (10, (rss_before = 1052.39453125, rss_after = 1166.7265625, d
elta_mib = 114.33203125, timing = 18.461356712))
(n, result) = (12, (rss_before = 1052.21484375, rss_after = 1187.4375, delt
a_mib = 135.22265625, timing = 87.79551839))
6-element Vector{@NamedTuple{rss_before::Float64, rss_after::Float64, delta
_mib::Float64, timing::Float64}}:
 (rss_before = 1054.09375, rss_after = 1168.4609375, delta_mib = 114.367187
5, timing = 0.004104796)
 (rss_before = 1050.1953125, rss_after = 1172.1640625, delta_mib = 121.9687
5, timing = 0.125163448)
 (rss_before = 1054.73046875, rss_after = 1185.44921875, delta_mib = 130.71
875, timing = 0.827041386)
 (rss_before = 1058.1796875, rss_after = 1184.7421875, delta_mib = 126.5625
, timing = 4.384329712)
 (rss_before = 1052.39453125, rss_after = 1166.7265625, delta_mib = 114.332
03125, timing = 18.461356712)
 (rss_before = 1052.21484375, rss_after = 1187.4375, delta_mib = 135.222656
25, timing = 87.79551839)
adjoint_ad_configs = [
    ("Interp user-Jacobian",
     "InterpolatingAdjoint(autodiff = false, autojacvec = false)", true),
    ("Interp AD-Jacobian",
     "InterpolatingAdjoint(autodiff = true, autojacvec = false)", false),
    ("Quad user-Jacobian",
     "QuadratureAdjoint(autodiff = false, autojacvec = false)", true),
    ("Quad AD-Jacobian",
     "QuadratureAdjoint(autodiff = true, autojacvec = false)", false),
    ("Gauss AD-Jacobian",
     "GaussAdjoint(autodiff = true, autojacvec = false)", false),
    ("GaussKronrod AD-Jacobian",
     "GaussKronrodAdjoint(autodiff = true, autojacvec = false)", false),
]

adjoint_ad_mem = map(adjoint_ad_configs) do (name, sensalg_str, needs_jac)
    results = map(mem_sizes) do n
        f_expr = needs_jac ? "ODEFunction(bfun, jac = brusselator_jac)" : "bfun"
        result = run_memory_benchmark(n, """
        sensalg = $(sensalg_str)
        f = $(f_expr)
        solver = Rodas5(autodiff = AutoFiniteDiff())
        diffeq_sen_l2(f, b_u0, tspan, b_p, bt, solver; sensalg = sensalg, tols...)
        t = @elapsed diffeq_sen_l2(f, b_u0, tspan, b_p, bt, solver;
            sensalg = sensalg, tols...)
        """)
        @show name, n, result
        result
    end
    (name = name, results = results)
end
(name, n, result) = ("Interp user-Jacobian", 2, (rss_before = 1055.8515625,
 rss_after = 1243.625, delta_mib = 187.7734375, timing = 0.006631972))
(name, n, result) = ("Interp user-Jacobian", 4, (rss_before = 1047.72265625
, rss_after = 1273.60546875, delta_mib = 225.8828125, timing = 0.110097776)
)
(name, n, result) = ("Interp user-Jacobian", 6, (rss_before = 1068.40625, r
ss_after = 1234.95703125, delta_mib = 166.55078125, timing = 1.008895529))
(name, n, result) = ("Interp user-Jacobian", 8, (rss_before = 1058.828125, 
rss_after = 1247.5703125, delta_mib = 188.7421875, timing = 5.846499621))
(name, n, result) = ("Interp user-Jacobian", 10, (rss_before = 1055.7382812
5, rss_after = 1246.3203125, delta_mib = 190.58203125, timing = 22.66086248
4))
(name, n, result) = ("Interp user-Jacobian", 12, (rss_before = 1056.0429687
5, rss_after = 1263.84375, delta_mib = 207.80078125, timing = 66.422538988)
)
(name, n, result) = ("Interp AD-Jacobian", 2, (rss_before = 1055.75, rss_af
ter = 1263.72265625, delta_mib = 207.97265625, timing = 0.003739789))
(name, n, result) = ("Interp AD-Jacobian", 4, (rss_before = 1051.2265625, r
ss_after = 1252.8125, delta_mib = 201.5859375, timing = 0.174536107))
(name, n, result) = ("Interp AD-Jacobian", 6, (rss_before = 1049.78515625, 
rss_after = 1288.4765625, delta_mib = 238.69140625, timing = 0.392971194))
(name, n, result) = ("Interp AD-Jacobian", 8, (rss_before = 1056.34765625, 
rss_after = 1276.18359375, delta_mib = 219.8359375, timing = 2.362474285))
(name, n, result) = ("Interp AD-Jacobian", 10, (rss_before = 1052.0078125, 
rss_after = 1286.21484375, delta_mib = 234.20703125, timing = 9.83631856))
(name, n, result) = ("Interp AD-Jacobian", 12, (rss_before = 1056.4140625, 
rss_after = 1278.4375, delta_mib = 222.0234375, timing = 29.941231494))
(name, n, result) = ("Quad user-Jacobian", 2, (rss_before = 1054.71875, rss
_after = 1289.71875, delta_mib = 235.0, timing = 0.002546275))
(name, n, result) = ("Quad user-Jacobian", 4, (rss_before = 1052.92578125, 
rss_after = 1300.359375, delta_mib = 247.43359375, timing = 0.010020213))
(name, n, result) = ("Quad user-Jacobian", 6, (rss_before = 1053.40625, rss
_after = 1310.44140625, delta_mib = 257.03515625, timing = 0.040376269))
(name, n, result) = ("Quad user-Jacobian", 8, (rss_before = 1054.55859375, 
rss_after = 1316.546875, delta_mib = 261.98828125, timing = 0.126771957))
(name, n, result) = ("Quad user-Jacobian", 10, (rss_before = 1053.41796875,
 rss_after = 1298.0703125, delta_mib = 244.65234375, timing = 0.382198548))
(name, n, result) = ("Quad user-Jacobian", 12, (rss_before = 1058.3203125, 
rss_after = 1301.08984375, delta_mib = 242.76953125, timing = 0.900633249))
(name, n, result) = ("Quad AD-Jacobian", 2, (rss_before = 1047.046875, rss_
after = 1318.94921875, delta_mib = 271.90234375, timing = 0.002255747))
(name, n, result) = ("Quad AD-Jacobian", 4, (rss_before = 1053.12890625, rs
s_after = 1314.1640625, delta_mib = 261.03515625, timing = 0.011100996))
(name, n, result) = ("Quad AD-Jacobian", 6, (rss_before = 1053.76171875, rs
s_after = 1336.9140625, delta_mib = 283.15234375, timing = 0.069339313))
(name, n, result) = ("Quad AD-Jacobian", 8, (rss_before = 1054.44140625, rs
s_after = 1315.26953125, delta_mib = 260.828125, timing = 0.347147345))
(name, n, result) = ("Quad AD-Jacobian", 10, (rss_before = 1057.5703125, rs
s_after = 1359.48828125, delta_mib = 301.91796875, timing = 1.281432732))
(name, n, result) = ("Quad AD-Jacobian", 12, (rss_before = 1053.83984375, r
ss_after = 1301.41015625, delta_mib = 247.5703125, timing = 3.96847054))
(name, n, result) = ("Gauss AD-Jacobian", 2, (rss_before = 1056.65625, rss_
after = 1270.53125, delta_mib = 213.875, timing = 0.002586566))
(name, n, result) = ("Gauss AD-Jacobian", 4, (rss_before = 1058.921875, rss
_after = 1269.734375, delta_mib = 210.8125, timing = 0.01208945))
(name, n, result) = ("Gauss AD-Jacobian", 6, (rss_before = 1056.6640625, rs
s_after = 1266.36328125, delta_mib = 209.69921875, timing = 0.0806165))
(name, n, result) = ("Gauss AD-Jacobian", 8, (rss_before = 1050.1640625, rs
s_after = 1271.05859375, delta_mib = 220.89453125, timing = 0.513967247))
(name, n, result) = ("Gauss AD-Jacobian", 10, (rss_before = 1058.91796875, 
rss_after = 1253.1484375, delta_mib = 194.23046875, timing = 1.523536295))
(name, n, result) = ("Gauss AD-Jacobian", 12, (rss_before = 1052.25390625, 
rss_after = 1281.953125, delta_mib = 229.69921875, timing = 4.669176973))
(name, n, result) = ("GaussKronrod AD-Jacobian", 2, (rss_before = 1058.8593
75, rss_after = 1261.9140625, delta_mib = 203.0546875, timing = 0.002911303
))
(name, n, result) = ("GaussKronrod AD-Jacobian", 4, (rss_before = 1060.1718
75, rss_after = 1275.9609375, delta_mib = 215.7890625, timing = 0.017846398
))
(name, n, result) = ("GaussKronrod AD-Jacobian", 6, (rss_before = 1056.1210
9375, rss_after = 1284.578125, delta_mib = 228.45703125, timing = 0.1116517
32))
(name, n, result) = ("GaussKronrod AD-Jacobian", 8, (rss_before = 1054.8984
375, rss_after = 1269.89453125, delta_mib = 214.99609375, timing = 0.515560
593))
(name, n, result) = ("GaussKronrod AD-Jacobian", 10, (rss_before = 1063.375
, rss_after = 1282.58203125, delta_mib = 219.20703125, timing = 1.943040913
))
(name, n, result) = ("GaussKronrod AD-Jacobian", 12, (rss_before = 1065.746
09375, rss_after = 1261.28125, delta_mib = 195.53515625, timing = 5.9435545
44))
6-element Vector{@NamedTuple{name::String, results::Vector{@NamedTuple{rss_
before::Float64, rss_after::Float64, delta_mib::Float64, timing::Float64}}}
}:
 (name = "Interp user-Jacobian", results = [(rss_before = 1055.8515625, rss
_after = 1243.625, delta_mib = 187.7734375, timing = 0.006631972), (rss_bef
ore = 1047.72265625, rss_after = 1273.60546875, delta_mib = 225.8828125, ti
ming = 0.110097776), (rss_before = 1068.40625, rss_after = 1234.95703125, d
elta_mib = 166.55078125, timing = 1.008895529), (rss_before = 1058.828125, 
rss_after = 1247.5703125, delta_mib = 188.7421875, timing = 5.846499621), (
rss_before = 1055.73828125, rss_after = 1246.3203125, delta_mib = 190.58203
125, timing = 22.660862484), (rss_before = 1056.04296875, rss_after = 1263.
84375, delta_mib = 207.80078125, timing = 66.422538988)])
 (name = "Interp AD-Jacobian", results = [(rss_before = 1055.75, rss_after 
= 1263.72265625, delta_mib = 207.97265625, timing = 0.003739789), (rss_befo
re = 1051.2265625, rss_after = 1252.8125, delta_mib = 201.5859375, timing =
 0.174536107), (rss_before = 1049.78515625, rss_after = 1288.4765625, delta
_mib = 238.69140625, timing = 0.392971194), (rss_before = 1056.34765625, rs
s_after = 1276.18359375, delta_mib = 219.8359375, timing = 2.362474285), (r
ss_before = 1052.0078125, rss_after = 1286.21484375, delta_mib = 234.207031
25, timing = 9.83631856), (rss_before = 1056.4140625, rss_after = 1278.4375
, delta_mib = 222.0234375, timing = 29.941231494)])
 (name = "Quad user-Jacobian", results = [(rss_before = 1054.71875, rss_aft
er = 1289.71875, delta_mib = 235.0, timing = 0.002546275), (rss_before = 10
52.92578125, rss_after = 1300.359375, delta_mib = 247.43359375, timing = 0.
010020213), (rss_before = 1053.40625, rss_after = 1310.44140625, delta_mib 
= 257.03515625, timing = 0.040376269), (rss_before = 1054.55859375, rss_aft
er = 1316.546875, delta_mib = 261.98828125, timing = 0.126771957), (rss_bef
ore = 1053.41796875, rss_after = 1298.0703125, delta_mib = 244.65234375, ti
ming = 0.382198548), (rss_before = 1058.3203125, rss_after = 1301.08984375,
 delta_mib = 242.76953125, timing = 0.900633249)])
 (name = "Quad AD-Jacobian", results = [(rss_before = 1047.046875, rss_afte
r = 1318.94921875, delta_mib = 271.90234375, timing = 0.002255747), (rss_be
fore = 1053.12890625, rss_after = 1314.1640625, delta_mib = 261.03515625, t
iming = 0.011100996), (rss_before = 1053.76171875, rss_after = 1336.9140625
, delta_mib = 283.15234375, timing = 0.069339313), (rss_before = 1054.44140
625, rss_after = 1315.26953125, delta_mib = 260.828125, timing = 0.34714734
5), (rss_before = 1057.5703125, rss_after = 1359.48828125, delta_mib = 301.
91796875, timing = 1.281432732), (rss_before = 1053.83984375, rss_after = 1
301.41015625, delta_mib = 247.5703125, timing = 3.96847054)])
 (name = "Gauss AD-Jacobian", results = [(rss_before = 1056.65625, rss_afte
r = 1270.53125, delta_mib = 213.875, timing = 0.002586566), (rss_before = 1
058.921875, rss_after = 1269.734375, delta_mib = 210.8125, timing = 0.01208
945), (rss_before = 1056.6640625, rss_after = 1266.36328125, delta_mib = 20
9.69921875, timing = 0.0806165), (rss_before = 1050.1640625, rss_after = 12
71.05859375, delta_mib = 220.89453125, timing = 0.513967247), (rss_before =
 1058.91796875, rss_after = 1253.1484375, delta_mib = 194.23046875, timing 
= 1.523536295), (rss_before = 1052.25390625, rss_after = 1281.953125, delta
_mib = 229.69921875, timing = 4.669176973)])
 (name = "GaussKronrod AD-Jacobian", results = [(rss_before = 1058.859375, 
rss_after = 1261.9140625, delta_mib = 203.0546875, timing = 0.002911303), (
rss_before = 1060.171875, rss_after = 1275.9609375, delta_mib = 215.7890625
, timing = 0.017846398), (rss_before = 1056.12109375, rss_after = 1284.5781
25, delta_mib = 228.45703125, timing = 0.111651732), (rss_before = 1054.898
4375, rss_after = 1269.89453125, delta_mib = 214.99609375, timing = 0.51556
0593), (rss_before = 1063.375, rss_after = 1282.58203125, delta_mib = 219.2
0703125, timing = 1.943040913), (rss_before = 1065.74609375, rss_after = 12
61.28125, delta_mib = 195.53515625, timing = 5.943554544)])
mem_params = n_to_param.(mem_sizes)

plt_mem1 = plot(title = "Brusselator Sensitivity Memory Scaling");
plot!(plt_mem1, mem_params, [r.delta_mib for r in forwarddiff_mem],
    lab = "Forward-Mode DSAAD", lw = lw, marksize = ms,
    linestyle = :auto, marker = :auto);
plot!(plt_mem1, mem_params, [r.delta_mib for r in numdiff_mem],
    lab = "Numerical Differentiation", lw = lw, marksize = ms,
    linestyle = :auto, marker = :auto);
for entry in adjoint_ad_mem
    plot!(plt_mem1, mem_params, [r.delta_mib for r in entry.results],
        lab = entry.name, lw = lw, marksize = ms,
        linestyle = :auto, marker = :auto)
end
xaxis!(plt_mem1, "Number of Parameters", :log10);
yaxis!(plt_mem1, "Memory (MiB)");
plot!(plt_mem1, legend = :outertopleft, size = (1200, 600))

Appendix

Appendix

These benchmarks are a part of the SciMLBenchmarks.jl repository, found at: https://github.com/SciML/SciMLBenchmarks.jl. For more information on high-performance scientific machine learning, check out the SciML Open Source Software Organization https://sciml.ai.

To locally run this benchmark, do the following commands:

using SciMLBenchmarks
SciMLBenchmarks.weave_file("benchmarks/AutomaticDifferentiation","BrussScaling.jmd")

Computer Information:

Julia Version 1.12.7
Commit 6d172b025e4 (2026-08-15 08:05 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 128 × AMD EPYC 7502 32-Core Processor
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, znver2)
  GC: Built with stock GC
Threads: 128 default, 1 interactive, 128 GC (on 128 virtual cores)
Environment:
  JULIA_NUM_THREADS = auto

Package Information:

Status `~/github-runners/amdci3-1/_work/SciMLBenchmarks.jl/SciMLBenchmarks.jl/benchmarks/AutomaticDifferentiation/Project.toml`
  [6e4b80f9] BenchmarkTools v1.8.0
  [0ca39b1e] Chairmarks v1.3.1
  [a93c6f00] DataFrames v1.8.2
  [1313f7d8] DataFramesMeta v0.15.6
  [a0c0ee7d] DifferentiationInterface v0.7.21
  [a82114a7] DifferentiationInterfaceTest v0.11.0
⌃ [7da242da] Enzyme v0.13.203
  [6a86dc24] FiniteDiff v2.33.0
  [f6369f11] ForwardDiff v1.4.6
⌃ [7ed4a6bd] LinearSolve v5.17.3
⌃ [da2b9cff] Mooncake v0.5.56
  [1dea7af3] OrdinaryDiffEq v7.8.1
⌃ [43230ef6] OrdinaryDiffEqRosenbrock v2.7.3
  [65888b18] ParameterizedFunctions v5.27.0
  [91a5bcdd] Plots v1.41.7
  [08abe8d2] PrettyTables v3.4.8
⌃ [37e2e3b7] ReverseDiff v1.17.0
  [31c91b34] SciMLBenchmarks v0.2.1 [loaded: `/home/crackauc/github-runners/amdci3-1/_work/SciMLBenchmarks.jl/SciMLBenchmarks.jl/src/SciMLBenchmarks.jl` (v0.2.1) expected `/home/crackauc/.julia/packages/SciMLBenchmarks/ceJyd/src/SciMLBenchmarks.jl` (v0.2.1)]
⌃ [1ed8b502] SciMLSensitivity v7.119.7
⌃ [90137ffa] StaticArrays v1.9.20
  [c3572dad] Sundials v6.7.1
  [9f7883ad] Tracker v0.2.39
  [e88e6eb3] Zygote v0.7.13
  [37e2e46d] LinearAlgebra v1.12.0
  [d6f4376e] Markdown v1.11.0
  [de0858da] Printf v1.11.0
  [8dfed614] Test v1.11.0
Info Packages marked with ⌃ have new versions available and may be upgradable.

And the full manifest:

Status `~/github-runners/amdci3-1/_work/SciMLBenchmarks.jl/SciMLBenchmarks.jl/benchmarks/AutomaticDifferentiation/Manifest.toml`
  [47edcb42] ADTypes v1.24.0
  [14f7f29c] AMD v0.5.4
  [621f4979] AbstractFFTs v1.5.0
  [6e696c72] AbstractPlutoDingetjes v1.4.1
  [1520ce14] AbstractTrees v0.4.5
  [7d9f7c33] Accessors v0.1.45
⌃ [79e6a3ab] Adapt v4.7.0
  [66dad0bd] AliasTables v1.1.3
  [9b6a8646] AllocCheck v0.2.6
  [ec485272] ArnoldiMethod v0.4.0
⌃ [4fba245c] ArrayInterface v7.30.1
⌃ [4c555306] ArrayLayouts v1.12.2
  [a9b6321e] Atomix v1.2.1
⌃ [ab4f0b2a] BFloat16s v0.6.1
⌃ [aae01518] BandedMatrices v1.12.0
  [6e4b80f9] BenchmarkTools v1.8.0
  [e2ed5e7c] Bijections v0.2.2
  [b2a6c25c] BinaryHeaps v1.1.0
  [caf10ac8] BipartiteGraphs v0.1.14
  [8e7c35d0] BlockArrays v1.10.0
  [70df07ce] BracketingNonlinearSolve v1.12.7
  [fa961155] CEnum v0.5.0
  [8be319e6] Chain v1.0.0
  [082447d4] ChainRules v1.73.0
  [d360d2e6] ChainRulesCore v1.26.1
  [0ca39b1e] Chairmarks v1.3.1
  [35d6a980] ColorSchemes v3.31.0
⌃ [3da002f7] ColorTypes v0.12.1
  [c3611d14] ColorVectorSpace v0.11.0
  [5ae59095] Colors v0.13.1
⌅ [861a8166] Combinatorics v1.0.2
  [38540f10] CommonSolve v0.2.14
  [bbf7d656] CommonSubexpressions v0.3.1
  [f70d9fcc] CommonWorldInvalidations v1.2.2
  [34da2185] Compat v4.18.1
  [b152e2b5] CompositeTypes v0.1.4
  [a33af91c] CompositionsBase v0.1.2
  [2569d6c7] ConcreteStructs v0.2.8
  [187b0558] ConstructionBase v1.6.0
  [d38c429a] Contour v0.6.3
  [a8cc5b0e] Crayons v4.2.0
  [9a962f9c] DataAPI v1.16.0
  [a93c6f00] DataFrames v1.8.2
  [1313f7d8] DataFramesMeta v0.15.6
  [864edb3b] DataStructures v0.19.6
  [e2d170a0] DataValueInterfaces v1.0.0
  [8bb1440f] DelimitedFiles v1.9.1
⌃ [2b5f629d] DiffEqBase v7.21.1
  [459566f4] DiffEqCallbacks v4.19.4
  [77a26b50] DiffEqNoiseProcess v5.36.3
  [163ba53b] DiffResults v1.1.0
  [b552c78f] DiffRules v1.16.0
  [a0c0ee7d] DifferentiationInterface v0.7.21
  [a82114a7] DifferentiationInterfaceTest v0.11.0
⌃ [8d63f2c5] DispatchDoctor v0.4.28
  [31c24e10] Distributions v0.25.131
  [ffbed154] DocStringExtensions v0.9.5
⌃ [5b8099bc] DomainSets v0.8.1
  [7c1d4256] DynamicPolynomials v0.6.8
  [4e289a0a] EnumX v1.0.7
⌃ [7da242da] Enzyme v0.13.203
  [f151be2c] EnzymeCore v0.8.21
  [e2ba6199] ExprTools v0.1.11
  [55351af7] ExproniconLite v0.10.14
  [c87230d0] FFMPEG v0.4.5
  [7034ab61] FastBroadcast v1.4.0
  [9aa1b823] FastClosures v0.3.2
  [a4df4552] FastPower v1.5.0
  [1a297f60] FillArrays v1.17.0
⌃ [64ca27bc] FindFirstFunctions v3.2.1
  [6a86dc24] FiniteDiff v2.33.0
⌅ [53c48c17] FixedPointNumbers v0.8.6
  [1fa38f19] Format v1.3.7
  [f6369f11] ForwardDiff v1.4.6
  [a85aefff] FunctionMaps v0.1.2
  [f62d2435] FunctionProperties v1.2.0
  [069b7b12] FunctionWrappers v1.1.3
  [77dc65aa] FunctionWrappersWrappers v1.13.0
  [d9f16b24] Functors v0.5.3
⌃ [46192b85] GPUArraysCore v0.2.0
⌅ [61eb1bfa] GPUCompiler v1.23.0
  [28b8d3ca] GR v0.73.27
  [a0844989] Gamma v1.2.0
  [86223c79] Graphs v1.15.0
  [076d061b] HashArrayMappedTries v0.2.0
⌅ [eafb193a] Highlights v0.5.3
  [34004b35] HypergeometricFunctions v0.3.30
  [7869d1d1] IRTools v0.4.20
  [3263718b] ImplicitDiscreteSolve v2.3.0
  [d25df0c9] Inflate v0.1.5
⌅ [842dd82b] InlineStrings v1.4.6
  [18e54dd8] IntegerMathUtils v0.1.4
  [8197267c] IntervalSets v0.7.14
  [3587e190] InverseFunctions v0.1.17
  [41ab1584] InvertedIndices v1.3.1
  [92d709cd] IrrationalConstants v0.2.6
  [82899510] IteratorInterfaceExtensions v1.0.0
  [1019f520] JLFzf v0.1.11
  [692b3bcd] JLLWrappers v1.8.0
⌅ [682c06a0] JSON v0.21.4
  [ae98c720] Jieko v0.2.1
⌃ [ccbc3e58] JumpProcesses v9.32.3
  [63c18a36] KernelAbstractions v0.9.42
  [ba0b0d4f] Krylov v0.10.10
  [2faa5264] LHLFactorization v2.2.2
  [929cbde3] LLVM v9.13.1
  [b964fa9f] LaTeXStrings v1.4.1
  [23fbe1c1] Latexify v0.16.12
  [87fe0de2] LineSearch v0.1.18
⌃ [7ed4a6bd] LinearSolve v5.17.3
  [2ab3a3ac] LogExpFunctions v1.0.1
  [e6f89c97] LoggingExtras v1.2.0
  [1914dd2f] MacroTools v0.5.16
  [bb5d69b7] MaybeInplace v0.1.8
  [442fdcdd] Measures v0.3.3
  [e1d29d7a] Missings v1.2.0
  [dbe65cb8] MistyClosures v2.1.0
⌃ [961ee093] ModelingToolkit v11.43.1
⌃ [7771a370] ModelingToolkitBase v1.71.1
  [6bb917b9] ModelingToolkitTearing v1.20.6
⌃ [da2b9cff] Mooncake v0.5.56
⌅ [2e0e35c7] Moshi v0.3.9
  [46d2c3a1] MuladdMacro v0.2.7
⌃ [102ac46a] MultivariatePolynomials v0.5.19
⌃ [ffc61752] Mustache v1.0.21 [loaded: v1.1.0]
⌃ [d8a4904e] MutableArithmetics v1.8.0
  [872c559c] NNlib v0.9.45
  [77ba4419] NaNMath v1.1.4
⌃ [8913a72c] NonlinearSolve v4.30.0
⌃ [be0214bd] NonlinearSolveBase v2.49.5
⌃ [5959db7a] NonlinearSolveFirstOrder v2.6.1
  [9a2c21bd] NonlinearSolveQuasiNewton v1.15.3
  [26075421] NonlinearSolveSpectralMethods v1.8.3
  [d8793406] ObjectFile v0.5.1
  [6fe1bfb0] OffsetArrays v1.17.0
  [3bd65402] Optimisers v0.4.9
⌅ [bac558e1] OrderedCollections v1.8.2 [loaded: v2.0.1]
  [1dea7af3] OrdinaryDiffEq v7.8.1
⌃ [6ad6398a] OrdinaryDiffEqBDF v2.4.9
⌃ [bbf590c4] OrdinaryDiffEqCore v4.17.2
  [50262376] OrdinaryDiffEqDefault v2.6.2
⌃ [4302a76b] OrdinaryDiffEqDifferentiation v3.12.0
  [127b3ac7] OrdinaryDiffEqNonlinearSolve v2.9.8
⌃ [43230ef6] OrdinaryDiffEqRosenbrock v2.7.3
  [b4bd8bb3] OrdinaryDiffEqRosenbrockTableaus v2.4.2
⌃ [2d112036] OrdinaryDiffEqSDIRK v2.9.3
  [b1df2697] OrdinaryDiffEqTsit5 v2.1.4
  [79d7bb75] OrdinaryDiffEqVerner v2.4.1
  [90014a1f] PDMats v0.11.41
  [65888b18] ParameterizedFunctions v5.27.0
⌅ [69de0a69] Parsers v2.8.8
  [ccf2f8ad] PlotThemes v3.3.0
⌃ [995b91a9] PlotUtils v1.4.4
  [91a5bcdd] Plots v1.41.7
  [e409e4f3] PoissonRandom v0.4.13
  [2dfb63ee] PooledArrays v1.4.3
  [d236fae5] PreallocationTools v1.7.1
  [aea7be01] PrecompileTools v1.3.4
  [21216c6a] Preferences v1.6.0
  [08abe8d2] PrettyTables v3.4.8
  [27ebfcd6] Primes v0.5.7
  [92933f4c] ProgressMeter v1.11.0
  [43287f4e] PtrArrays v1.4.0
⌃ [0c0d3e7f] PureKLU v1.5.0
  [1fd47b50] QuadGK v2.11.3
  [e6cf234a] RandomNumbers v1.6.0
  [988b38a3] ReadOnlyArrays v0.2.0
  [795d4caa] ReadOnlyDicts v1.0.1
  [c1ae055f] RealDot v0.1.0
  [3cdcf5f2] RecipesBase v1.3.4
  [01d81517] RecipesPipeline v0.6.12
⌃ [731186ca] RecursiveArrayTools v4.5.1
  [189a3867] Reexport v1.2.2
  [05181044] RelocatableFolders v1.0.1
  [ae029012] Requires v1.3.1
  [ae5879a3] ResettableStacks v1.4.0
  [9fe22ead] RespecializeParams v1.3.0
⌃ [37e2e3b7] ReverseDiff v1.17.0
  [79098fc4] Rmath v0.9.0
  [f2b01f46] Roots v3.0.8
  [7e49a35a] RuntimeGeneratedFunctions v0.5.26
  [9dfe8606] SCCNonlinearSolve v1.15.3
⌃ [0bca4576] SciMLBase v3.54.0
  [31c91b34] SciMLBenchmarks v0.2.1 [loaded: `/home/crackauc/github-runners/amdci3-1/_work/SciMLBenchmarks.jl/SciMLBenchmarks.jl/src/SciMLBenchmarks.jl` (v0.2.1) expected `/home/crackauc/.julia/packages/SciMLBenchmarks/ceJyd/src/SciMLBenchmarks.jl` (v0.2.1)]
  [19f34311] SciMLJacobianOperators v0.1.19
  [a6db7da4] SciMLLogging v2.1.0
⌃ [c0aeaf25] SciMLOperators v1.30.0
  [431bcebd] SciMLPublic v1.3.0
⌃ [1ed8b502] SciMLSensitivity v7.119.7
  [53ae85a6] SciMLStructures v1.10.5
  [7e506255] ScopedValues v1.6.2
  [6c6a2e73] Scratch v1.3.0
  [91c51154] SentinelArrays v1.4.10
  [efcf1570] Setfield v1.1.2
  [992d4aef] Showoff v1.1.1
  [727e6d20] SimpleNonlinearSolve v2.14.5
  [699a6c99] SimpleTraits v0.9.6
  [a2af1166] SortingAlgorithms v1.2.3
  [a57abbd0] SparseColumnPivotedQR v2.1.8
  [dc90abb0] SparseInverseSubset v0.1.3
  [0a514795] SparseMatrixColorings v0.4.28
  [276daf66] SpecialFunctions v2.9.0
  [860ef19b] StableRNGs v1.0.4
  [0c0c59c1] StarAlgebras v0.3.0
  [64909d44] StateSelection v1.11.1
⌃ [90137ffa] StaticArrays v1.9.20
  [1e83bf80] StaticArraysCore v1.4.4
  [10745b16] Statistics v1.11.5
  [82ae8749] StatsAPI v1.8.0
  [2913bbd2] StatsBase v0.34.13
  [4c63d2b9] StatsFuns v2.2.1
  [69024149] StringEncodings v0.3.7
⌅ [892a3eda] StringManipulation v0.5.0
  [09ab397b] StructArrays v0.7.3
  [53d494c1] StructIO v0.3.1
  [c3572dad] Sundials v6.7.1
  [2efcf032] SymbolicIndexingInterface v0.3.55
  [19f23fe9] SymbolicLimits v1.2.1
⌃ [d1185830] SymbolicUtils v4.46.6
⌃ [0c5d862f] Symbolics v7.39.2
  [9ce81f87] TableMetadataTools v0.1.0
  [3783bdb8] TableTraits v1.0.1
  [bd369af6] Tables v1.14.0
  [ed4db957] TaskLocalValues v0.1.3
  [62fd8b95] TensorCore v0.1.1
  [8ea1fca8] TermInterface v2.0.0
⌃ [a759f4b9] TimerOutputs v1.2.1
  [9f7883ad] Tracker v0.2.39
  [e689c965] Tracy v0.1.6
  [781d530d] TruncatedStacktraces v1.4.0
  [3a884ed6] UnPack v1.0.2
  [1cfade01] UnicodeFun v0.4.1
  [1986cc42] Unitful v1.29.0
  [013be700] UnsafeAtomics v0.3.2
  [41fe7b60] Unzip v0.2.0
  [d30d5f5c] WeakCacheSets v0.1.0
  [44d3d7a6] Weave v0.10.12
⌃ [ddb6d928] YAML v0.4.16 [loaded: v0.4.17]
  [e88e6eb3] Zygote v0.7.13
  [700de1a5] ZygoteRules v0.2.8
  [6e34b625] Bzip2_jll v1.0.9+0
  [83423d85] Cairo_jll v1.18.7+0
  [ee1fde0b] Dbus_jll v1.16.2+0
⌅ [7cc45869] Enzyme_jll v0.0.293+0
  [2702e6a9] EpollShim_jll v0.0.20230411+1
  [2e619515] Expat_jll v2.8.4+0
⌅ [b22a6f82] FFMPEG_jll v8.1.2+0
  [a3f928ae] Fontconfig_jll v2.17.1+0
  [d7e528f0] FreeType2_jll v2.14.3+1
  [559328eb] FriBidi_jll v1.0.17+0
  [0656b61e] GLFW_jll v3.5.1+0
  [d2c73de3] GR_jll v0.73.27+0
⌅ [b0724c58] GettextRuntime_jll v0.22.4+0
  [61579ee1] Ghostscript_jll v9.55.1+0
  [7746bdde] Glib_jll v2.88.3+0
  [3b182d85] Graphite2_jll v1.3.16+0
  [2e76f6c2] HarfBuzz_jll v100.14004.0+0
  [1d5cc7b8] IntelOpenMP_jll v2025.2.0+0
  [aacddb02] JpegTurbo_jll v3.2.0+1
  [c1c5ebd0] LAME_jll v3.100.3+0
  [88015f11] LERC_jll v4.2.0+0
  [dad2f222] LLVMExtra_jll v0.0.47+0
  [1d63c593] LLVMOpenMP_jll v23.1.1+0
  [ad6e5548] LibTracyClient_jll v0.13.1+0
⌅ [e9f186c6] Libffi_jll v3.4.7+0
  [7e76a0d4] Libglvnd_jll v1.7.1+1
  [94ce4f54] Libiconv_jll v1.18.0+0
  [4b2f31a3] Libmount_jll v2.42.0+0
  [89763e89] Libtiff_jll v4.7.3+0
  [38a345b3] Libuuid_jll v2.42.0+0
  [856f044c] MKL_jll v2025.2.0+0
  [e7412a2a] Ogg_jll v1.3.6+0
  [656ef2d0] OpenBLAS32_jll v0.3.34+0
  [efe28fd5] OpenSpecFun_jll v0.5.6+0
  [91d4177d] Opus_jll v1.6.1+0
  [36c8627f] Pango_jll v1.58.2+0
  [30392449] Pixman_jll v0.46.4+0
  [c0090381] Qt6Base_jll v6.10.2+2
  [629bc702] Qt6Declarative_jll v6.10.2+2
  [ce943373] Qt6ShaderTools_jll v6.10.2+1
  [6de9746b] Qt6Svg_jll v6.10.2+0
  [e99dba38] Qt6Wayland_jll v6.10.2+1
  [f50d1b31] Rmath_jll v0.5.2+0
  [ca45d3f4] SuiteSparse32_jll v7.12.1+1
  [fb77eaff] Sundials_jll v7.5.0+0
  [a44049a8] Vulkan_Loader_jll v1.3.243+0
  [a2964d1f] Wayland_jll v1.24.0+0
  [ffd25f8a] XZ_jll v5.8.4+0
  [f67eecfb] Xorg_libICE_jll v1.1.2+0
  [c834827a] Xorg_libSM_jll v1.2.6+0
  [4f6342f7] Xorg_libX11_jll v1.8.13+0
  [0c0b7dd1] Xorg_libXau_jll v1.0.13+0
  [935fb764] Xorg_libXcursor_jll v1.2.4+0
  [a3789734] Xorg_libXdmcp_jll v1.1.6+0
  [1082639a] Xorg_libXext_jll v1.3.8+0
  [d091e8ba] Xorg_libXfixes_jll v6.0.2+0
  [a51aa0fd] Xorg_libXi_jll v1.8.4+0
  [d1454406] Xorg_libXinerama_jll v1.1.7+0
  [ec84b674] Xorg_libXrandr_jll v1.5.6+0
  [ea2f1a96] Xorg_libXrender_jll v0.9.12+0
  [a65dc6b1] Xorg_libpciaccess_jll v0.19.0+0
  [c7cfdc94] Xorg_libxcb_jll v1.17.1+0
  [cc61e674] Xorg_libxkbfile_jll v1.2.0+0
  [e920d4aa] Xorg_xcb_util_cursor_jll v0.1.6+0
  [12413925] Xorg_xcb_util_image_jll v0.4.1+0
  [2def613f] Xorg_xcb_util_jll v0.4.1+0
  [975044d2] Xorg_xcb_util_keysyms_jll v0.4.1+0
  [0d47668e] Xorg_xcb_util_renderutil_jll v0.3.10+0
  [c22f9ab0] Xorg_xcb_util_wm_jll v0.4.2+0
  [35661453] Xorg_xkbcomp_jll v1.4.7+0
  [33bec58e] Xorg_xkeyboard_config_jll v2.47.0+2
  [c5fb5394] Xorg_xtrans_jll v1.6.0+0
  [3161d3a3] Zstd_jll v1.5.7+1
  [35ca27e7] eudev_jll v3.2.14+0
⌅ [214eeab7] fzf_jll v0.61.1+0
⌃ [a4ae2306] libaom_jll v3.14.1+0
  [0ac62f75] libass_jll v0.17.5+0
  [1183f4f0] libdecor_jll v0.2.2+0
  [8e53e030] libdrm_jll v2.4.134+0
  [2db6ffa8] libevdev_jll v1.13.4+0
  [f638f0a6] libfdk_aac_jll v2.0.4+0
  [36db933b] libinput_jll v1.28.1+0
  [b53b4c65] libpng_jll v1.6.58+0
  [9a156e7d] libva_jll v2.23.0+0
  [f27f6e37] libvorbis_jll v1.3.8+0
  [009596ad] mtdev_jll v1.1.7+0
  [1317d2d5] oneTBB_jll v2022.3.0+0
⌅ [1270edf5] x264_jll v10164.0.1+0
  [dfaa095f] x265_jll v4.1.0+0
  [d8fb68d0] xkbcommon_jll v1.13.0+0
  [0dad84c5] ArgTools v1.1.2
  [56f22d72] Artifacts v1.11.0
  [2a0f44e3] Base64 v1.11.0
  [ade2ca70] Dates v1.11.0
  [8ba89e20] Distributed v1.11.0
  [f43a241f] Downloads v1.7.0
  [7b1f6079] FileWatching v1.11.0
  [9fa8497b] Future v1.11.0
  [b77e0a4c] InteractiveUtils v1.11.0
  [ac6e5ff7] JuliaSyntaxHighlighting v1.12.0
  [4af54fe1] LazyArtifacts v1.11.0
  [b27032c2] LibCURL v0.6.4
  [76f85450] LibGit2 v1.11.0
  [8f399da3] Libdl v1.11.0
  [37e2e46d] LinearAlgebra v1.12.0
  [56ddb016] Logging v1.11.0
  [d6f4376e] Markdown v1.11.0
  [a63ad114] Mmap v1.11.0
  [ca575930] NetworkOptions v1.3.0
  [44cfe95a] Pkg v1.12.1
  [de0858da] Printf v1.11.0
  [9abbd945] Profile v1.11.0
  [3fa0cd96] REPL v1.11.0
  [9a3f8284] Random v1.11.0
  [ea8e919c] SHA v0.7.0
  [9e88b42a] Serialization v1.11.0
  [6462fe0b] Sockets v1.11.0
  [2f01184e] SparseArrays v1.12.0
  [f489334b] StyledStrings v1.11.0
  [4607b0f0] SuiteSparse
  [fa267f1f] TOML v1.0.3
  [a4e569a6] Tar v1.10.0
  [8dfed614] Test v1.11.0
  [cf7118a7] UUIDs v1.11.0
  [4ec0a83e] Unicode v1.11.0
  [e66e0078] CompilerSupportLibraries_jll v1.3.1+2
  [deac9b47] LibCURL_jll v8.15.0+0
  [e37daf67] LibGit2_jll v1.9.0+0
  [29816b5a] LibSSH2_jll v1.11.3+1
  [14a3606d] MozillaCACerts_jll v2025.11.4
  [4536629a] OpenBLAS_jll v0.3.29+0
  [05823500] OpenLibm_jll v0.8.7+0
  [458c3c95] OpenSSL_jll v3.5.6+0
  [efcefdf7] PCRE2_jll v10.44.0+1
  [bea87d4a] SuiteSparse_jll v7.8.3+2
  [83775a58] Zlib_jll v1.3.1+2
  [8e850b90] libblastrampoline_jll v5.15.0+0
  [8e850ede] nghttp2_jll v1.64.0+1
  [3f19e933] p7zip_jll v17.7.0+0
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m`