Accessing the solution of differential equations: understanding alloc

General Usage
1

I’m preallocating an array for the solution of an ODE at each time step, for post-processing. However, julia --track-allocation=user foo.jl still reports large allocations. To simplify with @btime:

using DifferentialEquations
using BenchmarkTools

function f(du,u,p,t)
    for i = 1:10
        du[i] = -u[i]*sin(i) + 1*im
    end
end

u0 = [cos(i)+im*sin(i) for i = 1:10]
tf = 10.0
tspan = (0.0, tf)
prob = ODEProblem(f, u0, tspan)

sol = solve(
    prob,
    Vern8(),
    saveat = 0.0:.01:tf,
)

sol_tmp = similar(u0)
ns = Int(tf/.01)

@btime for sol_i in sol.u
    copyto!(sol_tmp, sol_i)
end

@btime for i in 1:length(sol.u)
    copyto!(sol_tmp, sol.u[i])
end

@btime copyto!(sol_tmp, sol.u[1])

returns

  66.102 μs (1492 allocations: 38.95 KiB)
  110.124 μs (1983 allocations: 46.64 KiB)
  65.983 ns (0 allocations: 0 bytes)

note that accessing via iterator (1st approach) is ~1.6x faster here than accessing by index (2nd approach); a difficult thing for me to understand is why @btime copyto!(sol_tmp, sol.u[1]) does not allocate while it allocates in loop. Note that in loops, normal arrays do not allocate with copyto!:

a = [rand(3,3) for i = 1:1000];
b = rand(3,3)
@btime for i in 1:1000
          copyto!($b,$a[i])
       end

gives

4.603 μs (0 allocations: 0 bytes)
2

your bechmarking is using global variables, allocations will be caused by this.
Also, if you’re using @btime, you do not need to write a loop yourself

3

Would it be the same with @inbounds on index iteration?

4

thanks! i’ll look into global var. effects, but the for loop is not to do repeated measurements but to iterate through the sol.u array, in real program there’s post processing after that, the minimum example’s just to illustrate the allocations.