PointNeighbors.jl API

DictionaryCellList{NDIMS}()

A simple cell list implementation where a cell index (i, j) or (i, j, k) is mapped to a Vector{Int} by a Dict. By using a dictionary, which only stores non-empty cells, the domain is potentially infinite.

This implementation is very simple, but it neither uses an optimized hash function for integer tuples, nor does it use a contiguous memory layout. Consequently, this cell list is not GPU-compatible.

Arguments

• NDIMS: Number of dimensions.

FullGridCellList(; min_corner, max_corner, search_radius = 0.0,
                 periodicity = false, backend = DynamicVectorOfVectors{Int32},
                 max_points_per_cell = 100)

A simple cell list implementation where each (empty or non-empty) cell of a rectangular (axis-aligned) domain is assigned a list of points. This cell list only works when all points are inside the specified domain at all times.

Only set min_corner and max_corner and use the default values for the other arguments to create an empty "template" cell list that can be used to create an empty "template" neighborhood search. See copy_neighborhood_search for more details.

Keywords

• min_corner: Coordinates of the domain corner in negative coordinate directions.
• max_corner: Coordinates of the domain corner in positive coordinate directions.
• search_radius = 0.0: Search radius of the neighborhood search, which will determine the cell size. Use the default of 0.0 to create a template (see above).
• periodicity = false: Set to true when using a PeriodicBox with the neighborhood search. When using copy_neighborhood_search, this option can be ignored an will be set automatically depending on the periodicity of the neighborhood search.
• backend = DynamicVectorOfVectors{Int32}: Type of the data structure to store the actual cell lists. Can be
  • Vector{Vector{Int32}}: Scattered memory, but very memory-efficient.
  • DynamicVectorOfVectors{Int32}: Contiguous memory, optimizing cache-hits.
• max_points_per_cell = 100: Maximum number of points per cell. This will be used to allocate the DynamicVectorOfVectors. It is not used with the Vector{Vector{Int32}} backend.

GridNeighborhoodSearch{NDIMS}(; search_radius = 0.0, n_points = 0,
                              periodic_box = nothing,
                              cell_list = DictionaryCellList{NDIMS}(),
                              update_strategy = nothing)

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

\[\left( \left\lfloor \frac{x}{d} \right\rfloor, \left\lfloor \frac{y}{d} \right\rfloor \right) \quad \text{or} \quad \left( \left\lfloor \frac{x}{d} \right\rfloor, \left\lfloor \frac{y}{d} \right\rfloor, \left\lfloor \frac{z}{d} \right\rfloor \right),\]

where $x, y, z$ are the space coordinates and $d$ is the search radius.

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 points in any way, since not sorting makes our implementation a lot faster (although less parallelizable).

Arguments

• 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.
• n_points = 0: Total number of points. The default of 0 is useful together with copy_neighborhood_search.
• periodic_box = nothing: In order to use a (rectangular) periodic domain, pass a PeriodicBox.
• cell_list: The cell list that maps a cell index to a list of points inside the cell. By default, a DictionaryCellList is used.
• update_strategy = nothing: Strategy to parallelize update!. Available options are:
  • nothing: Automatically choose the best available option.
  • ParallelUpdate(): This is not available for all cell list implementations, but is the default when available.
  • SemiParallelUpdate(): This is available for all cell list implementations and is the default when ParallelUpdate is not available.
  • SerialUpdate()

References

• M. Chalela, E. Sillero, L. Pereyra, M.A. Garcia, J.B. Cabral, M. Lares, M. Merchán. "GriSPy: A Python package for fixed-radius nearest neighbors search". In: Astronomy and Computing 34 (2021). doi: 10.1016/j.ascom.2020.100443
• Markus Ihmsen, Nadir Akinci, Markus Becker, Matthias Teschner. "A Parallel SPH Implementation on Multi-Core CPUs". In: Computer Graphics Forum 30.1 (2011), pages 99–112. doi: 10.1111/J.1467-8659.2010.01832.X

ParallelUpdate()

Fully parallel update by using atomic operations to avoid race conditions when adding points into the same cell. This is not available for all cell list implementations, but is the default when available.

See GridNeighborhoodSearch for usage information.

PeriodicBox(; min_corner, max_corner)

Define a rectangular (axis-aligned) 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.

PolyesterBackend()

Pass as first argument to the @threaded macro to make the loop multithreaded with Polyester.@batch.

PrecomputedNeighborhoodSearch{NDIMS}(; search_radius = 0.0, n_points = 0,
                                     periodic_box = nothing, update_strategy = nothing)

Neighborhood search with precomputed neighbor lists. A list of all neighbors is computed for each point during initialization and update. This neighborhood search maximizes the performance of neighbor loops at the cost of a much slower update!.

