CellLists

VortexPasta.CellLists Module

CellLists

Module implementing the cell lists algorithm over $N$-dimensional periodic domains.

See the Wikipedia article for some details.

source

Types

VortexPasta.CellLists.PeriodicCellList Type

PeriodicCellList(
    [backend = CPU()],
    rs_cut::NTuple{N, Real}, periods::NTuple{N, Real},
    [nsubdiv = Val(1)];
) -> PeriodicCellList{N}

Construct a cell list for dealing with pair interactions in N dimensions.

The cutoff distances rs_cut (which can be different in each direction) don't need to exactly divide the domain period L into equal pieces. In that case the cutoff distances will be adjusted (slightly increased) to fit the domain period.

Note that this cell-lists implementation can identify pairs which are slightly beyond the given cut-off distance. For performance reasons, such pairs are still returned, and it's up to the user to see whether such pairs should still be computed or not depending on whether a strict cut-off distance is wanted.

Optionally, one can choose to subdivide each cell (of size ≈ rcut) onto nsubdiv subcells. This can significantly improve performance, since it allows to discard some spurious pair interactions (i.e. beyond the chosen cutoff radius) as described here. In practice, a value of 2 or 3 can significantly improve performance compared to no subdivision (1). For convenience, it can be passed as nsubdiv = Val(M) or as nsubdiv = static(M).

Infinite non-periodic domains (in the sense of period = Infinity()) are not supported.

GPU usage

To run on GPUs, pass a KernelAbstractions backend such as CUDABackend or ROCBackend as the first argument.

If multiple GPUs are available, make sure to activate the device where computations should be performed before constructing a PeriodicCellList. For example:

using KernelAbstractions: KernelAbstractions as KA
backend = ROCBackend()
device_id = 2  # assuming there are two or more available GPUs
KA.device!(backend, device_id)
cl = PeriodicCellList(backend, ...)

source

Functions

VortexPasta.CellLists.set_elements! Function

CellLists.set_elements!([get_coordinate::Function], cl::PeriodicCellList, xp::AbstractVector; folded = Val(false))

Set all elements of the cell list.

In general xp is a vector of spatial locations.

Optionally, a get_coordinate function can be passed, which will be applied to each xp[j] in order to obtain a coordinate. By default this is identity.

Moreover, one may pass folded = Val(true) in case all points in xp have been folded into the main periodic cell (i.e. x⃗ ∈ [0, L]ᵈ), which can speed-up computations. This is a dangerous option: if points haven't been really folded, this lead to wrong results or crashes.

This function resets the cell list, removing all previously existent points.

source

VortexPasta.CellLists.foreach_pair Function

CellLists.foreach_pair(
    f::Function, cl::PeriodicCellList, xp_dest::AbstractVector;
    batch_size = nothing,
    folded = Val(false),
)

Iterate over all point pairs within the chosen cut-off distance.

A pair is given by a destination point xp_dest[i] and a source point xp[j], where xp is the vector of coordinates that was previously passed to CellLists.set_elements!. Note that the vector of destination points xp_dest can be different than that of source points.

The first argument should be a user-defined function f(x⃗, i::Integer, j::Integer) taking a destination point x⃗ = xp_dest[i] and a pair of destination/source indices (i, j). Typically, one wants to modify one or more destination vectors at index vp[i], which can be done within this function.

On CPUs this function is parallelised using threads. Parallelisation is done at the level of all destination points xp_dest, ensuring that only a single thread can write to a location vp[i] (in other words, atomics are not needed). This function also works on GPUs.

One may set folded = Val(true) if points in xp_dest have been previously folded. See set_elements! for more details.

This function should be called after CellLists.set_elements!.

See also CellLists.foreach_source.

Iterating over batches

Similarly to foreach_source, one can iterate over batches of (up to) W source points at a time (for each single destination point x⃗ = xp_dest[i]). For this one should pass batch_size = Val(W) as a keyword argument. Note that in this case, the user-defined function f is expected to have the signature f(x⃗, i::Integer, js::NTuple{W}, m::Integer). See foreach_source for details on the last two arguments.

Example usage

julia> using StaticArrays: SVector

julia> using LinearAlgebra: ⋅  # dot product

julia> r_cut = 0.4; rs_cut = (r_cut, r_cut, r_cut);

julia> Ls = (2π, 2π, 2π);

julia> nsubdiv = Val(2);

julia> cl = PeriodicCellList(CPU(), rs_cut, Ls, nsubdiv)
PeriodicCellList{3} with:
 - backend:         CPU(false)
 - number of cells: (31, 31, 31)
 - period:          (6.283185307179586, 6.283185307179586, 6.283185307179586)
 - cell size:       (0.2026833970057931, 0.2026833970057931, 0.2026833970057931)
 - effective r_cut: (0.4053667940115862, 0.4053667940115862, 0.4053667940115862)
 - number of subdivisions: 2

julia> xp = [rand(SVector{3, Float64}) .* Ls for _ in 1:1000];  # create source points in [0, L]³

julia> CellLists.set_elements!(cl, xp);  # process source points

julia> xp_dest = [rand(SVector{3, Float64}) .* Ls for _ in 1:200];  # create destination points in [0, L]³

julia> vp_dest = zeros(length(xp_dest));  # vector where "interactions" will be written to

julia> CellLists.foreach_pair(cl, xp_dest) do x⃗, i, j
           # Compute some "influence" of xp[j] on x⃗ == xp_dest[i]
           # Note that it is safe to write to vp_dest[i] without atomics, even when running in parallel.
           @inbounds vp_dest[i] += x⃗ ⋅ xp[j]
       end

source

VortexPasta.CellLists.foreach_source Function

CellLists.foreach_source(f::Function, cl::PeriodicCellList, x⃗_dest; batch_size = nothing, folded = Val(false))

Iterate over source points which are close to a given destination point x⃗_dest.

This is similar to CellLists.foreach_pair but only works on a single destination point. On the CPU, foreach_source is not parallelised. One may want to parallelise at the level of multiple destination points to match the behaviour of foreach_pair.

The first argument should be a user-defined function f(j::Integer) taking the index j of a source point.

Performance-wise, this function seems to give similar performance to foreach_pair on CPUs. On GPUs one may prefer to use foreach_pair which is more adapted for GPU computing.

One may set folded = Val(true) if x⃗_dest has been previously folded. See set_elements! for more details.

Iterating over batches

For performance reasons (e.g. to enable SIMD), one may want to iterate over batches of source points in cl. For this, one should pass the keyword argument batch_size = Val(W) where W is the wanted batch size.

When batching is enabled, the user-defined function f is expected to have the signature f(js::NTuple{W}, m::Integer), where js are the indices of W nearby source points, and m ∈ 1:W is the number of "valid" points. In other words, indices js[(m + 1):W] should be discarded (or masked out) in the f function. Note that, even though these indices should be ignored to get correct results, these are guaranteed to be valid indices in 1:Np (where Np is the number of source points), and thus can be safely used to access source coordinates (which is convenient in the context of SIMD).

source

Base.empty! Method

Base.empty!(cl::PeriodicCellList) -> cl

Remove all elements from the cell list.

source

VortexPasta.CellLists.nearby_elements Function

nearby_elements(cl::PeriodicCellList{N}, x⃗)

Return an iterator over the elements that are sufficiently close to the point x⃗.

The iterator returns the elements which are likely to be within the cutoff radius of the point x⃗. More precisely, it returns elements in the same cell as x⃗ as well as in neighbouring cells.

Here x⃗ should be a coordinate, usually represented by an SVector{N} or an NTuple{N}.

source