From 73d62782087e77347df50586348b146bc1871d78 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Fri, 22 May 2026 19:14:29 +0200 Subject: [PATCH 1/4] Update precomputed NHS with padding to allow for less frequent updates --- src/nhs_precomputed.jl | 65 ++++++++++++++++++++++++++++-------------- 1 file changed, 44 insertions(+), 21 deletions(-) diff --git a/src/nhs_precomputed.jl b/src/nhs_precomputed.jl index c7ddceca..07e9eba2 100644 --- a/src/nhs_precomputed.jl +++ b/src/nhs_precomputed.jl @@ -2,6 +2,7 @@ PrecomputedNeighborhoodSearch{NDIMS}(; search_radius = 0.0, n_points = 0, periodic_box = nothing, update_strategy = nothing, update_neighborhood_search = GridNeighborhoodSearch{NDIMS}(), + update_neighborhood_search_padding = 0, backend = DynamicVectorOfVectors{Int32}, transpose_backend = false, max_neighbors = max_neighbors(NDIMS)) @@ -40,6 +41,9 @@ to strip the internal neighborhood search, which is not needed anymore. If the precomputed NHS is to be used on the GPU, make sure to either freeze it after initialization and never update it again, or pass a GPU-compatible neighborhood search here. +- `update_neighborhood_search_padding = 0`: Relative padding used for the fixed + search radius of the internal [`GridNeighborhoodSearch`](@ref) + that computes the neighbor lists. - `backend = DynamicVectorOfVectors{Int32}`: Type of the data structure to store the neighbor lists. Can be - `Vector{Vector{Int32}}`: Scattered memory, but very memory-efficient. @@ -66,31 +70,37 @@ to strip the internal neighborhood search, which is not needed anymore. """ struct PrecomputedNeighborhoodSearch{NDIMS, NL, ELTYPE, PB, NHS} <: AbstractNeighborhoodSearch - neighbor_lists :: NL - search_radius :: ELTYPE - periodic_box :: PB - neighborhood_search :: NHS - sort_neighbor_lists :: Bool + neighbor_lists :: NL + search_radius :: ELTYPE + periodic_box :: PB + neighborhood_search :: NHS + sort_neighbor_lists :: Bool + update_neighborhood_search_padding :: Float64 function PrecomputedNeighborhoodSearch{NDIMS}(neighbor_lists, search_radius, periodic_box, update_neighborhood_search, - sort_neighbor_lists) where {NDIMS} + sort_neighbor_lists, + update_neighborhood_search_padding) where {NDIMS} return new{NDIMS, typeof(neighbor_lists), typeof(search_radius), typeof(periodic_box), typeof(update_neighborhood_search)}(neighbor_lists, search_radius, periodic_box, update_neighborhood_search, - sort_neighbor_lists) + sort_neighbor_lists, + update_neighborhood_search_padding) end end function PrecomputedNeighborhoodSearch{NDIMS}(; search_radius = 0.0, n_points = 0, periodic_box = nothing, update_strategy = nothing, + update_neighborhood_search_padding = 0, update_neighborhood_search = GridNeighborhoodSearch{NDIMS}(; - search_radius, + search_radius = search_radius * + (1 + + update_neighborhood_search_padding), n_points, periodic_box, update_strategy), @@ -102,7 +112,8 @@ function PrecomputedNeighborhoodSearch{NDIMS}(; search_radius = 0.0, n_points = PrecomputedNeighborhoodSearch{NDIMS}(neighbor_lists, search_radius, periodic_box, update_neighborhood_search, - sort_neighbor_lists) + sort_neighbor_lists, + update_neighborhood_search_padding) end # Default values for maximum neighbor count @@ -136,6 +147,7 @@ function initialize!(search::PrecomputedNeighborhoodSearch, initialize!(neighborhood_search, x, y; parallelization_backend) initialize_neighbor_lists!(neighbor_lists, neighborhood_search, x, y, + search.search_radius, parallelization_backend, search.sort_neighbor_lists) return search @@ -144,7 +156,7 @@ end function update!(search::PrecomputedNeighborhoodSearch, x::AbstractMatrix, y::AbstractMatrix; points_moving = (true, true), parallelization_backend = default_backend(x), - eachindex_y = axes(y, 2)) + eachindex_y = axes(y, 2), search_radius = search.search_radius) (; neighborhood_search, neighbor_lists) = search if eachindex_y != axes(y, 2) @@ -157,6 +169,7 @@ function update!(search::PrecomputedNeighborhoodSearch, # Skip update if both point sets are static if any(points_moving) initialize_neighbor_lists!(neighbor_lists, neighborhood_search, x, y, + search_radius, parallelization_backend, search.sort_neighbor_lists) end @@ -164,7 +177,8 @@ function update!(search::PrecomputedNeighborhoodSearch, end function initialize_neighbor_lists!(neighbor_lists, neighborhood_search, x, y, - parallelization_backend, sort_neighbor_lists) + search_radius, parallelization_backend, + sort_neighbor_lists) # Initialize neighbor lists empty!(neighbor_lists) resize!(neighbor_lists, size(x, 2)) @@ -174,14 +188,16 @@ function initialize_neighbor_lists!(neighbor_lists, neighborhood_search, x, y, # Fill neighbor lists foreach_point_neighbor(x, y, neighborhood_search; - parallelization_backend) do point, neighbor, _, _ - push!(neighbor_lists[point], neighbor) + parallelization_backend) do point, neighbor, _, distance + if distance <= search_radius + push!(neighbor_lists[point], neighbor) + end end end function initialize_neighbor_lists!(neighbor_lists::DynamicVectorOfVectors, - neighborhood_search, x, y, parallelization_backend, - sort_neighbor_lists) + neighborhood_search, x, y, search_radius, + parallelization_backend, sort_neighbor_lists) resize!(neighbor_lists, size(x, 2)) # `Base.empty!.(neighbor_lists)`, but for all backends @@ -191,8 +207,10 @@ function initialize_neighbor_lists!(neighbor_lists::DynamicVectorOfVectors, # Fill neighbor lists foreach_point_neighbor(x, y, neighborhood_search; - parallelization_backend) do point, neighbor, _, _ - pushat!(neighbor_lists, point, neighbor) + parallelization_backend) do point, neighbor, _, distance + if distance <= search_radius + pushat!(neighbor_lists, point, neighbor) + end end if sort_neighbor_lists @@ -238,19 +256,23 @@ end function copy_neighborhood_search(nhs::PrecomputedNeighborhoodSearch, search_radius, n_points; eachpoint = 1:n_points) update_neighborhood_search = copy_neighborhood_search(nhs.neighborhood_search, - search_radius, n_points; - eachpoint) + search_radius * (1 + + nhs.update_neighborhood_search_padding), + n_points; eachpoint) # For `Vector{Vector}` backend use `max_neighbors(NDIMS)` as fallback. # This should never be used because this backend doesn't require a `max_neighbors`. max_neighbors_ = max_inner_length(nhs.neighbor_lists, max_neighbors(ndims(nhs))) transpose_backend = transposed_backend(nhs.neighbor_lists) + update_neighborhood_search_padding = nhs.update_neighborhood_search_padding return PrecomputedNeighborhoodSearch{ndims(nhs)}(; search_radius, n_points, periodic_box = nhs.periodic_box, update_neighborhood_search, backend = typeof(nhs.neighbor_lists), transpose_backend, - max_neighbors = max_neighbors_) + max_neighbors = max_neighbors_, + sort_neighbor_lists = nhs.sort_neighbor_lists, + update_neighborhood_search_padding) end @inline function freeze_neighborhood_search(search::PrecomputedNeighborhoodSearch) @@ -261,5 +283,6 @@ end search.search_radius, search.periodic_box, nothing, - search.sort_neighbor_lists) + search.sort_neighbor_lists, + search.update_neighborhood_search_padding) end From 4c638d6b55c91180359afbccf9125b4bbc15bd71 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Mon, 25 May 2026 17:32:20 +0200 Subject: [PATCH 2/4] Reforat --- src/nhs_precomputed.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/nhs_precomputed.jl b/src/nhs_precomputed.jl index 07e9eba2..cce1a70a 100644 --- a/src/nhs_precomputed.jl +++ b/src/nhs_precomputed.jl @@ -168,8 +168,7 @@ function update!(search::PrecomputedNeighborhoodSearch, # Skip update if both point sets are static if any(points_moving) - initialize_neighbor_lists!(neighbor_lists, neighborhood_search, x, y, - search_radius, + initialize_neighbor_lists!(neighbor_lists, neighborhood_search, x, y, search_radius, parallelization_backend, search.sort_neighbor_lists) end From 9925bbea3ebe66aa4686d2827f26ab4b793b6a4c Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Mon, 25 May 2026 17:41:33 +0200 Subject: [PATCH 3/4] Add additional docs and tests --- src/nhs_precomputed.jl | 32 ++++++++++++++++++++++++++++---- test/nhs_precomputed.jl | 21 +++++++++++++++++++++ 2 files changed, 49 insertions(+), 4 deletions(-) diff --git a/src/nhs_precomputed.jl b/src/nhs_precomputed.jl index cce1a70a..ac54ec90 100644 --- a/src/nhs_precomputed.jl +++ b/src/nhs_precomputed.jl @@ -42,8 +42,16 @@ to strip the internal neighborhood search, which is not needed anymore. either freeze it after initialization and never update it again, or pass a GPU-compatible neighborhood search here. - `update_neighborhood_search_padding = 0`: Relative padding used for the fixed - search radius of the internal [`GridNeighborhoodSearch`](@ref) - that computes the neighbor lists. + search radius of the default internal + [`GridNeighborhoodSearch`](@ref) that computes the neighbor + lists. For example, `update_neighborhood_search_padding = 0.05` + creates the internal search with a 5% larger radius than + `search_radius`. This allows [`update!`](@ref) to rebuild the + precomputed neighbor lists with a larger `search_radius`, so the + lists can be updated less frequently when points move only + slowly. The search radius of the internal update neighborhood + search must be at least as large as the `search_radius` passed + to `update!` to build the precomputed neighbor lists. - `backend = DynamicVectorOfVectors{Int32}`: Type of the data structure to store the neighbor lists. Can be - `Vector{Vector{Int32}}`: Scattered memory, but very memory-efficient. @@ -146,9 +154,11 @@ function initialize!(search::PrecomputedNeighborhoodSearch, # Initialize grid NHS initialize!(neighborhood_search, x, y; parallelization_backend) + check_update_neighborhood_search_radius(search, search.search_radius) + initialize_neighbor_lists!(neighbor_lists, neighborhood_search, x, y, - search.search_radius, - parallelization_backend, search.sort_neighbor_lists) + search.search_radius, parallelization_backend, + search.sort_neighbor_lists) return search end @@ -168,6 +178,8 @@ function update!(search::PrecomputedNeighborhoodSearch, # Skip update if both point sets are static if any(points_moving) + check_update_neighborhood_search_radius(search, search_radius) + initialize_neighbor_lists!(neighbor_lists, neighborhood_search, x, y, search_radius, parallelization_backend, search.sort_neighbor_lists) end @@ -175,6 +187,18 @@ function update!(search::PrecomputedNeighborhoodSearch, return search end +function check_update_neighborhood_search_radius(search::PrecomputedNeighborhoodSearch, + neighbor_list_search_radius) + update_search_radius = search_radius(search.neighborhood_search) + + if update_search_radius < neighbor_list_search_radius + throw(ArgumentError("the search radius of the internal update neighborhood search " * + "($update_search_radius) must be at least as large as the " * + "search radius used to build the precomputed neighbor lists " * + "($neighbor_list_search_radius)")) + end +end + function initialize_neighbor_lists!(neighbor_lists, neighborhood_search, x, y, search_radius, parallelization_backend, sort_neighbor_lists) diff --git a/test/nhs_precomputed.jl b/test/nhs_precomputed.jl index fe084352..872a24d0 100644 --- a/test/nhs_precomputed.jl +++ b/test/nhs_precomputed.jl @@ -35,4 +35,25 @@ pointer_ = pointer(neighbor_lists.backend.parent) @test unsafe_load(pointer_, 1) == 101 @test unsafe_load(pointer_, 2) == 201 + + @testset "`update_neighborhood_search_padding`" begin + search_radius = 1.0 + update_neighborhood_search_padding = 0.05 + nhs = PrecomputedNeighborhoodSearch{2}(; search_radius, + update_neighborhood_search_padding, + n_points = 2) + + @test PointNeighbors.search_radius(nhs.neighborhood_search) == 1.05 + + coordinates = [0.0 1.0 + 0.0 0.0] + + initialize!(nhs, coordinates, coordinates) + @test update!(nhs, coordinates, coordinates; search_radius = 1.05) === nhs + + error = ArgumentError("the search radius of the internal update neighborhood " * + "search (1.05) must be at least as large as the search " * + "radius used to build the precomputed neighbor lists (1.06)") + @test_throws error update!(nhs, coordinates, coordinates; search_radius = 1.06) + end end From b4dcfad151b0923a8f1c7e86be31a031e7eaf99b Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Mon, 25 May 2026 19:27:03 +0200 Subject: [PATCH 4/4] Fix --- src/gpu.jl | 3 ++- test/nhs_precomputed.jl | 18 ++++++++++++++++++ 2 files changed, 20 insertions(+), 1 deletion(-) diff --git a/src/gpu.jl b/src/gpu.jl index 1de3dbb4..31728da3 100644 --- a/src/gpu.jl +++ b/src/gpu.jl @@ -31,7 +31,8 @@ function Adapt.adapt_structure(to, nhs::PrecomputedNeighborhoodSearch) return PrecomputedNeighborhoodSearch{ndims(nhs)}(neighbor_lists, search_radius, periodic_box, neighborhood_search, - nhs.sort_neighbor_lists) + nhs.sort_neighbor_lists, + nhs.update_neighborhood_search_padding) end function Adapt.adapt_structure(to, cell_list::SpatialHashingCellList{NDIMS}) where {NDIMS} diff --git a/test/nhs_precomputed.jl b/test/nhs_precomputed.jl index 872a24d0..e4c3801f 100644 --- a/test/nhs_precomputed.jl +++ b/test/nhs_precomputed.jl @@ -56,4 +56,22 @@ "radius used to build the precomputed neighbor lists (1.06)") @test_throws error update!(nhs, coordinates, coordinates; search_radius = 1.06) end + + @testset "`Adapt.adapt_structure` preserves all fields" begin + nhs = PrecomputedNeighborhoodSearch{3}(; search_radius = 1.0f0, + n_points = 2, + update_neighborhood_search_padding = 0.05) + frozen_nhs = PointNeighbors.freeze_neighborhood_search(nhs) + + adapted_nhs = @trixi_test_nowarn PointNeighbors.Adapt.adapt_structure(Array, + frozen_nhs) + + @test adapted_nhs.neighbor_lists == frozen_nhs.neighbor_lists + @test adapted_nhs.search_radius == frozen_nhs.search_radius + @test adapted_nhs.periodic_box == frozen_nhs.periodic_box + @test adapted_nhs.neighborhood_search == frozen_nhs.neighborhood_search + @test adapted_nhs.sort_neighbor_lists == frozen_nhs.sort_neighbor_lists + @test adapted_nhs.update_neighborhood_search_padding == + frozen_nhs.update_neighborhood_search_padding + end end