Allocations in Navier-Stokes boundary condition implementations

[jlchan]

The implementation of Navier-Stokes wall boundary conditions causes additional allocations. For example, in elixir_navierstokes_convergence.jl, calc_single_boundary_flux results in a significant number of allocations.

This appears to be due to allocations related to the velocity boundary function:

velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2:3])

Calling velocity_bc_top_bottom.boundary_value_function(x, t, equations) results in allocations due to https://discourse.julialang.org/t/memory-allocation-when-evaluating-a-function-inside-struct/78784/2

[jlchan]

This isn't unique to parabolic BCs either - it appears in all of the boundary condition implementations which require wrapping a function. For example:

equations = CompressibleEulerEquations2D(1.4)
initial_condition = initial_condition_convergence_test
boundary_condition_convergence_test = BoundaryConditionDirichlet(initial_condition)

yields

julia> @btime boundary_condition_convergence_test.boundary_value_function($(SVector(1.,1.)), $.1, $equations)
  64.840 ns (4 allocations: 128 bytes)

[sloede]

I think I ran into similar issues while benchmarking the pointer wrappers for P4est.jl. If you run these benchmarks on the REPL, you get different results than when you run them inside a function. I don't know exactly when and where it happens, but if you wrap your function call like so

julia> function foo(bc, equations)
         @btime $bc.boundary_value_function($(SVector(1.,1.)), $.1, $equations)
       end
foo (generic function with 2 methods)

julia> foo(boundary_condition_convergence_test, equations)
  6.850 ns (0 allocations: 0 bytes)
4-element SVector{4, Float64} with indices SOneTo(4):
 1.9690983005625053
 1.9690983005625053
 1.9690983005625053
 3.8773481172781468

or like so

julia> function foo2(bc, equations)
         @time bc.boundary_value_function((SVector(1.,1.)), .1, equations)
       end
foo2 (generic function with 1 method)

julia> foo2(boundary_condition_convergence_test, equations)
  0.000001 seconds
4-element SVector{4, Float64} with indices SOneTo(4):
 1.9690983005625053
 1.9690983005625053
 1.9690983005625053
 3.8773481172781468

you can see that actually no allocations occur.

[jlchan]

Hm. Julia Slack suggested interpolating $(bc.boundary_value_function(...)), and I saw no allocations after that. However, there are still allocations reported by @trixi_timeit (more for DGMultiMesh, but they seem to be present for TreeMesh as well)

[sloede]

And you are sure they are caused by the BC function unwrapping and not something else?

[jlchan]

After looking with Cthulu.jl, I don't think they are anymore. I think I found a few type instabilities that may be the culprit, so I'll close this issue.