# Memory allocation when solving ODE with Vern7 and ContinuousCallback

Hello,

I went through a memory allocation problem in a ODE system, when using Vern7() algorithm and ContinuousCallback. The minimal working example is the following:

``````using LinearAlgebra
using SparseArrays
using OrdinaryDiffEq
using DiffEqCallbacks
using BenchmarkTools

A = Diagonal(Array{ComplexF64}(1:20))
A += spdiagm(1=>sqrt.(0.01 .* Array(1:19)))
A -= 0.01im * Diagonal(Array{ComplexF64}(1:20))
u0 = rand(ComplexF64, 20)
prm = (U = A, rand_n = Ref(rand()))

my_condition_discrete(u, t, integrator) = real(dot(u, u)) < integrator.p.rand_n[]
my_condition_continuous(u, t, integrator) = integrator.p.random_n[] - real(dot(u, u))
my_affect!(integrator) = u_modified!(integrator, false)
cb1 = DiscreteCallback(my_condition_discrete, my_affect!, save_positions=(false, false))
cb2 = ContinuousCallback(my_condition_continuous, my_affect!, nothing, save_positions=(false, false))

prob1 = ODEProblem((du,u,p,t)->mul!(du, p.U, u), u0, (0.0, 100.0), prm, callback=cb1)
prob2 = ODEProblem((du,u,p,t)->mul!(du, p.U, u), u0, (0.0, 100.0), prm, callback=cb2)

@benchmark sol1 = solve(\$prob1, \$Vern7(), saveat=\$[100.0])

BenchmarkTools.Trial: 9021 samples with 1 evaluation.
Range (min β¦ max):  383.100 ΞΌs β¦   5.543 ms  β GC (min β¦ max): 0.00% β¦ 0.00%
Time  (median):     475.300 ΞΌs               β GC (median):    0.00%
Time  (mean Β± Ο):   546.050 ΞΌs Β± 198.079 ΞΌs  β GC (mean Β± Ο):  0.07% Β± 0.91%

ββ
βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ β
383 ΞΌs           Histogram: frequency by time          1.2 ms <

Memory estimate: 10.34 KiB, allocs estimate: 50.

@benchmark sol2 = solve(\$prob2, \$Vern7(), saveat=\$[100.0])

BenchmarkTools.Trial: 2566 samples with 1 evaluation.
Range (min β¦ max):  1.274 ms β¦  12.947 ms  β GC (min β¦ max): 0.00% β¦ 79.46%
Time  (median):     1.709 ms               β GC (median):    0.00%
Time  (mean Β± Ο):   1.931 ms Β± 783.272 ΞΌs  β GC (mean Β± Ο):  2.63% Β±  7.29%

βββββββββββ
βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ β
1.27 ms         Histogram: frequency by time        4.71 ms <

Memory estimate: 784.83 KiB, allocs estimate: 2032.
``````

While it doesnβt happen for example with the `Tsit5()` algorithm:

``````@benchmark sol1 = solve(\$prob1, \$Tsit5(), saveat=\$[100.0])

BenchmarkTools.Trial: 10000 samples with 1 evaluation.
Range (min β¦ max):  309.600 ΞΌs β¦   5.215 ms  β GC (min β¦ max): 0.00% β¦ 89.28%
Time  (median):     392.600 ΞΌs               β GC (median):    0.00%
Time  (mean Β± Ο):   436.892 ΞΌs Β± 136.050 ΞΌs  β GC (mean Β± Ο):  0.11% Β±  0.89%

ββ
βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ β
310 ΞΌs           Histogram: frequency by time          859 ΞΌs <

Memory estimate: 9.44 KiB, allocs estimate: 48.

@benchmark sol2 = solve(\$prob2, \$Tsit5(), saveat=\$[100.0])

BenchmarkTools.Trial: 4593 samples with 1 evaluation.
Range (min β¦ max):  744.500 ΞΌs β¦   3.187 ms  β GC (min β¦ max): 0.00% β¦ 0.00%
Time  (median):     968.600 ΞΌs               β GC (median):    0.00%
Time  (mean Β± Ο):     1.079 ms Β± 306.638 ΞΌs  β GC (mean Β± Ο):  0.00% Β± 0.00%

βββ
βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ β
744 ΞΌs           Histogram: frequency by time         2.07 ms <

Memory estimate: 9.69 KiB, allocs estimate: 48.
``````
1 Like

Yes itβs something odd with type inference. This is a good example of it and Iβll see if someone from the compiler team can help figure out what it is. Iβve seen this for years but havenβt figured out how to handle it yet, but this is a cleaner example than the others I had.

2 Likes
``````using LinearAlgebra
using SparseArrays
using OrdinaryDiffEq
using DiffEqCallbacks
using BenchmarkTools

A = Diagonal(Array{ComplexF64}(1:20))
A += spdiagm(1=>sqrt.(0.01 .* Array(1:19)))
A -= 0.01im * Diagonal(Array{ComplexF64}(1:20))
u0 = rand(ComplexF64, 20)
prm = (U = A, rand_n = Ref(rand()))

my_condition_discrete(u, t, integrator) = real(dot(u, u)) < integrator.p.rand_n[]
my_condition_continuous(u, t, integrator) = integrator.p.rand_n[] - real(dot(u, u))
my_affect!(integrator) = u_modified!(integrator, false)
cb1 = DiscreteCallback(my_condition_discrete, my_affect!, save_positions=(false, false))
cb2 = ContinuousCallback(my_condition_continuous, my_affect!, nothing, save_positions=(false, false))

prob1 = ODEProblem((du,u,p,t)->mul!(du, p.U, u), u0, (0.0, 100.0), prm, callback=cb1)
prob2 = ODEProblem((du,u,p,t)->mul!(du, p.U, u), u0, (0.0, 100.0), prm, callback=cb2)

@benchmark sol1 = solve(\$prob1, \$Vern7(lazy=false), saveat=\$[100.0])

@benchmark sol2 = solve(\$prob2, \$Vern7(lazy=false), saveat=\$[100.0])
``````

The actual reason for this is because `Vern7` is uses a lazy interpolation which only extends to every `f` call if the interpolation is used, but if you use a ContinuousCallback it will always extend, in which case simply turning off the lazy interpolation would be more efficient.