Skip to content
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.

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

Replace all references to "particle" by "point" #31

Merged
merged 3 commits into from
Jun 21, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
2 changes: 1 addition & 1 deletion src/PointNeighbors.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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!

Expand Down
6 changes: 3 additions & 3 deletions src/cell_lists/dictionary.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
50 changes: 24 additions & 26 deletions src/neighborhood_search.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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`.

Expand All @@ -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

Expand All @@ -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,
Expand All @@ -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
Expand All @@ -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)
Expand Down
92 changes: 46 additions & 46 deletions src/nhs_grid.jl
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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
Expand All @@ -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.
Expand All @@ -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() ||
Expand Down Expand Up @@ -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

Expand All @@ -124,23 +124,23 @@ 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
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
Expand All @@ -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

Expand Down Expand Up @@ -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,
Expand All @@ -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

Expand All @@ -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

Expand All @@ -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)]
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down
Loading