From e6aef57176d0880c3b4a24d1f8b65a563cc19a35 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Thu, 9 Nov 2023 22:32:54 +0100 Subject: [PATCH 1/4] [skip ci] Add `NeighborListNeighborhoodSearch` --- src/general/semidiscretization.jl | 18 ++++- src/neighborhood_search/grid_nhs.jl | 4 + src/neighborhood_search/neighbor_list_nhs.jl | 75 +++++++++++++++++++ .../neighborhood_search.jl | 1 + 4 files changed, 97 insertions(+), 1 deletion(-) create mode 100644 src/neighborhood_search/neighbor_list_nhs.jl diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 0ed168cc3..04f935e5c 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -101,6 +101,21 @@ function create_neighborhood_search(system, neighbor, ::Val{GridNeighborhoodSear return search end +function create_neighborhood_search(system, neighbor, ::Val{NeighborListNeighborhoodSearch}, + min_corner, max_corner) + radius = compact_support(system, neighbor) + grid_nhs = GridNeighborhoodSearch{ndims(system)}(radius, nparticles(neighbor), + min_corner=min_corner, + max_corner=max_corner) + + search = NeighborListNeighborhoodSearch(grid_nhs, nparticles(system)) + + # Initialize neighborhood search + initialize!(search, initial_coordinates(system), initial_coordinates(neighbor)) + + return search +end + @inline function compact_support(system, neighbor) (; smoothing_kernel, smoothing_length) = system return compact_support(smoothing_kernel, smoothing_length) @@ -343,10 +358,11 @@ function update_nhs(u_ode, semi) # Update NHS for each pair of systems foreach_enumerate(systems) do (system_index, system) foreach_enumerate(systems) do (neighbor_index, neighbor) + u = wrap_u(u_ode, system_index, system, semi) u_neighbor = wrap_u(u_ode, neighbor_index, neighbor, semi) neighborhood_search = neighborhood_searches[system_index][neighbor_index] - update!(neighborhood_search, nhs_coords(system, neighbor, u_neighbor)) + update!(neighborhood_search, current_coordinates(u, system), nhs_coords(system, neighbor, u_neighbor)) end end end diff --git a/src/neighborhood_search/grid_nhs.jl b/src/neighborhood_search/grid_nhs.jl index 81fc280fc..2b9825eae 100644 --- a/src/neighborhood_search/grid_nhs.jl +++ b/src/neighborhood_search/grid_nhs.jl @@ -134,6 +134,10 @@ function update!(neighborhood_search::GridNeighborhoodSearch{NDIMS}, update!(neighborhood_search, i -> extract_svector(x, Val(NDIMS), i)) end +function update!(neighborhood_search::GridNeighborhoodSearch, x, y) + update!(neighborhood_search, y) +end + # Modify the existing hash table by moving particles into their new cells function update!(neighborhood_search::GridNeighborhoodSearch, coords_fun) (; hashtable, cell_buffer, cell_buffer_indices) = neighborhood_search diff --git a/src/neighborhood_search/neighbor_list_nhs.jl b/src/neighborhood_search/neighbor_list_nhs.jl new file mode 100644 index 000000000..2830741c6 --- /dev/null +++ b/src/neighborhood_search/neighbor_list_nhs.jl @@ -0,0 +1,75 @@ +struct NeighborListNeighborhoodSearch{ELTYPE, NHS, PB} + search_radius :: ELTYPE + periodic_box :: PB + grid_nhs :: NHS + neighbor_lists :: Vector{Vector{Int}} + + function NeighborListNeighborhoodSearch(grid_nhs, n_particles) + (; search_radius, periodic_box) = grid_nhs + + neighbor_lists = Vector{Vector{Int}}(undef, n_particles) + + new{typeof(search_radius), + typeof(grid_nhs), typeof(periodic_box)}(search_radius, periodic_box, + grid_nhs, neighbor_lists) + end +end + +@inline function Base.ndims(neighborhood_search::NeighborListNeighborhoodSearch) + return ndims(neighborhood_search.grid_nhs) +end + +function initialize!(search::NeighborListNeighborhoodSearch, coords, neighbor_coords) + initialize!(search.grid_nhs, neighbor_coords) + + build_neighbor_lists!(search, coords) +end + +function update!(search::NeighborListNeighborhoodSearch, coords, neighbor_coords) + update!(search.grid_nhs, neighbor_coords) + + build_neighbor_lists!(search, coords) +end + +@inline function eachneighbor(particle, particle_coords, search::NeighborListNeighborhoodSearch) + return search.neighbor_lists[particle] +end + +function build_neighbor_lists!(search::NeighborListNeighborhoodSearch, coords) + (; grid_nhs, neighbor_lists) = search + + for particle in eachindex(neighbor_lists) + particle_coords = extract_svector(coords, Val(ndims(search)), particle) + neighbor_lists[particle] = collect(eachneighbor(particle_coords, grid_nhs)) + end + + return search +end + +@inline function for_particle_neighbor_inner(f, system_coords, neighbor_system_coords, + neighborhood_search::NeighborListNeighborhoodSearch, + particle) + (; search_radius, periodic_box) = neighborhood_search + + particle_coords = extract_svector(system_coords, Val(ndims(neighborhood_search)), + particle) + + for neighbor in eachneighbor(particle, particle_coords, neighborhood_search) + neighbor_coords = extract_svector(neighbor_system_coords, + Val(ndims(neighborhood_search)), neighbor) + + pos_diff = particle_coords - neighbor_coords + distance2 = dot(pos_diff, pos_diff) + + pos_diff, distance2 = compute_periodic_distance(pos_diff, distance2, search_radius, + periodic_box) + + if distance2 <= search_radius^2 + distance = sqrt(distance2) + + # Inline to avoid loss of performance + # compared to not using `for_particle_neighbor`. + @inline f(particle, neighbor, pos_diff, distance) + end + end +end diff --git a/src/neighborhood_search/neighborhood_search.jl b/src/neighborhood_search/neighborhood_search.jl index 0b378454e..623e6dff9 100644 --- a/src/neighborhood_search/neighborhood_search.jl +++ b/src/neighborhood_search/neighborhood_search.jl @@ -108,3 +108,4 @@ end include("trivial_nhs.jl") include("grid_nhs.jl") +include("neighbor_list_nhs.jl") From fe91de93e6985e52bff425a01b15fc1296905818 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Thu, 9 Nov 2023 22:39:46 +0100 Subject: [PATCH 2/4] Only consider actual neighbors --- src/neighborhood_search/neighbor_list_nhs.jl | 38 ++++++++++++++------ 1 file changed, 28 insertions(+), 10 deletions(-) diff --git a/src/neighborhood_search/neighbor_list_nhs.jl b/src/neighborhood_search/neighbor_list_nhs.jl index 2830741c6..29b3e1525 100644 --- a/src/neighborhood_search/neighbor_list_nhs.jl +++ b/src/neighborhood_search/neighbor_list_nhs.jl @@ -22,25 +22,44 @@ end function initialize!(search::NeighborListNeighborhoodSearch, coords, neighbor_coords) initialize!(search.grid_nhs, neighbor_coords) - build_neighbor_lists!(search, coords) + build_neighbor_lists!(search, coords, neighbor_coords) end function update!(search::NeighborListNeighborhoodSearch, coords, neighbor_coords) update!(search.grid_nhs, neighbor_coords) - build_neighbor_lists!(search, coords) + build_neighbor_lists!(search, coords, neighbor_coords) end -@inline function eachneighbor(particle, particle_coords, search::NeighborListNeighborhoodSearch) +@inline function eachneighbor(particle, search::NeighborListNeighborhoodSearch) return search.neighbor_lists[particle] end -function build_neighbor_lists!(search::NeighborListNeighborhoodSearch, coords) - (; grid_nhs, neighbor_lists) = search +function build_neighbor_lists!(search::NeighborListNeighborhoodSearch, coords, + neighbor_coords) + (; grid_nhs, neighbor_lists, search_radius, periodic_box) = search for particle in eachindex(neighbor_lists) particle_coords = extract_svector(coords, Val(ndims(search)), particle) - neighbor_lists[particle] = collect(eachneighbor(particle_coords, grid_nhs)) + + neighbors = [] + + for neighbor in eachneighbor(particle_coords, grid_nhs) + neighbor_particle_coords = extract_svector(neighbor_coords, + Val(ndims(search)), neighbor) + + pos_diff = particle_coords - neighbor_particle_coords + distance2 = dot(pos_diff, pos_diff) + + pos_diff, distance2 = compute_periodic_distance(pos_diff, distance2, + search_radius, periodic_box) + + if distance2 <= search_radius^2 + append!(neighbors, neighbor) + end + end + + neighbor_lists[particle] = neighbors end return search @@ -51,10 +70,9 @@ end particle) (; search_radius, periodic_box) = neighborhood_search - particle_coords = extract_svector(system_coords, Val(ndims(neighborhood_search)), - particle) - - for neighbor in eachneighbor(particle, particle_coords, neighborhood_search) + for neighbor in eachneighbor(particle, neighborhood_search) + particle_coords = extract_svector(system_coords, Val(ndims(neighborhood_search)), + particle) neighbor_coords = extract_svector(neighbor_system_coords, Val(ndims(neighborhood_search)), neighbor) From 5cb4b004c367b062d062dfe1043b2f5d0a84b1ca Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Mon, 13 Nov 2023 12:14:52 +0100 Subject: [PATCH 3/4] Use contiguous vector as neighbor list --- src/neighborhood_search/neighbor_list_nhs.jl | 51 +++++++++++--------- 1 file changed, 29 insertions(+), 22 deletions(-) diff --git a/src/neighborhood_search/neighbor_list_nhs.jl b/src/neighborhood_search/neighbor_list_nhs.jl index 29b3e1525..82b772d47 100644 --- a/src/neighborhood_search/neighbor_list_nhs.jl +++ b/src/neighborhood_search/neighbor_list_nhs.jl @@ -1,17 +1,20 @@ struct NeighborListNeighborhoodSearch{ELTYPE, NHS, PB} - search_radius :: ELTYPE - periodic_box :: PB - grid_nhs :: NHS - neighbor_lists :: Vector{Vector{Int}} + search_radius :: ELTYPE + periodic_box :: PB + grid_nhs :: NHS + neighbor_lists :: Vector{Int} + neighbor_list_start :: Vector{Int} function NeighborListNeighborhoodSearch(grid_nhs, n_particles) (; search_radius, periodic_box) = grid_nhs - neighbor_lists = Vector{Vector{Int}}(undef, n_particles) + neighbor_lists = Int[] + neighbor_list_start = zeros(Int, n_particles + 1) new{typeof(search_radius), typeof(grid_nhs), typeof(periodic_box)}(search_radius, periodic_box, - grid_nhs, neighbor_lists) + grid_nhs, neighbor_lists, + neighbor_list_start) end end @@ -32,36 +35,40 @@ function update!(search::NeighborListNeighborhoodSearch, coords, neighbor_coords end @inline function eachneighbor(particle, search::NeighborListNeighborhoodSearch) - return search.neighbor_lists[particle] + (; neighbor_lists, neighbor_list_start) = search + + return (neighbor_lists[i] for i in neighbor_list_start[particle]:(neighbor_list_start[particle + 1] - 1)) end function build_neighbor_lists!(search::NeighborListNeighborhoodSearch, coords, neighbor_coords) - (; grid_nhs, neighbor_lists, search_radius, periodic_box) = search + (; grid_nhs, neighbor_lists, neighbor_list_start, search_radius, periodic_box) = search - for particle in eachindex(neighbor_lists) - particle_coords = extract_svector(coords, Val(ndims(search)), particle) + resize!(neighbor_lists, 0) + + for particle in 1:(length(neighbor_list_start) - 1) + neighbor_list_start[particle] = length(neighbor_lists) + 1 - neighbors = [] + particle_coords = extract_svector(coords, Val(ndims(search)), particle) for neighbor in eachneighbor(particle_coords, grid_nhs) - neighbor_particle_coords = extract_svector(neighbor_coords, - Val(ndims(search)), neighbor) + # neighbor_particle_coords = extract_svector(neighbor_coords, + # Val(ndims(search)), neighbor) - pos_diff = particle_coords - neighbor_particle_coords - distance2 = dot(pos_diff, pos_diff) + # pos_diff = particle_coords - neighbor_particle_coords + # distance2 = dot(pos_diff, pos_diff) - pos_diff, distance2 = compute_periodic_distance(pos_diff, distance2, - search_radius, periodic_box) + # pos_diff, distance2 = compute_periodic_distance(pos_diff, distance2, + # search_radius, periodic_box) - if distance2 <= search_radius^2 - append!(neighbors, neighbor) - end + # if distance2 <= search_radius^2 + append!(neighbor_lists, neighbor) + # end end - - neighbor_lists[particle] = neighbors end + neighbor_list_start[end] = length(neighbor_lists) + 1 + return search end From 83f1e19e38e9d055fe4c9a4a28eb2c8d11377ad4 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Mon, 13 Nov 2023 14:36:07 +0100 Subject: [PATCH 4/4] [skip ci] Reformat --- src/general/semidiscretization.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/general/semidiscretization.jl b/src/general/semidiscretization.jl index 04f935e5c..ef3f455ca 100644 --- a/src/general/semidiscretization.jl +++ b/src/general/semidiscretization.jl @@ -362,7 +362,8 @@ function update_nhs(u_ode, semi) u_neighbor = wrap_u(u_ode, neighbor_index, neighbor, semi) neighborhood_search = neighborhood_searches[system_index][neighbor_index] - update!(neighborhood_search, current_coordinates(u, system), nhs_coords(system, neighbor, u_neighbor)) + update!(neighborhood_search, current_coordinates(u, system), + nhs_coords(system, neighbor, u_neighbor)) end end end