From a2d9f64a7dde587bca635891193c4f3dff0118f9 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Fri, 21 Jun 2024 12:38:01 +0200 Subject: [PATCH 1/4] Replace all references to "particle" by "point" (#31) * Replace all references to "particle" by "point" * Reformat * Fix typo --- src/PointNeighbors.jl | 2 +- src/cell_lists/dictionary.jl | 6 +-- src/neighborhood_search.jl | 50 ++++++++++---------- src/nhs_grid.jl | 92 ++++++++++++++++++------------------ src/nhs_neighbor_lists.jl | 33 +++++++------ src/nhs_trivial.jl | 20 ++++---- test/neighborhood_search.jl | 40 ++++++++-------- test/nhs_grid.jl | 64 ++++++++++++------------- test/nhs_trivial.jl | 2 +- 9 files changed, 150 insertions(+), 159 deletions(-) diff --git a/src/PointNeighbors.jl b/src/PointNeighbors.jl index 520099c2..f8e4de39 100644 --- a/src/PointNeighbors.jl +++ b/src/PointNeighbors.jl @@ -13,7 +13,7 @@ include("cell_lists/cell_lists.jl") include("nhs_grid.jl") include("nhs_neighbor_lists.jl") -export for_particle_neighbor, foreach_neighbor +export foreach_point_neighbor, foreach_neighbor export TrivialNeighborhoodSearch, GridNeighborhoodSearch, PrecomputedNeighborhoodSearch export initialize!, update!, initialize_grid!, update_grid! diff --git a/src/cell_lists/dictionary.jl b/src/cell_lists/dictionary.jl index b65b59dd..d872160b 100644 --- a/src/cell_lists/dictionary.jl +++ b/src/cell_lists/dictionary.jl @@ -16,13 +16,13 @@ function Base.empty!(cell_list::DictionaryCellList) return cell_list end -function push_cell!(cell_list::DictionaryCellList, cell, particle) +function push_cell!(cell_list::DictionaryCellList, cell, point) (; hashtable) = cell_list if haskey(hashtable, cell) - append!(hashtable[cell], particle) + append!(hashtable[cell], point) else - hashtable[cell] = [particle] + hashtable[cell] = [point] end return cell_list diff --git a/src/neighborhood_search.jl b/src/neighborhood_search.jl index 02ccc53b..87ab1296 100644 --- a/src/neighborhood_search.jl +++ b/src/neighborhood_search.jl @@ -15,7 +15,7 @@ See also [`update!`](@ref). @inline initialize!(search::AbstractNeighborhoodSearch, x, y) = search """ - update!(search::AbstractNeighborhoodSearch, x, y; particles_moving = (true, true)) + update!(search::AbstractNeighborhoodSearch, x, y; points_moving = (true, true)) Update an already initialized neighborhood search with the two coordinate arrays `x` and `y`. @@ -29,15 +29,15 @@ If incremental updates are not possible for an implementation, `update!` will fa to a regular `initialize!`. Some neighborhood searches might not need to update when only `x` changed since the last -`update!` or `initialize!` and `y` did not change. Pass `particles_moving = (true, false)` +`update!` or `initialize!` and `y` did not change. Pass `points_moving = (true, false)` in this case to avoid unnecessary updates. -The first flag in `particles_moving` indicates if points in `x` are moving. +The first flag in `points_moving` indicates if points in `x` are moving. The second flag indicates if points in `y` are moving. See also [`initialize!`](@ref). """ @inline function update!(search::AbstractNeighborhoodSearch, x, y; - particles_moving = (true, true)) + points_moving = (true, true)) return search end @@ -56,43 +56,42 @@ end # Otherwise, unspecialized code will cause a lot of allocations # and heavily impact performance. # See https://docs.julialang.org/en/v1/manual/performance-tips/#Be-aware-of-when-Julia-avoids-specializing -function for_particle_neighbor(f::T, system_coords, neighbor_coords, neighborhood_search; - particles = axes(system_coords, 2), - parallel = true) where {T} - for_particle_neighbor(f, system_coords, neighbor_coords, neighborhood_search, particles, - Val(parallel)) +function foreach_point_neighbor(f::T, system_coords, neighbor_coords, neighborhood_search; + points = axes(system_coords, 2), + parallel = true) where {T} + foreach_point_neighbor(f, system_coords, neighbor_coords, neighborhood_search, points, + Val(parallel)) end -@inline function for_particle_neighbor(f, system_coords, neighbor_coords, - neighborhood_search, particles, parallel::Val{true}) - @threaded for particle in particles - foreach_neighbor(f, system_coords, neighbor_coords, neighborhood_search, particle) +@inline function foreach_point_neighbor(f, system_coords, neighbor_coords, + neighborhood_search, points, parallel::Val{true}) + @threaded for point in points + foreach_neighbor(f, system_coords, neighbor_coords, neighborhood_search, point) end return nothing end -@inline function for_particle_neighbor(f, system_coords, neighbor_coords, - neighborhood_search, particles, parallel::Val{false}) - for particle in particles - foreach_neighbor(f, system_coords, neighbor_coords, neighborhood_search, particle) +@inline function foreach_point_neighbor(f, system_coords, neighbor_coords, + neighborhood_search, points, parallel::Val{false}) + for point in points + foreach_neighbor(f, system_coords, neighbor_coords, neighborhood_search, point) end return nothing end @inline function foreach_neighbor(f, system_coords, neighbor_system_coords, - neighborhood_search, particle; + neighborhood_search, point; search_radius = neighborhood_search.search_radius) (; periodic_box) = neighborhood_search - particle_coords = extract_svector(system_coords, Val(ndims(neighborhood_search)), - particle) - for neighbor in eachneighbor(particle_coords, neighborhood_search) + point_coords = extract_svector(system_coords, Val(ndims(neighborhood_search)), point) + for neighbor in eachneighbor(point_coords, neighborhood_search) neighbor_coords = extract_svector(neighbor_system_coords, Val(ndims(neighborhood_search)), neighbor) - pos_diff = particle_coords - neighbor_coords + pos_diff = point_coords - neighbor_coords distance2 = dot(pos_diff, pos_diff) pos_diff, distance2 = compute_periodic_distance(pos_diff, distance2, search_radius, @@ -102,8 +101,8 @@ end distance = sqrt(distance2) # Inline to avoid loss of performance - # compared to not using `for_particle_neighbor`. - @inline f(particle, neighbor, pos_diff, distance) + # compared to not using `foreach_point_neighbor`. + @inline f(point, neighbor, pos_diff, distance) end end end @@ -113,8 +112,7 @@ end return pos_diff, distance2 end -@inline function compute_periodic_distance(pos_diff, distance2, search_radius, - periodic_box) +@inline function compute_periodic_distance(pos_diff, distance2, search_radius, periodic_box) if distance2 > search_radius^2 # Use periodic `pos_diff` pos_diff -= periodic_box.size .* round.(pos_diff ./ periodic_box.size) diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index 038eb661..af4b690d 100644 --- a/src/nhs_grid.jl +++ b/src/nhs_grid.jl @@ -1,10 +1,10 @@ @doc raw""" - GridNeighborhoodSearch{NDIMS}(search_radius, n_particles; periodic_box_min_corner=nothing, + GridNeighborhoodSearch{NDIMS}(search_radius, n_points; periodic_box_min_corner=nothing, periodic_box_max_corner=nothing, threaded_nhs_update=true) Simple grid-based neighborhood search with uniform search radius. The domain is divided into a regular grid. -For each (non-empty) grid cell, a list of particles in this cell is stored. +For each (non-empty) grid cell, a list of points in this cell is stored. Instead of representing a finite domain by an array of cells, a potentially infinite domain is represented by storing cell lists in a hash table (using Julia's `Dict` data structure), indexed by the cell index tuple @@ -14,18 +14,18 @@ indexed by the cell index tuple ``` where ``x, y, z`` are the space coordinates and ``d`` is the search radius. -To find particles within the search radius around a point, only particles in the neighboring +To find points within the search radius around a position, only points in the neighboring cells are considered. See also (Chalela et al., 2021), (Ihmsen et al. 2011, Section 4.4). -As opposed to (Ihmsen et al. 2011), we do not sort the particles in any way, +As opposed to (Ihmsen et al. 2011), we do not sort the points in any way, since not sorting makes our implementation a lot faster (although less parallelizable). # Arguments - `NDIMS`: Number of dimensions. - `search_radius`: The uniform search radius. -- `n_particles`: Total number of particles. +- `n_points`: Total number of points. # Keywords - `periodic_box_min_corner`: In order to use a (rectangular) periodic domain, pass the @@ -36,7 +36,7 @@ since not sorting makes our implementation a lot faster (although less paralleli directions. - `threaded_nhs_update=true`: Can be used to deactivate thread parallelization in the neighborhood search update. This can be one of the largest sources of variations between simulations - with different thread numbers due to particle ordering changes. + with different thread numbers due to neighbor ordering changes. ## References - M. Chalela, E. Sillero, L. Pereyra, M.A. Garcia, J.B. Cabral, M. Lares, M. Merchán. @@ -58,14 +58,14 @@ struct GridNeighborhoodSearch{NDIMS, ELTYPE, CL, PB} <: AbstractNeighborhoodSear cell_buffer_indices :: Vector{Int} # Store which entries of `cell_buffer` are initialized threaded_nhs_update :: Bool - function GridNeighborhoodSearch{NDIMS}(search_radius, n_particles; + function GridNeighborhoodSearch{NDIMS}(search_radius, n_points; periodic_box_min_corner = nothing, periodic_box_max_corner = nothing, threaded_nhs_update = true) where {NDIMS} ELTYPE = typeof(search_radius) cell_list = DictionaryCellList{NDIMS}() - cell_buffer = Array{NTuple{NDIMS, Int}, 2}(undef, n_particles, Threads.nthreads()) + cell_buffer = Array{NTuple{NDIMS, Int}, 2}(undef, n_points, Threads.nthreads()) cell_buffer_indices = zeros(Int, Threads.nthreads()) if search_radius < eps() || @@ -105,7 +105,7 @@ end @inline Base.ndims(neighborhood_search::GridNeighborhoodSearch{NDIMS}) where {NDIMS} = NDIMS -@inline function nparticles(neighborhood_search::GridNeighborhoodSearch) +@inline function npoints(neighborhood_search::GridNeighborhoodSearch) return size(neighborhood_search.cell_buffer, 1) end @@ -124,12 +124,12 @@ function initialize_grid!(neighborhood_search::GridNeighborhoodSearch, coords_fu empty!(cell_list) - for particle in 1:nparticles(neighborhood_search) - # Get cell index of the particle's cell - cell = cell_coords(coords_fun(particle), neighborhood_search) + for point in 1:npoints(neighborhood_search) + # Get cell index of the point's cell + cell = cell_coords(coords_fun(point), neighborhood_search) - # Add particle to corresponding cell - push_cell!(cell_list, cell, particle) + # Add point to corresponding cell + push_cell!(cell_list, cell, point) end return neighborhood_search @@ -137,10 +137,10 @@ end function update!(neighborhood_search::GridNeighborhoodSearch, x::AbstractMatrix, y::AbstractMatrix; - particles_moving = (true, true)) - # The coordinates of the first set of particles are irrelevant for this NHS. + points_moving = (true, true)) + # The coordinates of the first set of points are irrelevant for this NHS. # Only update when the second set is moving. - particles_moving[2] || return neighborhood_search + points_moving[2] || return neighborhood_search update_grid!(neighborhood_search, y) end @@ -151,40 +151,40 @@ function update_grid!(neighborhood_search::GridNeighborhoodSearch{NDIMS}, update_grid!(neighborhood_search, i -> extract_svector(y, Val(NDIMS), i)) end -# Modify the existing hash table by moving particles into their new cells +# Modify the existing hash table by moving points into their new cells function update_grid!(neighborhood_search::GridNeighborhoodSearch, coords_fun) (; cell_list, cell_buffer, cell_buffer_indices, threaded_nhs_update) = neighborhood_search # Reset `cell_buffer` by moving all pointers to the beginning cell_buffer_indices .= 0 - # Find all cells containing particles that now belong to another cell + # Find all cells containing points that now belong to another cell mark_changed_cell!(neighborhood_search, cell_list, coords_fun, Val(threaded_nhs_update)) - # Iterate over all marked cells and move the particles into their new cells. + # Iterate over all marked cells and move the points into their new cells. for thread in 1:Threads.nthreads() # Only the entries `1:cell_buffer_indices[thread]` are initialized for `thread`. for i in 1:cell_buffer_indices[thread] cell = cell_buffer[i, thread] - particles = cell_list[cell] + points = cell_list[cell] - # Find all particles whose coordinates do not match this cell - moved_particle_indices = (i for i in eachindex(particles) - if cell_coords(coords_fun(particles[i]), - neighborhood_search) != cell) + # Find all points whose coordinates do not match this cell + moved_point_indices = (i for i in eachindex(points) + if cell_coords(coords_fun(points[i]), + neighborhood_search) != cell) - # Add moved particles to new cell - for i in moved_particle_indices - particle = particles[i] - new_cell_coords = cell_coords(coords_fun(particle), neighborhood_search) + # Add moved points to new cell + for i in moved_point_indices + point = points[i] + new_cell_coords = cell_coords(coords_fun(point), neighborhood_search) - # Add particle to corresponding cell or create cell if it does not exist - push_cell!(cell_list, new_cell_coords, particle) + # Add point to corresponding cell or create cell if it does not exist + push_cell!(cell_list, new_cell_coords, point) end - # Remove moved particles from this cell - deleteat_cell!(cell_list, cell, moved_particle_indices) + # Remove moved points from this cell + deleteat_cell!(cell_list, cell, moved_point_indices) end end @@ -213,8 +213,8 @@ end @inline function mark_changed_cell!(neighborhood_search, cell, coords_fun) (; cell_list, cell_buffer, cell_buffer_indices) = neighborhood_search - for particle in cell_list[cell] - if cell_coords(coords_fun(particle), neighborhood_search) != cell + for point in cell_list[cell] + if cell_coords(coords_fun(point), neighborhood_search) != cell # Mark this cell and continue with the next one. # # `cell_buffer` is preallocated, @@ -233,8 +233,8 @@ end # Generator of all neighboring cells to consider neighboring_cells = ((x + i) for i in -1:1) - # Merge all lists of particles in the neighboring cells into one iterator - Iterators.flatten(particles_in_cell(cell, neighborhood_search) + # Merge all lists of points in the neighboring cells into one iterator + Iterators.flatten(points_in_cell(cell, neighborhood_search) for cell in neighboring_cells) end @@ -245,8 +245,8 @@ end # Generator of all neighboring cells to consider neighboring_cells = ((x + i, y + j) for i in -1:1, j in -1:1) - # Merge all lists of particles in the neighboring cells into one iterator - Iterators.flatten(particles_in_cell(cell, neighborhood_search) + # Merge all lists of points in the neighboring cells into one iterator + Iterators.flatten(points_in_cell(cell, neighborhood_search) for cell in neighboring_cells) end @@ -257,12 +257,12 @@ end # Generator of all neighboring cells to consider neighboring_cells = ((x + i, y + j, z + k) for i in -1:1, j in -1:1, k in -1:1) - # Merge all lists of particles in the neighboring cells into one iterator - Iterators.flatten(particles_in_cell(cell, neighborhood_search) + # Merge all lists of points in the neighboring cells into one iterator + Iterators.flatten(points_in_cell(cell, neighborhood_search) for cell in neighboring_cells) end -@inline function particles_in_cell(cell_index, neighborhood_search) +@inline function points_in_cell(cell_index, neighborhood_search) (; cell_list) = neighborhood_search return cell_list[periodic_cell_index(cell_index, neighborhood_search)] @@ -300,9 +300,9 @@ end return Tuple(floor_to_int.(offset_coords ./ cell_size)) end -# When particles end up with coordinates so big that the cell coordinates +# When points end up with coordinates so big that the cell coordinates # exceed the range of Int, then `floor(Int, i)` will fail with an InexactError. -# In this case, we can just use typemax(Int), since we can assume that particles +# In this case, we can just use typemax(Int), since we can assume that points # that far away will not interact with anything, anyway. # This usually indicates an instability, but we don't want the simulation to crash, # since adaptive time integration methods may detect the instability and reject the @@ -322,9 +322,9 @@ end # Create a copy of a neighborhood search but with a different search radius function copy_neighborhood_search(nhs::GridNeighborhoodSearch, search_radius, x, y) if nhs.periodic_box === nothing - search = GridNeighborhoodSearch{ndims(nhs)}(search_radius, nparticles(nhs)) + search = GridNeighborhoodSearch{ndims(nhs)}(search_radius, npoints(nhs)) else - search = GridNeighborhoodSearch{ndims(nhs)}(search_radius, nparticles(nhs), + search = GridNeighborhoodSearch{ndims(nhs)}(search_radius, npoints(nhs), periodic_box_min_corner = nhs.periodic_box.min_corner, periodic_box_max_corner = nhs.periodic_box.max_corner) end diff --git a/src/nhs_neighbor_lists.jl b/src/nhs_neighbor_lists.jl index 244df129..07854b65 100644 --- a/src/nhs_neighbor_lists.jl +++ b/src/nhs_neighbor_lists.jl @@ -1,10 +1,10 @@ @doc raw""" - PrecomputedNeighborhoodSearch{NDIMS}(search_radius, n_particles; + PrecomputedNeighborhoodSearch{NDIMS}(search_radius, n_points; periodic_box_min_corner = nothing, periodic_box_max_corner = nothing) Neighborhood search with precomputed neighbor lists. A list of all neighbors is computed -for each particle during initialization and update. +for each point during initialization and update. This neighborhood search maximizes the performance of neighbor loops at the cost of a much slower [`update!`](@ref). @@ -14,7 +14,7 @@ initialization and update. # Arguments - `NDIMS`: Number of dimensions. - `search_radius`: The uniform search radius. -- `n_particles`: Total number of particles. +- `n_points`: Total number of points. # Keywords - `periodic_box_min_corner`: In order to use a (rectangular) periodic domain, pass the @@ -29,12 +29,12 @@ struct PrecomputedNeighborhoodSearch{NDIMS, NHS, NL, PB} neighbor_lists :: NL periodic_box :: PB - function PrecomputedNeighborhoodSearch{NDIMS}(search_radius, n_particles; + function PrecomputedNeighborhoodSearch{NDIMS}(search_radius, n_points; periodic_box_min_corner = nothing, periodic_box_max_corner = nothing) where { NDIMS } - nhs = GridNeighborhoodSearch{NDIMS}(search_radius, n_particles, + nhs = GridNeighborhoodSearch{NDIMS}(search_radius, n_points, periodic_box_min_corner = periodic_box_min_corner, periodic_box_max_corner = periodic_box_max_corner) @@ -62,14 +62,14 @@ end function update!(search::PrecomputedNeighborhoodSearch, x::AbstractMatrix, y::AbstractMatrix; - particles_moving = (true, true)) + points_moving = (true, true)) (; neighborhood_search, neighbor_lists) = search # Update grid NHS - update!(neighborhood_search, x, y, particles_moving = particles_moving) + update!(neighborhood_search, x, y, points_moving = points_moving) # Skip update if both point sets are static - if any(particles_moving) + if any(points_moving) initialize_neighbor_lists!(neighbor_lists, neighborhood_search, x, y) end end @@ -83,24 +83,23 @@ function initialize_neighbor_lists!(neighbor_lists, neighborhood_search, x, y) end # Fill neighbor lists - for_particle_neighbor(x, y, neighborhood_search) do particle, neighbor, _, _ - push!(neighbor_lists[particle], neighbor) + foreach_point_neighbor(x, y, neighborhood_search) do point, neighbor, _, _ + push!(neighbor_lists[point], neighbor) end end @inline function foreach_neighbor(f, system_coords, neighbor_system_coords, neighborhood_search::PrecomputedNeighborhoodSearch, - particle; search_radius = nothing) + point; search_radius = nothing) (; periodic_box, neighbor_lists) = neighborhood_search (; search_radius) = neighborhood_search.neighborhood_search - particle_coords = extract_svector(system_coords, Val(ndims(neighborhood_search)), - particle) - for neighbor in neighbor_lists[particle] + point_coords = extract_svector(system_coords, Val(ndims(neighborhood_search)), point) + for neighbor in neighbor_lists[point] neighbor_coords = extract_svector(neighbor_system_coords, Val(ndims(neighborhood_search)), neighbor) - pos_diff = particle_coords - neighbor_coords + pos_diff = point_coords - neighbor_coords distance2 = dot(pos_diff, pos_diff) pos_diff, distance2 = compute_periodic_distance(pos_diff, distance2, search_radius, @@ -109,7 +108,7 @@ end distance = sqrt(distance2) # Inline to avoid loss of performance - # compared to not using `for_particle_neighbor`. - @inline f(particle, neighbor, pos_diff, distance) + # compared to not using `foreach_point_neighbor`. + @inline f(point, neighbor, pos_diff, distance) end end diff --git a/src/nhs_trivial.jl b/src/nhs_trivial.jl index 699b663d..d99091e3 100644 --- a/src/nhs_trivial.jl +++ b/src/nhs_trivial.jl @@ -1,15 +1,12 @@ @doc raw""" - TrivialNeighborhoodSearch{NDIMS}(search_radius, eachparticle) + TrivialNeighborhoodSearch{NDIMS}(search_radius, eachpoint) -Trivial neighborhood search that simply loops over all particles. -The search radius still needs to be passed in order to sort out particles outside the -search radius in the internal function `for_particle_neighbor`, but it's not used in the -internal function `eachneighbor`. +Trivial neighborhood search that simply loops over all points. # Arguments - `NDIMS`: Number of dimensions. - `search_radius`: The uniform search radius. -- `eachparticle`: `UnitRange` of all particle indices. Usually just `1:n_particles`. +- `eachpoint`: `UnitRange` of all point indices. Usually just `1:n_points`. # Keywords - `periodic_box_min_corner`: In order to use a (rectangular) periodic domain, pass the @@ -21,10 +18,10 @@ internal function `eachneighbor`. """ struct TrivialNeighborhoodSearch{NDIMS, ELTYPE, EP, PB} <: AbstractNeighborhoodSearch search_radius :: ELTYPE - eachparticle :: EP + eachpoint :: EP periodic_box :: PB - function TrivialNeighborhoodSearch{NDIMS}(search_radius, eachparticle; + function TrivialNeighborhoodSearch{NDIMS}(search_radius, eachpoint; periodic_box_min_corner = nothing, periodic_box_max_corner = nothing) where { NDIMS @@ -42,8 +39,7 @@ struct TrivialNeighborhoodSearch{NDIMS, ELTYPE, EP, PB} <: AbstractNeighborhoodS end new{NDIMS, typeof(search_radius), - typeof(eachparticle), typeof(periodic_box)}(search_radius, eachparticle, - periodic_box) + typeof(eachpoint), typeof(periodic_box)}(search_radius, eachpoint, periodic_box) end end @@ -56,11 +52,11 @@ end @inline initialize!(search::TrivialNeighborhoodSearch, x, y) = search @inline function update!(search::TrivialNeighborhoodSearch, x, y; - particles_moving = (true, true)) + points_moving = (true, true)) return search end -@inline eachneighbor(coords, search::TrivialNeighborhoodSearch) = search.eachparticle +@inline eachneighbor(coords, search::TrivialNeighborhoodSearch) = search.eachpoint # Create a copy of a neighborhood search but with a different search radius function copy_neighborhood_search(nhs::TrivialNeighborhoodSearch, search_radius, x, y) diff --git a/test/neighborhood_search.jl b/test/neighborhood_search.jl index 3cc256ba..2027328b 100644 --- a/test/neighborhood_search.jl +++ b/test/neighborhood_search.jl @@ -35,17 +35,17 @@ coords = coordinates[i] NDIMS = size(coords, 1) - n_particles = size(coords, 2) + n_points = size(coords, 2) search_radius = 0.1 neighborhood_searches = [ - TrivialNeighborhoodSearch{NDIMS}(search_radius, 1:n_particles, + TrivialNeighborhoodSearch{NDIMS}(search_radius, 1:n_points, periodic_box_min_corner = periodic_boxes[i][1], periodic_box_max_corner = periodic_boxes[i][2]), - GridNeighborhoodSearch{NDIMS}(search_radius, n_particles, + GridNeighborhoodSearch{NDIMS}(search_radius, n_points, periodic_box_min_corner = periodic_boxes[i][1], periodic_box_max_corner = periodic_boxes[i][2]), - PrecomputedNeighborhoodSearch{NDIMS}(search_radius, n_particles, + PrecomputedNeighborhoodSearch{NDIMS}(search_radius, n_points, periodic_box_min_corner = periodic_boxes[i][1], periodic_box_max_corner = periodic_boxes[i][2]), ] @@ -63,10 +63,10 @@ neighbors = [Int[] for _ in axes(coords, 2)] - for_particle_neighbor(coords, coords, nhs, - particles = axes(coords, 2)) do particle, neighbor, - pos_diff, distance - append!(neighbors[particle], neighbor) + foreach_point_neighbor(coords, coords, nhs, + points = axes(coords, 2)) do point, neighbor, + pos_diff, distance + append!(neighbors[point], neighbor) end # All of these tests are designed to yield the same neighbor lists. @@ -97,28 +97,28 @@ coords = point_cloud(cloud_size, seed = seed) NDIMS = length(cloud_size) - n_particles = size(coords, 2) + n_points = size(coords, 2) search_radius = 2.5 # Use different coordinates for `initialize!` and then `update!` with the # correct coordinates to make sure that `update!` is working as well. coords_initialize = point_cloud(cloud_size, seed = 1) - # Compute expected neighbor lists by brute-force looping over all particles + # Compute expected neighbor lists by brute-force looping over all points # as potential neighbors (`TrivialNeighborhoodSearch`). trivial_nhs = TrivialNeighborhoodSearch{NDIMS}(search_radius, axes(coords, 2)) neighbors_expected = [Int[] for _ in axes(coords, 2)] - for_particle_neighbor(coords, coords, trivial_nhs, - parallel = false) do particle, neighbor, - pos_diff, distance - append!(neighbors_expected[particle], neighbor) + foreach_point_neighbor(coords, coords, trivial_nhs, + parallel = false) do point, neighbor, + pos_diff, distance + append!(neighbors_expected[point], neighbor) end neighborhood_searches = [ - GridNeighborhoodSearch{NDIMS}(search_radius, n_particles), - PrecomputedNeighborhoodSearch{NDIMS}(search_radius, n_particles), + GridNeighborhoodSearch{NDIMS}(search_radius, n_points), + PrecomputedNeighborhoodSearch{NDIMS}(search_radius, n_points), ] neighborhood_searches_names = [ @@ -141,10 +141,10 @@ neighbors = [Int[] for _ in axes(coords, 2)] - for_particle_neighbor(coords, coords, nhs, - parallel = false) do particle, neighbor, - pos_diff, distance - append!(neighbors[particle], neighbor) + foreach_point_neighbor(coords, coords, nhs, + parallel = false) do point, neighbor, + pos_diff, distance + append!(neighbors[point], neighbor) end @test sort.(neighbors) == neighbors_expected diff --git a/test/nhs_grid.jl b/test/nhs_grid.jl index 90801ccd..68aca46c 100644 --- a/test/nhs_grid.jl +++ b/test/nhs_grid.jl @@ -15,54 +15,53 @@ @testset "Rectangular Point Cloud 2D" begin #### Setup - # Rectangle of equidistantly spaced particles + # Rectangle of equidistantly spaced points # from (x, y) = (-0.25, -0.25) to (x, y) = (0.35, 0.35). range = -0.25:0.1:0.35 coordinates1 = hcat(collect.(Iterators.product(range, range))...) - nparticles = size(coordinates1, 2) + npoints = size(coordinates1, 2) - particle_position1 = [0.05, 0.05] - particle_spacing = 0.1 - radius = particle_spacing + point_position1 = [0.05, 0.05] + radius = 0.1 # Create neighborhood search - nhs1 = GridNeighborhoodSearch{2}(radius, nparticles) + nhs1 = GridNeighborhoodSearch{2}(radius, npoints) initialize_grid!(nhs1, coordinates1) - # Get each neighbor for `particle_position1` - neighbors1 = sort(collect(PointNeighbors.eachneighbor(particle_position1, nhs1))) + # Get each neighbor for `point_position1` + neighbors1 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs1))) - # Move particles + # Move points coordinates2 = coordinates1 .+ [1.4, -3.5] # Update neighborhood search update_grid!(nhs1, coordinates2) # Get each neighbor for updated NHS - neighbors2 = sort(collect(PointNeighbors.eachneighbor(particle_position1, nhs1))) + neighbors2 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs1))) # Change position - particle_position2 = particle_position1 .+ [1.4, -3.5] + point_position2 = point_position1 .+ [1.4, -3.5] - # Get each neighbor for `particle_position2` - neighbors3 = sort(collect(PointNeighbors.eachneighbor(particle_position2, nhs1))) + # Get each neighbor for `point_position2` + neighbors3 = sort(collect(PointNeighbors.eachneighbor(point_position2, nhs1))) # Double search radius nhs2 = GridNeighborhoodSearch{2}(2 * radius, size(coordinates1, 2)) initialize!(nhs2, coordinates1, coordinates1) # Get each neighbor in double search radius - neighbors4 = sort(collect(PointNeighbors.eachneighbor(particle_position1, nhs2))) + neighbors4 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs2))) - # Move particles + # Move points coordinates2 = coordinates1 .+ [0.4, -0.4] # Update neighborhood search update!(nhs2, coordinates2, coordinates2) # Get each neighbor in double search radius - neighbors5 = sort(collect(PointNeighbors.eachneighbor(particle_position1, nhs2))) + neighbors5 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs2))) #### Verification against lists of potential neighbors built by hand @test neighbors1 == [17, 18, 19, 24, 25, 26, 31, 32, 33] @@ -80,26 +79,25 @@ @testset "Rectangular Point Cloud 3D" begin #### Setup - # Rectangle of equidistantly spaced particles + # Rectangle of equidistantly spaced points # from (x, y, z) = (-0.25, -0.25, -0.25) to (x, y, z) = (0.35, 0.35, 0.35). range = -0.25:0.1:0.35 coordinates1 = hcat(collect.(Iterators.product(range, range, range))...) - nparticles = size(coordinates1, 2) + npoints = size(coordinates1, 2) - particle_position1 = [0.05, 0.05, 0.05] - particle_spacing = 0.1 - radius = particle_spacing + point_position1 = [0.05, 0.05, 0.05] + radius = 0.1 # Create neighborhood search - nhs1 = GridNeighborhoodSearch{3}(radius, nparticles) + nhs1 = GridNeighborhoodSearch{3}(radius, npoints) coords_fun(i) = coordinates1[:, i] initialize_grid!(nhs1, coords_fun) - # Get each neighbor for `particle_position1` - neighbors1 = sort(collect(PointNeighbors.eachneighbor(particle_position1, nhs1))) + # Get each neighbor for `point_position1` + neighbors1 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs1))) - # Move particles + # Move points coordinates2 = coordinates1 .+ [1.4, -3.5, 0.8] # Update neighborhood search @@ -107,13 +105,13 @@ update_grid!(nhs1, coords_fun2) # Get each neighbor for updated NHS - neighbors2 = sort(collect(PointNeighbors.eachneighbor(particle_position1, nhs1))) + neighbors2 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs1))) # Change position - particle_position2 = particle_position1 .+ [1.4, -3.5, 0.8] + point_position2 = point_position1 .+ [1.4, -3.5, 0.8] - # Get each neighbor for `particle_position2` - neighbors3 = sort(collect(PointNeighbors.eachneighbor(particle_position2, nhs1))) + # Get each neighbor for `point_position2` + neighbors3 = sort(collect(PointNeighbors.eachneighbor(point_position2, nhs1))) #### Verification against lists of potential neighbors built by hand @test neighbors1 == @@ -129,8 +127,8 @@ @testset verbose=true "Periodicity" begin # These setups are the same as in `test/neighborhood_search.jl`, - # but instead of testing the actual neighbors with `for_particle_neighbor`, - # we only test the potential neighbors (particles in neighboring cells) here. + # but instead of testing the actual neighbors with `foreach_point_neighbor`, + # we only test the potential neighbors (points in neighboring cells) here. # Names, coordinates and corresponding periodic boxes for each test names = [ @@ -180,8 +178,8 @@ @testset "Offset Domain Triggering Split Cells" begin # This test used to trigger a "split cell bug", where the left and right # boundary cells were only partially contained in the domain. - # The left particle was placed inside a ghost cells, which caused it to not - # see the right particle, even though it was within the search distance. + # The left point was placed inside a ghost cells, which caused it to not + # see the right point, even though it was within the search distance. # The domain size is an integer multiple of the cell size, but the NHS did not # offset the grid based on the domain position. # See https://github.com/trixi-framework/TrixiParticles.jl/pull/211 diff --git a/test/nhs_trivial.jl b/test/nhs_trivial.jl index fb514da1..096dbf0c 100644 --- a/test/nhs_trivial.jl +++ b/test/nhs_trivial.jl @@ -1,5 +1,5 @@ @testset verbose=true "TrivialNeighborhoodSearch" begin - # Setup with 5 particles + # Setup with 5 points nhs = TrivialNeighborhoodSearch{2}(1.0, Base.OneTo(5)) # Get each neighbor for arbitrary coordinates From fbee028f4c763b9dd9eb5d8fa1cb054995f785af Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Fri, 21 Jun 2024 14:06:29 +0200 Subject: [PATCH 2/4] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index a3270299..fecbf1cc 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "PointNeighbors" uuid = "1c4d5385-0a27-49de-8e2c-43b175c8985c" authors = ["Erik Faulhaber "] -version = "0.2.4-dev" +version = "0.3.0" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" From 021f57f9bdb90bbcbd442c8a385b4706f007d3a8 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Fri, 21 Jun 2024 15:54:50 +0200 Subject: [PATCH 3/4] Revise periodic box API (#30) * Revise periodic box API * Allow `Vector`s in constructor of `PeriodicBox` * Fix tests * Fix docs for `TrivialNeighborhoodSearch` --- src/PointNeighbors.jl | 2 +- src/neighborhood_search.jl | 18 ++++++++++++--- src/nhs_grid.jl | 45 ++++++++++++------------------------- src/nhs_neighbor_lists.jl | 21 +++++------------ src/nhs_trivial.jl | 35 +++++------------------------ test/neighborhood_search.jl | 15 +++++-------- test/nhs_grid.jl | 22 ++++++++++-------- 7 files changed, 61 insertions(+), 97 deletions(-) diff --git a/src/PointNeighbors.jl b/src/PointNeighbors.jl index f8e4de39..68c0e41c 100644 --- a/src/PointNeighbors.jl +++ b/src/PointNeighbors.jl @@ -13,7 +13,7 @@ include("cell_lists/cell_lists.jl") include("nhs_grid.jl") include("nhs_neighbor_lists.jl") -export foreach_point_neighbor, foreach_neighbor +export foreach_point_neighbor, foreach_neighbor, PeriodicBox export TrivialNeighborhoodSearch, GridNeighborhoodSearch, PrecomputedNeighborhoodSearch export initialize!, update!, initialize_grid!, update_grid! diff --git a/src/neighborhood_search.jl b/src/neighborhood_search.jl index 87ab1296..0dc4f9e4 100644 --- a/src/neighborhood_search.jl +++ b/src/neighborhood_search.jl @@ -41,14 +41,26 @@ See also [`initialize!`](@ref). return search end +""" + PeriodicBox(; min_corner, max_corner) + +Define a rectangular periodic domain. + +# Keywords +- `min_corner`: Coordinates of the domain corner in negative coordinate directions. +- `max_corner`: Coordinates of the domain corner in positive coordinate directions. +""" struct PeriodicBox{NDIMS, ELTYPE} min_corner :: SVector{NDIMS, ELTYPE} max_corner :: SVector{NDIMS, ELTYPE} size :: SVector{NDIMS, ELTYPE} - function PeriodicBox(min_corner, max_corner) - new{length(min_corner), eltype(min_corner)}(min_corner, max_corner, - max_corner - min_corner) + function PeriodicBox(; min_corner, max_corner) + min_corner_ = SVector(Tuple(min_corner)) + max_corner_ = SVector(Tuple(max_corner)) + + new{length(min_corner), eltype(min_corner)}(min_corner_, max_corner_, + max_corner_ - min_corner_) end end diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index af4b690d..6058af26 100644 --- a/src/nhs_grid.jl +++ b/src/nhs_grid.jl @@ -1,6 +1,6 @@ @doc raw""" - GridNeighborhoodSearch{NDIMS}(search_radius, n_points; periodic_box_min_corner=nothing, - periodic_box_max_corner=nothing, threaded_nhs_update=true) + GridNeighborhoodSearch{NDIMS}(search_radius, n_points; + periodic_box = nothing, threaded_nhs_update = true) Simple grid-based neighborhood search with uniform search radius. The domain is divided into a regular grid. @@ -28,15 +28,12 @@ since not sorting makes our implementation a lot faster (although less paralleli - `n_points`: Total number of points. # Keywords -- `periodic_box_min_corner`: In order to use a (rectangular) periodic domain, pass the - coordinates of the domain corner in negative coordinate - directions. -- `periodic_box_max_corner`: In order to use a (rectangular) periodic domain, pass the - coordinates of the domain corner in positive coordinate - directions. -- `threaded_nhs_update=true`: Can be used to deactivate thread parallelization in the neighborhood search update. - This can be one of the largest sources of variations between simulations - with different thread numbers due to neighbor ordering changes. +- `periodic_box = nothing`: In order to use a (rectangular) periodic domain, pass a + [`PeriodicBox`](@ref). +- `threaded_nhs_update = true`: Can be used to deactivate thread parallelization in the + neighborhood search update. This can be one of the largest + sources of variations between simulations with different + thread numbers due to neighbor ordering changes. ## References - M. Chalela, E. Sillero, L. Pereyra, M.A. Garcia, J.B. Cabral, M. Lares, M. Merchán. @@ -59,8 +56,7 @@ struct GridNeighborhoodSearch{NDIMS, ELTYPE, CL, PB} <: AbstractNeighborhoodSear threaded_nhs_update :: Bool function GridNeighborhoodSearch{NDIMS}(search_radius, n_points; - periodic_box_min_corner = nothing, - periodic_box_max_corner = nothing, + periodic_box = nothing, threaded_nhs_update = true) where {NDIMS} ELTYPE = typeof(search_radius) cell_list = DictionaryCellList{NDIMS}() @@ -68,16 +64,11 @@ struct GridNeighborhoodSearch{NDIMS, ELTYPE, CL, PB} <: AbstractNeighborhoodSear cell_buffer = Array{NTuple{NDIMS, Int}, 2}(undef, n_points, Threads.nthreads()) cell_buffer_indices = zeros(Int, Threads.nthreads()) - if search_radius < eps() || - (periodic_box_min_corner === nothing && periodic_box_max_corner === nothing) - + if search_radius < eps() || isnothing(periodic_box) # No periodicity - periodic_box = nothing n_cells = ntuple(_ -> -1, Val(NDIMS)) cell_size = ntuple(_ -> search_radius, Val(NDIMS)) - elseif periodic_box_min_corner !== nothing && periodic_box_max_corner !== nothing - periodic_box = PeriodicBox(periodic_box_min_corner, periodic_box_max_corner) - + else # Round up search radius so that the grid fits exactly into the domain without # splitting any cells. This might impact performance slightly, since larger # cells mean that more potential neighbors are considered than necessary. @@ -91,9 +82,6 @@ struct GridNeighborhoodSearch{NDIMS, ELTYPE, CL, PB} <: AbstractNeighborhoodSear "in each dimension when used with periodicity. " * "Please use no NHS for very small problems.")) end - else - throw(ArgumentError("`periodic_box_min_corner` and `periodic_box_max_corner` " * - "must either be both `nothing` or both an array or tuple")) end new{NDIMS, ELTYPE, typeof(cell_list), @@ -103,7 +91,7 @@ struct GridNeighborhoodSearch{NDIMS, ELTYPE, CL, PB} <: AbstractNeighborhoodSear end end -@inline Base.ndims(neighborhood_search::GridNeighborhoodSearch{NDIMS}) where {NDIMS} = NDIMS +@inline Base.ndims(::GridNeighborhoodSearch{NDIMS}) where {NDIMS} = NDIMS @inline function npoints(neighborhood_search::GridNeighborhoodSearch) return size(neighborhood_search.cell_buffer, 1) @@ -321,13 +309,8 @@ end # Create a copy of a neighborhood search but with a different search radius function copy_neighborhood_search(nhs::GridNeighborhoodSearch, search_radius, x, y) - if nhs.periodic_box === nothing - search = GridNeighborhoodSearch{ndims(nhs)}(search_radius, npoints(nhs)) - else - search = GridNeighborhoodSearch{ndims(nhs)}(search_radius, npoints(nhs), - periodic_box_min_corner = nhs.periodic_box.min_corner, - periodic_box_max_corner = nhs.periodic_box.max_corner) - end + search = GridNeighborhoodSearch{ndims(nhs)}(search_radius, npoints(nhs), + periodic_box = nhs.periodic_box) # Initialize neighborhood search initialize!(search, x, y) diff --git a/src/nhs_neighbor_lists.jl b/src/nhs_neighbor_lists.jl index 07854b65..21c9e8a9 100644 --- a/src/nhs_neighbor_lists.jl +++ b/src/nhs_neighbor_lists.jl @@ -1,7 +1,6 @@ @doc raw""" PrecomputedNeighborhoodSearch{NDIMS}(search_radius, n_points; - periodic_box_min_corner = nothing, - periodic_box_max_corner = nothing) + periodic_box = nothing) Neighborhood search with precomputed neighbor lists. A list of all neighbors is computed for each point during initialization and update. @@ -17,12 +16,8 @@ initialization and update. - `n_points`: Total number of points. # Keywords -- `periodic_box_min_corner`: In order to use a (rectangular) periodic domain, pass the - coordinates of the domain corner in negative coordinate - directions. -- `periodic_box_max_corner`: In order to use a (rectangular) periodic domain, pass the - coordinates of the domain corner in positive coordinate - directions. +- `periodic_box = nothing`: In order to use a (rectangular) periodic domain, pass a + [`PeriodicBox`](@ref). """ struct PrecomputedNeighborhoodSearch{NDIMS, NHS, NL, PB} neighborhood_search :: NHS @@ -30,19 +25,15 @@ struct PrecomputedNeighborhoodSearch{NDIMS, NHS, NL, PB} periodic_box :: PB function PrecomputedNeighborhoodSearch{NDIMS}(search_radius, n_points; - periodic_box_min_corner = nothing, - periodic_box_max_corner = nothing) where { - NDIMS - } + periodic_box = nothing) where {NDIMS} nhs = GridNeighborhoodSearch{NDIMS}(search_radius, n_points, - periodic_box_min_corner = periodic_box_min_corner, - periodic_box_max_corner = periodic_box_max_corner) + periodic_box = periodic_box) neighbor_lists = Vector{Vector{Int}}() new{NDIMS, typeof(nhs), typeof(neighbor_lists), - typeof(nhs.periodic_box)}(nhs, neighbor_lists, nhs.periodic_box) + typeof(periodic_box)}(nhs, neighbor_lists, periodic_box) end end diff --git a/src/nhs_trivial.jl b/src/nhs_trivial.jl index d99091e3..eeaff5b0 100644 --- a/src/nhs_trivial.jl +++ b/src/nhs_trivial.jl @@ -1,20 +1,16 @@ @doc raw""" - TrivialNeighborhoodSearch{NDIMS}(search_radius, eachpoint) + TrivialNeighborhoodSearch{NDIMS}(search_radius, eachpoint, periodic_box = nothing) Trivial neighborhood search that simply loops over all points. # Arguments - `NDIMS`: Number of dimensions. - `search_radius`: The uniform search radius. -- `eachpoint`: `UnitRange` of all point indices. Usually just `1:n_points`. +- `eachparticle`: Iterator for all point indices. Usually just `1:n_points`. # Keywords -- `periodic_box_min_corner`: In order to use a (rectangular) periodic domain, pass the - coordinates of the domain corner in negative coordinate - directions. -- `periodic_box_max_corner`: In order to use a (rectangular) periodic domain, pass the - coordinates of the domain corner in positive coordinate - directions. +- `periodic_box = nothing`: In order to use a (rectangular) periodic domain, pass a + [`PeriodicBox`](@ref). """ struct TrivialNeighborhoodSearch{NDIMS, ELTYPE, EP, PB} <: AbstractNeighborhoodSearch search_radius :: ELTYPE @@ -22,32 +18,13 @@ struct TrivialNeighborhoodSearch{NDIMS, ELTYPE, EP, PB} <: AbstractNeighborhoodS periodic_box :: PB function TrivialNeighborhoodSearch{NDIMS}(search_radius, eachpoint; - periodic_box_min_corner = nothing, - periodic_box_max_corner = nothing) where { - NDIMS - } - if search_radius < eps() || - (periodic_box_min_corner === nothing && periodic_box_max_corner === nothing) - - # No periodicity - periodic_box = nothing - elseif periodic_box_min_corner !== nothing && periodic_box_max_corner !== nothing - periodic_box = PeriodicBox(periodic_box_min_corner, periodic_box_max_corner) - else - throw(ArgumentError("`periodic_box_min_corner` and `periodic_box_max_corner` " * - "must either be both `nothing` or both an array or tuple")) - end - + periodic_box = nothing) where {NDIMS} new{NDIMS, typeof(search_radius), typeof(eachpoint), typeof(periodic_box)}(search_radius, eachpoint, periodic_box) end end -@inline function Base.ndims(neighborhood_search::TrivialNeighborhoodSearch{NDIMS}) where { - NDIMS - } - return NDIMS -end +@inline Base.ndims(::TrivialNeighborhoodSearch{NDIMS}) where {NDIMS} = NDIMS @inline initialize!(search::TrivialNeighborhoodSearch, x, y) = search diff --git a/test/neighborhood_search.jl b/test/neighborhood_search.jl index 2027328b..3111c8ba 100644 --- a/test/neighborhood_search.jl +++ b/test/neighborhood_search.jl @@ -24,11 +24,11 @@ ] periodic_boxes = [ - ([-0.1, -0.2], [0.2, 0.4]), + PeriodicBox(min_corner = [-0.1, -0.2], max_corner = [0.2, 0.4]), # The `GridNeighborhoodSearch` is forced to round up the cell sizes in this test # to avoid split cells. - ([-0.1, -0.2], [0.205, 0.43]), - ([-0.1, -0.2, 0.05], [0.2, 0.4, 0.35]), + PeriodicBox(min_corner = [-0.1, -0.2], max_corner = [0.205, 0.43]), + PeriodicBox(min_corner = [-0.1, -0.2, 0.05], max_corner = [0.2, 0.4, 0.35]), ] @testset verbose=true "$(names[i])" for i in eachindex(names) @@ -40,14 +40,11 @@ neighborhood_searches = [ TrivialNeighborhoodSearch{NDIMS}(search_radius, 1:n_points, - periodic_box_min_corner = periodic_boxes[i][1], - periodic_box_max_corner = periodic_boxes[i][2]), + periodic_box = periodic_boxes[i]), GridNeighborhoodSearch{NDIMS}(search_radius, n_points, - periodic_box_min_corner = periodic_boxes[i][1], - periodic_box_max_corner = periodic_boxes[i][2]), + periodic_box = periodic_boxes[i]), PrecomputedNeighborhoodSearch{NDIMS}(search_radius, n_points, - periodic_box_min_corner = periodic_boxes[i][1], - periodic_box_max_corner = periodic_boxes[i][2]), + periodic_box = periodic_boxes[i]), ] neighborhood_searches_names = [ "`TrivialNeighborhoodSearch`", diff --git a/test/nhs_grid.jl b/test/nhs_grid.jl index 68aca46c..c6a59b5c 100644 --- a/test/nhs_grid.jl +++ b/test/nhs_grid.jl @@ -148,19 +148,18 @@ ] periodic_boxes = [ - ([-0.1, -0.2], [0.2, 0.4]), + PeriodicBox(min_corner = [-0.1, -0.2], max_corner = [0.2, 0.4]), # The `GridNeighborhoodSearch` is forced to round up the cell sizes in this test # to avoid split cells. - ([-0.1, -0.2], [0.205, 0.43]), - ([-0.1, -0.2, 0.05], [0.2, 0.4, 0.35]), + PeriodicBox(min_corner = [-0.1, -0.2], max_corner = [0.205, 0.43]), + PeriodicBox(min_corner = [-0.1, -0.2, 0.05], max_corner = [0.2, 0.4, 0.35]), ] @testset verbose=true "$(names[i])" for i in eachindex(names) coords = coordinates[i] nhs = GridNeighborhoodSearch{size(coords, 1)}(0.1, size(coords, 2), - periodic_box_min_corner = periodic_boxes[i][1], - periodic_box_max_corner = periodic_boxes[i][2]) + periodic_box = periodic_boxes[i]) initialize_grid!(nhs, coords) @@ -189,13 +188,18 @@ # 5 x 1 cells nhs = GridNeighborhoodSearch{2}(1.0, size(coords, 2), - periodic_box_min_corner = [-1.5, 0.0], - periodic_box_max_corner = [2.5, 3.0]) + periodic_box = PeriodicBox(min_corner = [ + -1.5, + 0.0, + ], + max_corner = [ + 2.5, + 3.0, + ])) initialize_grid!(nhs, coords) - neighbors = [sort(unique(collect(PointNeighbors.eachneighbor(coords[:, - i], + neighbors = [sort(unique(collect(PointNeighbors.eachneighbor(coords[:, i], nhs)))) for i in 1:2] From 5d52dfd7ed6363801cf1763aa506fd62ac71de0f Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Mon, 24 Jun 2024 10:35:47 +0200 Subject: [PATCH 4/4] Implement `copy_neighborhood_search` to allow empty "template" neighborhood search (#32) * Rename `threaded_nhs_update` to `threaded_update` * Add `threaded_update` option to `PrecomputedNeighborhoodSearch` * Properly implement `copy_neighborhood_search` * Create empty NHS templates by default * Fix tests * Add tests for copied templates * Add doctest * Fix rebase * Break doctests to verify that they are running * Revert "Break doctests to verify that they are running" This reverts commit 7e4275dc962eaea627a380185d3b6ce892b3bd9a. * Fix typo * Implement suggestions --- src/PointNeighbors.jl | 5 +- src/neighborhood_search.jl | 31 ++++++++++++ src/nhs_grid.jl | 50 +++++++++---------- ...s_neighbor_lists.jl => nhs_precomputed.jl} | 34 +++++++++---- src/nhs_trivial.jl | 20 ++++++-- test/neighborhood_search.jl | 45 +++++++++++++---- test/nhs_grid.jl | 20 ++++---- test/nhs_trivial.jl | 2 +- test/runtests.jl | 2 +- 9 files changed, 146 insertions(+), 63 deletions(-) rename src/{nhs_neighbor_lists.jl => nhs_precomputed.jl} (67%) diff --git a/src/PointNeighbors.jl b/src/PointNeighbors.jl index 68c0e41c..4e294c3a 100644 --- a/src/PointNeighbors.jl +++ b/src/PointNeighbors.jl @@ -11,10 +11,11 @@ include("neighborhood_search.jl") include("nhs_trivial.jl") include("cell_lists/cell_lists.jl") include("nhs_grid.jl") -include("nhs_neighbor_lists.jl") +include("nhs_precomputed.jl") -export foreach_point_neighbor, foreach_neighbor, PeriodicBox +export foreach_point_neighbor, foreach_neighbor export TrivialNeighborhoodSearch, GridNeighborhoodSearch, PrecomputedNeighborhoodSearch export initialize!, update!, initialize_grid!, update_grid! +export PeriodicBox, copy_neighborhood_search end # module PointNeighbors diff --git a/src/neighborhood_search.jl b/src/neighborhood_search.jl index 0dc4f9e4..c44ac1b0 100644 --- a/src/neighborhood_search.jl +++ b/src/neighborhood_search.jl @@ -41,6 +41,37 @@ See also [`initialize!`](@ref). return search end +""" + copy_neighborhood_search(search::AbstractNeighborhoodSearch, search_radius, n_points; + eachpoint = 1:n_points) + +Create a new **uninitialized** neighborhood search of the same type and with the same +configuration options as `search`, but with a different search radius and number of points. + +The [`TrivialNeighborhoodSearch`](@ref) also requires an iterator `eachpoint`, which most +of the time will be `1:n_points`. If the `TrivialNeighborhoodSearch` is never going to be +used, the keyword argument `eachpoint` can be ignored. + +This is useful when a simulation code requires multiple neighborhood searches of the same +kind. One can then just pass an empty neighborhood search as a template and use +this function inside the simulation code to generate similar neighborhood searches with +different search radii and different numbers of points. +```jldoctest; filter = r"GridNeighborhoodSearch{2,.*" +# Template +nhs = GridNeighborhoodSearch{2}() + +# Inside the simulation code, generate similar neighborhood searches +nhs1 = copy_neighborhood_search(nhs, 1.0, 100) + +# output +GridNeighborhoodSearch{2, Float64, ...}(...) +``` +""" +@inline function copy_neighborhood_search(search::AbstractNeighborhoodSearch, + search_radius, n_points; eachpoint = 1:n_points) + return nothing +end + """ PeriodicBox(; min_corner, max_corner) diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index 6058af26..1dbc0fd3 100644 --- a/src/nhs_grid.jl +++ b/src/nhs_grid.jl @@ -1,6 +1,6 @@ @doc raw""" - GridNeighborhoodSearch{NDIMS}(search_radius, n_points; - periodic_box = nothing, threaded_nhs_update = true) + GridNeighborhoodSearch{NDIMS}(; search_radius = 0.0, n_points = 0, + periodic_box = nothing, threaded_update = true) Simple grid-based neighborhood search with uniform search radius. The domain is divided into a regular grid. @@ -23,17 +23,19 @@ As opposed to (Ihmsen et al. 2011), we do not sort the points in any way, since not sorting makes our implementation a lot faster (although less parallelizable). # Arguments -- `NDIMS`: Number of dimensions. -- `search_radius`: The uniform search radius. -- `n_points`: Total number of points. +- `NDIMS`: Number of dimensions. # Keywords +- `search_radius = 0.0`: The fixed search radius. The default of `0.0` is useful together + with [`copy_neighborhood_search`](@ref). +- `n_points = 0`: Total number of points. The default of `0` is useful together + with [`copy_neighborhood_search`](@ref). - `periodic_box = nothing`: In order to use a (rectangular) periodic domain, pass a [`PeriodicBox`](@ref). -- `threaded_nhs_update = true`: Can be used to deactivate thread parallelization in the - neighborhood search update. This can be one of the largest - sources of variations between simulations with different - thread numbers due to neighbor ordering changes. +- `threaded_update = true`: Can be used to deactivate thread parallelization in the + neighborhood search update. This can be one of the largest + sources of variations between simulations with different + thread numbers due to neighbor ordering changes. ## References - M. Chalela, E. Sillero, L. Pereyra, M.A. Garcia, J.B. Cabral, M. Lares, M. Merchán. @@ -53,11 +55,11 @@ struct GridNeighborhoodSearch{NDIMS, ELTYPE, CL, PB} <: AbstractNeighborhoodSear cell_size :: NTuple{NDIMS, ELTYPE} # Required to calculate cell index cell_buffer :: Array{NTuple{NDIMS, Int}, 2} # Multithreaded buffer for `update!` cell_buffer_indices :: Vector{Int} # Store which entries of `cell_buffer` are initialized - threaded_nhs_update :: Bool + threaded_update :: Bool - function GridNeighborhoodSearch{NDIMS}(search_radius, n_points; + function GridNeighborhoodSearch{NDIMS}(; search_radius = 0.0, n_points = 0, periodic_box = nothing, - threaded_nhs_update = true) where {NDIMS} + threaded_update = true) where {NDIMS} ELTYPE = typeof(search_radius) cell_list = DictionaryCellList{NDIMS}() @@ -87,7 +89,7 @@ struct GridNeighborhoodSearch{NDIMS, ELTYPE, CL, PB} <: AbstractNeighborhoodSear new{NDIMS, ELTYPE, typeof(cell_list), typeof(periodic_box)}(cell_list, search_radius, periodic_box, n_cells, cell_size, cell_buffer, cell_buffer_indices, - threaded_nhs_update) + threaded_update) end end @@ -141,14 +143,14 @@ end # Modify the existing hash table by moving points into their new cells function update_grid!(neighborhood_search::GridNeighborhoodSearch, coords_fun) - (; cell_list, cell_buffer, cell_buffer_indices, threaded_nhs_update) = neighborhood_search + (; cell_list, cell_buffer, cell_buffer_indices, threaded_update) = neighborhood_search # Reset `cell_buffer` by moving all pointers to the beginning cell_buffer_indices .= 0 # Find all cells containing points that now belong to another cell mark_changed_cell!(neighborhood_search, cell_list, coords_fun, - Val(threaded_nhs_update)) + Val(threaded_update)) # Iterate over all marked cells and move the points into their new cells. for thread in 1:Threads.nthreads() @@ -180,7 +182,7 @@ function update_grid!(neighborhood_search::GridNeighborhoodSearch, coords_fun) end @inline function mark_changed_cell!(neighborhood_search, cell_list, coords_fun, - threaded_nhs_update::Val{true}) + threaded_update::Val{true}) # `collect` the keyset to be able to loop over it with `@threaded` @threaded for cell in collect(eachcell(cell_list)) mark_changed_cell!(neighborhood_search, cell, coords_fun) @@ -188,7 +190,7 @@ end end @inline function mark_changed_cell!(neighborhood_search, cell_list, coords_fun, - threaded_nhs_update::Val{false}) + threaded_update::Val{false}) for cell in eachcell(cell_list) mark_changed_cell!(neighborhood_search, cell, coords_fun) end @@ -307,13 +309,9 @@ end return floor(Int, i) end -# Create a copy of a neighborhood search but with a different search radius -function copy_neighborhood_search(nhs::GridNeighborhoodSearch, search_radius, x, y) - search = GridNeighborhoodSearch{ndims(nhs)}(search_radius, npoints(nhs), - periodic_box = nhs.periodic_box) - - # Initialize neighborhood search - initialize!(search, x, y) - - return search +function copy_neighborhood_search(nhs::GridNeighborhoodSearch, search_radius, n_points; + eachpoint = 1:n_points) + return GridNeighborhoodSearch{ndims(nhs)}(; search_radius, n_points, + periodic_box = nhs.periodic_box, + threaded_update = nhs.threaded_update) end diff --git a/src/nhs_neighbor_lists.jl b/src/nhs_precomputed.jl similarity index 67% rename from src/nhs_neighbor_lists.jl rename to src/nhs_precomputed.jl index 21c9e8a9..ddf8ac1c 100644 --- a/src/nhs_neighbor_lists.jl +++ b/src/nhs_precomputed.jl @@ -1,6 +1,6 @@ @doc raw""" - PrecomputedNeighborhoodSearch{NDIMS}(search_radius, n_points; - periodic_box = nothing) + PrecomputedNeighborhoodSearch{NDIMS}(; search_radius = 0.0, n_points = 0, + periodic_box = nothing, threaded_update = true) Neighborhood search with precomputed neighbor lists. A list of all neighbors is computed for each point during initialization and update. @@ -11,23 +11,31 @@ A [`GridNeighborhoodSearch`](@ref) is used internally to compute the neighbor li initialization and update. # Arguments -- `NDIMS`: Number of dimensions. -- `search_radius`: The uniform search radius. -- `n_points`: Total number of points. +- `NDIMS`: Number of dimensions. # Keywords +- `search_radius = 0.0`: The fixed search radius. The default of `0.0` is useful together + with [`copy_neighborhood_search`](@ref). +- `n_points = 0`: Total number of points. The default of `0` is useful together + with [`copy_neighborhood_search`](@ref). - `periodic_box = nothing`: In order to use a (rectangular) periodic domain, pass a [`PeriodicBox`](@ref). +- `threaded_update = true`: Can be used to deactivate thread parallelization in the + neighborhood search update. This can be one of the largest + sources of variations between simulations with different + thread numbers due to neighbor ordering changes. """ struct PrecomputedNeighborhoodSearch{NDIMS, NHS, NL, PB} neighborhood_search :: NHS neighbor_lists :: NL periodic_box :: PB - function PrecomputedNeighborhoodSearch{NDIMS}(search_radius, n_points; - periodic_box = nothing) where {NDIMS} - nhs = GridNeighborhoodSearch{NDIMS}(search_radius, n_points, - periodic_box = periodic_box) + function PrecomputedNeighborhoodSearch{NDIMS}(; search_radius = 0.0, n_points = 0, + periodic_box = nothing, + threaded_update = true) where {NDIMS} + nhs = GridNeighborhoodSearch{NDIMS}(; search_radius, n_points, + periodic_box = periodic_box, + threaded_update = threaded_update) neighbor_lists = Vector{Vector{Int}}() @@ -103,3 +111,11 @@ end @inline f(point, neighbor, pos_diff, distance) end end + +function copy_neighborhood_search(nhs::PrecomputedNeighborhoodSearch, + search_radius, n_points; eachpoint = 1:n_points) + threaded_update = nhs.neighborhood_search.threaded_update + return PrecomputedNeighborhoodSearch{ndims(nhs)}(; search_radius, n_points, + periodic_box = nhs.periodic_box, + threaded_update) +end diff --git a/src/nhs_trivial.jl b/src/nhs_trivial.jl index eeaff5b0..13b24d58 100644 --- a/src/nhs_trivial.jl +++ b/src/nhs_trivial.jl @@ -1,14 +1,18 @@ @doc raw""" - TrivialNeighborhoodSearch{NDIMS}(search_radius, eachpoint, periodic_box = nothing) + TrivialNeighborhoodSearch{NDIMS}(; search_radius = 0.0, eachpoint = 1:0, + periodic_box = nothing) Trivial neighborhood search that simply loops over all points. # Arguments -- `NDIMS`: Number of dimensions. -- `search_radius`: The uniform search radius. -- `eachparticle`: Iterator for all point indices. Usually just `1:n_points`. +- `NDIMS`: Number of dimensions. # Keywords +- `search_radius = 0.0`: The fixed search radius. The default of `0.0` is useful together + with [`copy_neighborhood_search`](@ref). +- `eachpoint = 1:0`: Iterator for all point indices. Usually just `1:n_points`. + The default of `1:0` is useful together with + [`copy_neighborhood_search`](@ref). - `periodic_box = nothing`: In order to use a (rectangular) periodic domain, pass a [`PeriodicBox`](@ref). """ @@ -17,7 +21,7 @@ struct TrivialNeighborhoodSearch{NDIMS, ELTYPE, EP, PB} <: AbstractNeighborhoodS eachpoint :: EP periodic_box :: PB - function TrivialNeighborhoodSearch{NDIMS}(search_radius, eachpoint; + function TrivialNeighborhoodSearch{NDIMS}(; search_radius = 0.0, eachpoint = 1:0, periodic_box = nothing) where {NDIMS} new{NDIMS, typeof(search_radius), typeof(eachpoint), typeof(periodic_box)}(search_radius, eachpoint, periodic_box) @@ -39,3 +43,9 @@ end function copy_neighborhood_search(nhs::TrivialNeighborhoodSearch, search_radius, x, y) return nhs end + +function copy_neighborhood_search(nhs::TrivialNeighborhoodSearch, + search_radius, n_points; eachpoint = 1:n_points) + return TrivialNeighborhoodSearch{ndims(nhs)}(; search_radius, eachpoint, + periodic_box = nhs.periodic_box) +end diff --git a/test/neighborhood_search.jl b/test/neighborhood_search.jl index 3111c8ba..7705257d 100644 --- a/test/neighborhood_search.jl +++ b/test/neighborhood_search.jl @@ -39,21 +39,34 @@ search_radius = 0.1 neighborhood_searches = [ - TrivialNeighborhoodSearch{NDIMS}(search_radius, 1:n_points, + TrivialNeighborhoodSearch{NDIMS}(; search_radius, eachpoint = 1:n_points, periodic_box = periodic_boxes[i]), - GridNeighborhoodSearch{NDIMS}(search_radius, n_points, + GridNeighborhoodSearch{NDIMS}(; search_radius, n_points, periodic_box = periodic_boxes[i]), - PrecomputedNeighborhoodSearch{NDIMS}(search_radius, n_points, + PrecomputedNeighborhoodSearch{NDIMS}(; search_radius, n_points, periodic_box = periodic_boxes[i]), ] - neighborhood_searches_names = [ + + names = [ "`TrivialNeighborhoodSearch`", "`GridNeighborhoodSearch`", "`PrecomputedNeighborhoodSearch`", ] + # Also test copied templates + template_nhs = [ + TrivialNeighborhoodSearch{NDIMS}(periodic_box = periodic_boxes[i]), + GridNeighborhoodSearch{NDIMS}(periodic_box = periodic_boxes[i]), + PrecomputedNeighborhoodSearch{NDIMS}(periodic_box = periodic_boxes[i]), + ] + copied_nhs = copy_neighborhood_search.(template_nhs, search_radius, n_points) + append!(neighborhood_searches, copied_nhs) + + names_copied = [name * " copied" for name in names] + append!(names, names_copied) + # Run this for every neighborhood search - @testset "$(neighborhood_searches_names[j])" for j in eachindex(neighborhood_searches_names) + @testset "$(names[j])" for j in eachindex(names) nhs = neighborhood_searches[j] initialize!(nhs, coords, coords) @@ -103,7 +116,8 @@ # Compute expected neighbor lists by brute-force looping over all points # as potential neighbors (`TrivialNeighborhoodSearch`). - trivial_nhs = TrivialNeighborhoodSearch{NDIMS}(search_radius, axes(coords, 2)) + trivial_nhs = TrivialNeighborhoodSearch{NDIMS}(; search_radius, + eachpoint = axes(coords, 2)) neighbors_expected = [Int[] for _ in axes(coords, 2)] @@ -114,16 +128,27 @@ end neighborhood_searches = [ - GridNeighborhoodSearch{NDIMS}(search_radius, n_points), - PrecomputedNeighborhoodSearch{NDIMS}(search_radius, n_points), + GridNeighborhoodSearch{NDIMS}(; search_radius, n_points), + PrecomputedNeighborhoodSearch{NDIMS}(; search_radius, n_points), ] - neighborhood_searches_names = [ + names = [ "`GridNeighborhoodSearch`", "`PrecomputedNeighborhoodSearch`", ] - @testset "$(neighborhood_searches_names[i])" for i in eachindex(neighborhood_searches_names) + # Also test copied templates + template_nhs = [ + GridNeighborhoodSearch{NDIMS}(), + PrecomputedNeighborhoodSearch{NDIMS}(), + ] + copied_nhs = copy_neighborhood_search.(template_nhs, search_radius, n_points) + append!(neighborhood_searches, copied_nhs) + + names_copied = [name * " copied" for name in names] + append!(names, names_copied) + + @testset "$(names[i])" for i in eachindex(names) nhs = neighborhood_searches[i] # Initialize with `seed = 1` diff --git a/test/nhs_grid.jl b/test/nhs_grid.jl index c6a59b5c..d7451c4c 100644 --- a/test/nhs_grid.jl +++ b/test/nhs_grid.jl @@ -19,13 +19,13 @@ # from (x, y) = (-0.25, -0.25) to (x, y) = (0.35, 0.35). range = -0.25:0.1:0.35 coordinates1 = hcat(collect.(Iterators.product(range, range))...) - npoints = size(coordinates1, 2) + n_points = size(coordinates1, 2) point_position1 = [0.05, 0.05] - radius = 0.1 + search_radius = 0.1 # Create neighborhood search - nhs1 = GridNeighborhoodSearch{2}(radius, npoints) + nhs1 = GridNeighborhoodSearch{2}(; search_radius, n_points) initialize_grid!(nhs1, coordinates1) @@ -48,7 +48,8 @@ neighbors3 = sort(collect(PointNeighbors.eachneighbor(point_position2, nhs1))) # Double search radius - nhs2 = GridNeighborhoodSearch{2}(2 * radius, size(coordinates1, 2)) + nhs2 = GridNeighborhoodSearch{2}(search_radius = 2 * search_radius, + n_points = size(coordinates1, 2)) initialize!(nhs2, coordinates1, coordinates1) # Get each neighbor in double search radius @@ -83,13 +84,13 @@ # from (x, y, z) = (-0.25, -0.25, -0.25) to (x, y, z) = (0.35, 0.35, 0.35). range = -0.25:0.1:0.35 coordinates1 = hcat(collect.(Iterators.product(range, range, range))...) - npoints = size(coordinates1, 2) + n_points = size(coordinates1, 2) point_position1 = [0.05, 0.05, 0.05] - radius = 0.1 + search_radius = 0.1 # Create neighborhood search - nhs1 = GridNeighborhoodSearch{3}(radius, npoints) + nhs1 = GridNeighborhoodSearch{3}(; search_radius, n_points) coords_fun(i) = coordinates1[:, i] initialize_grid!(nhs1, coords_fun) @@ -158,7 +159,8 @@ @testset verbose=true "$(names[i])" for i in eachindex(names) coords = coordinates[i] - nhs = GridNeighborhoodSearch{size(coords, 1)}(0.1, size(coords, 2), + nhs = GridNeighborhoodSearch{size(coords, 1)}(search_radius = 0.1, + n_points = size(coords, 2), periodic_box = periodic_boxes[i]) initialize_grid!(nhs, coords) @@ -187,7 +189,7 @@ 0.0 0.0] # 5 x 1 cells - nhs = GridNeighborhoodSearch{2}(1.0, size(coords, 2), + nhs = GridNeighborhoodSearch{2}(search_radius = 1.0, n_points = size(coords, 2), periodic_box = PeriodicBox(min_corner = [ -1.5, 0.0, diff --git a/test/nhs_trivial.jl b/test/nhs_trivial.jl index 096dbf0c..441a1578 100644 --- a/test/nhs_trivial.jl +++ b/test/nhs_trivial.jl @@ -1,6 +1,6 @@ @testset verbose=true "TrivialNeighborhoodSearch" begin # Setup with 5 points - nhs = TrivialNeighborhoodSearch{2}(1.0, Base.OneTo(5)) + nhs = TrivialNeighborhoodSearch{2}(search_radius = 1.0, eachpoint = Base.OneTo(5)) # Get each neighbor for arbitrary coordinates neighbors = collect(PointNeighbors.eachneighbor([1.0, 2.0], nhs)) diff --git a/test/runtests.jl b/test/runtests.jl index ad5c7bab..0eb0f4b5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,4 +4,4 @@ include("test_util.jl") include("nhs_trivial.jl") include("nhs_grid.jl") include("neighborhood_search.jl") -end +end;