Memory leak in threaded calls to CVODE

[termi-official]

I think calling CVODE from threads causes memory leaks. A somewhat minimal reproducer is given by the following code sample:

using ModelingToolkit, OrdinaryDiffEq, Sundials, Evolutionary 

function objective(p, solver=CVODE_BDF())
    @parameters t a b c d e f g h i
    @variables O(t), P1(t), P2(t), P3(t)

    D = Differential(t)
    
    t_offset = 60000.0
    eqs = [
        D(O) ~ d*P2 + h*P3 + a*P1*exp(i*ifelse(t_offset <= t, ifelse(t <= t_offset+430.0, 70.0, 0.0), 0.0)) - b*O - c*O - g*O,
        D(P1) ~ ((-a*c*e) / (b*d))*P1 + b*O + e*P2*exp(-i*ifelse(t_offset <= t, ifelse(t <= t_offset+430.0, 70.0, 0.0), 0.0)) - a*P1*exp(i*ifelse(t_offset <= t, ifelse(t <= t_offset+430.0, 70.0, 0.0), 0.0)),
        D(P2) ~ c*O - ((-a*c*e) / (b*d))*P1 - d*P2 - e*P2*exp(-i*ifelse(t_offset <= t, ifelse(t <= t_offset+430.0, 70.0, 0.0), 0.0)),
        D(P3) ~ g*O - h*P3
    ]

    @named sys = ODESystem(eqs)
    sys = structural_simplify(sys)
    
    tspan = (0.0, t_offset+1100.0)
    u₀ = [
        O => 0.0,
        P1 => 1.0,
        P2 => 0.0,
        P3 => 0.0
    ]
    prob = ODEProblem(sys, u₀, tspan, p)
    solution = solve(prob, solver; dtmax=10.0)
    if solution.retcode != ReturnCode.Success
        return Inf
    end

    t = [17.26973684211, 31.08552631579, 54.11184210526, 81.74342105263, 117.43421052631, 207.23684210526, 316.61184210526, 427.13815789474, 572.20394736842, 1046.54605263158] .+ t_offset
    soleval = solution(t)
    x = [0.8545454545454545, 1.0, 0.7818181818181817, 0.5909090909090908, 0.3909090909090909, 0.23636363636363633, 0.09999999999999999, 0.09090909090909091, 0.06363636363636363, 0.01818181818181818]

    sum2 = sum((soleval[1, :] - x).^2)
    return sqrt(sum2/length(x))
end

Tmax = 5.0
μ = 10
initial_population = [[rand()+1e-16 for _ ∈ 1:9] for _ ∈ 1:μ]
plower = [1e-16 for i ∈ 1:9]
pupper = [1e0 for i ∈ 1:9]
strat = CMAES(;μ=μ, sigma0=0.1, metrics=[Evolutionary.AbsDiff(1e-11)], c_1=0.7, c_mu=0.3, c_c=0.95, c_sigma=0.95)
p_opt = Evolutionary.optimize(p -> objective(p), BoxConstraints(plower, pupper), strat, initial_population, Evolutionary.Options(show_trace=true, parallelization=:thread, time_limit=Tmax))
GC.gc() # Manually trigger garbage collection
# 7.5 GB usage

# Threaded Sundials
Tmax = 60.0
p_opt = Evolutionary.optimize(p -> objective(p), BoxConstraints(plower, pupper), strat, initial_population, Evolutionary.Options(show_trace=true, parallelization=:thread, time_limit=Tmax))
GC.gc() # Manually trigger garbage collection
# 13.1 GB usage

Tmax = 60.0
p_opt = Evolutionary.optimize(p -> objective(p), BoxConstraints(plower, pupper), strat, initial_population, Evolutionary.Options(show_trace=true, time_limit=Tmax))
GC.gc() # Manually trigger garbage collection
# 13.1 GB usage

Tmax = 60.0
p_opt = Evolutionary.optimize(p -> objective(p, SDIRK2()), BoxConstraints(plower, pupper), strat, initial_population, Evolutionary.Options(show_trace=true, parallelization=:thread, time_limit=Tmax))
GC.gc() # Manually trigger garbage collection
# 13.1 GB usage

However, I am not super sure if the SUNDIALS interface is thread-safe. From the SUNDIALS docs I think it should be (see https://sundials.readthedocs.io/en/latest/sundials/SUNContext_link.html#implications-for-task-based-programming-and-multi-threading).