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/src/nhs_precomputed.jl b/src/nhs_precomputed.jl index c7ddceca..ac54ec90 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,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. @@ -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), @@ -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 @@ -135,8 +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, - parallelization_backend, search.sort_neighbor_lists) + search.search_radius, parallelization_backend, + search.sort_neighbor_lists) return search end @@ -144,7 +166,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) @@ -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)) @@ -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 @@ -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 @@ -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) @@ -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 diff --git a/test/nhs_precomputed.jl b/test/nhs_precomputed.jl index fe084352..e4c3801f 100644 --- a/test/nhs_precomputed.jl +++ b/test/nhs_precomputed.jl @@ -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