Skip to content

Update precomputed NHS with padding to allow for less frequent updates #159

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
efaulhaber wants to merge 4 commits into
base: main
Choose a base branch
from
Draft

Update precomputed NHS with padding to allow for less frequent updates #159

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
73d6278
Update precomputed NHS with padding to allow for less frequent updates
efaulhaber May 22, 2026
4c638d6
Reforat
efaulhaber May 25, 2026
9925bbe
Add additional docs and tests
efaulhaber May 25, 2026
b4dcfad
Fix
efaulhaber May 25, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion src/gpu.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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}
Expand Down
92 changes: 69 additions & 23 deletions src/nhs_precomputed.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down Expand Up @@ -40,6 +41,17 @@ 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 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.
Expand All @@ -66,31 +78,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),
Expand All @@ -102,7 +120,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
Expand Down Expand Up @@ -135,16 +154,19 @@ 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,
parallelization_backend, search.sort_neighbor_lists)
search.search_radius, parallelization_backend,
search.sort_neighbor_lists)

return search
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)
Expand All @@ -156,15 +178,30 @@ 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,
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

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,
parallelization_backend, sort_neighbor_lists)
search_radius, parallelization_backend,
sort_neighbor_lists)
# Initialize neighbor lists
empty!(neighbor_lists)
resize!(neighbor_lists, size(x, 2))
Expand All @@ -174,14 +211,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
Expand All @@ -191,8 +230,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
Expand Down Expand Up @@ -238,19 +279,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)
Expand All @@ -261,5 +306,6 @@ end
search.search_radius,
search.periodic_box,
nothing,
search.sort_neighbor_lists)
search.sort_neighbor_lists,
search.update_neighborhood_search_padding)
end
39 changes: 39 additions & 0 deletions test/nhs_precomputed.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,4 +35,43 @@
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

@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
Loading