[BoundaryValueDiffEq.jl] Reducing allocations

Hi,

I started to implement my Python pipeline for solving boundary value problems in Julia using BoundaryValueDiffEq.jl. I tried to optimize the functions but I am wondering if it is possible to increase performance further by reducing the number of allocations. In the minimal reproducible example below I can’t get below 440 k allocations. Is there any additional optimization I could try or are the numbers of allocations I am getting as good as it gets?

import BoundaryValueDiffEq as BVP


function distribution_fun_test!(dydt::AbstractVector{T}, y::AbstractVector{T},  p::Any, #
                                t::Float64) where {T<:Real}

    y1, y2, y3, y4 = y

    dydt[1] = 0.001 * y2
    dydt[2] = 10 * y1
    dydt[3] = -0.03 * y2 + 0.0001 * y4
    dydt[4] = 10 * y3
                    
    uptake = 0.1 * y3 / (1000 * y3 + 1)
    dydt[3] += uptake
end


function distribution_bc_test!(residual::AbstractVector{}, u::Any, p::Any, t::Vector{Float64})

    residual[1] = (u[4, end] - 0.001)
    residual[2] = u[1, end]
    residual[3] =   (- u[3, end])
    residual[4] =   (u[1, 1])

end

function test_solve()
    L = 100.0
    mesh_size = 500
    tspan::Tuple{Float64, Float64} = (0.0, Float64(L))

    y_initial::Vector{Float64} = [1, 0.001/0.03, 1, 0.001]

    bvp_problem = BVP.BVProblem(distribution_fun_test!, 
                                distribution_bc_test!, 
                                y_initial, tspan)
        
    @time BVP.solve(bvp_problem, BVP.MIRK5(); adaptive = false, dt=Float64(L/mesh_size))
    @time BVP.solve(bvp_problem, BVP.MIRK5(); adaptive = false, dt=Float64(L/mesh_size))
end

test_solve()

The output of the timing I am getting is:
3.606669 seconds (4.35 M allocations: 216.794 MiB, 1.00% gc time, 96.31% compilation time: 100% of which was recompilation) 0.132779 seconds (440.05 k allocations: 26.994 MiB)

I am using BoundaryValueDiffEq v5.18.0 and Julia v1.12.

Thanks a lot for your help in advance!
Max

Welcome to the Julia community.

It looks like… um… I think the performance bottleneck is in the BoundaryValueDiffEq package (rather than the code you wrote). The type instability is also caused by this package. :thinking:

Thanks! And thanks for your reply!

It’s a pitty that there seems to be nothing I can change in my functions.

Maybe there is a way to optimize how I set up BVProblem instead?
I tried setting the the parameter nlls to false or true for type stability as suggested in the docs - but this seemed to increase allocations and increased execution time.
I also tried returning a staticarray for both functions instead of in place changes (I supplied y_array as a staticarray, too). But this increased allocations massively.

As these now are quite specific points for the boundaryvaluediffeq package - should I post the question somewhere else?

You can profile exactly where the allocations are coming from by using @profview_allocs (if you’re using VSCode)