A GridNeighborhoodSearch is used internally to compute the neighbor lists during initialization and update.

Arguments

• 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.
• n_points = 0: Total number of points. The default of 0 is useful together with copy_neighborhood_search.
• periodic_box = nothing: In order to use a (rectangular) periodic domain, pass a PeriodicBox.
• update_strategy: Strategy to parallelize update! of the internally used GridNeighborhoodSearch. See GridNeighborhoodSearch for available options.

SemiParallelUpdate()

Loop over all cells in parallel to mark cells with points that now belong to a different cell. Then, move points of affected cells serially to avoid race conditions. This is available for all cell list implementations and is the default when ParallelUpdate is not available.

See GridNeighborhoodSearch for usage information.

SerialUpdate()

Deactivate parallelization in the neighborhood search update. Parallel neighborhood search update can be one of the largest sources of error variations between simulations with different thread numbers due to neighbor ordering changes.

See GridNeighborhoodSearch for usage information.

ThreadsDynamicBackend()

Pass as first argument to the @threaded macro to make the loop multithreaded with Threads.@threads :dynamic.

ThreadsStaticBackend()

Pass as first argument to the @threaded macro to make the loop multithreaded with Threads.@threads :static.

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.

Keywords

• search_radius = 0.0: The fixed search radius. The default of 0.0 is useful together with copy_neighborhood_search.
• 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.
• periodic_box = nothing: In order to use a (rectangular) periodic domain, pass a PeriodicBox.

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 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.

# Template
nhs = GridNeighborhoodSearch{2}()

# Inside the simulation code, generate similar neighborhood searches
nhs1 = copy_neighborhood_search(nhs, 1.0, 100)

# output
GridNeighborhoodSearch{2, Float64, ...}(...)

foreach_point_neighbor(f, system_coords, neighbor_coords, neighborhood_search;
                       points = axes(system_coords, 2), parallel = true)

Loop for each point in system_coords over all points in neighbor_coords whose distances to that point are smaller than the search radius and execute the function f(i, j, x, y, d), where

• i is the column index of the point in system_coords,
• j the column index of the neighbor in neighbor_coords,
• x an SVector of the coordinates of the point (system_coords[:, i]),
• y an SVector of the coordinates of the neighbor (neighbor_coords[:, j]),
• d the distance between x and y.

The neighborhood_search must have been initialized or updated with system_coords as first coordinate array and neighbor_coords as second coordinate array.

Note that system_coords and neighbor_coords can be identical.

Arguments

• f: The function explained above.
• system_coords: A matrix where the i-th column contains the coordinates of point i.
• neighbor_coords: A matrix where the j-th column contains the coordinates of point j.
• neighborhood_search: A neighborhood search initialized or updated with system_coords as first coordinate array and neighbor_coords as second coordinate array.

Keywords

• points: Loop over these point indices. By default all columns of system_coords.
• parallel=true: Run the outer loop over points thread-parallel.

See also initialize!, update!.

initialize!(search::AbstractNeighborhoodSearch, x, y)

Initialize a neighborhood search with the two coordinate arrays x and y.

In general, the purpose of a neighborhood search is to find for one point in x all points in y whose distances to that point are smaller than the search radius. x and y are expected to be matrices, where the i-th column contains the coordinates of point i. Note that x and y can be identical.

See also update!.

update!(search::AbstractNeighborhoodSearch, x, y; points_moving = (true, true))

Update an already initialized neighborhood search with the two coordinate arrays x and y.

Like initialize!, but reusing the existing data structures of the already initialized neighborhood search. When the points only moved a small distance since the last update! or initialize!, this is significantly faster than initialize!.

Not all implementations support incremental updates. If incremental updates are not possible for an implementation, update! will fall back 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 points_moving = (true, false) in this case to avoid unnecessary updates. 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!.

@threaded x for ... end

Run either a threaded CPU loop or launch a kernel on the GPU, depending on the type of x. Semantically the same as Threads.@threads when iterating over a AbstractUnitRange but without guarantee that the underlying implementation uses Threads.@threads or works for more general for loops.

The first argument must either be a parallelization backend (see below) or an array from which the backend can be derived to determine if the loop must be run threaded on the CPU or launched as a kernel on the GPU. Passing KernelAbstractions.CPU() will run the GPU kernel on the CPU.

Possible parallelization backends are:

• PolyesterBackend to use Polyester.@batch
• ThreadsDynamicBackend to use Threads.@threads :dynamic
• ThreadsStaticBackend to use Threads.@threads :static
• KernelAbstractions.Backend to execute the loop as a GPU kernel

In particular, the underlying threading capabilities might be provided by other packages such as Polyester.jl.