Steffen Haug

One particular challenge of FRNN is that you don't know ahead of time how many neighbors a point has, so even just storing the set of neighbors compactly is a challenge. With GPUs not being particularily well suited for dynamically reallocating and growing buffers on the fly, to put it mildly, you need to allocate space up front, but on the other hand you don't actually know how much space you need.

Doing a full global sort of the particles sounds like lunacy, but in practice it can perform quite well, and it solves the issue of memory allocation. With a finite set of cells, you can sort the particles by cell index using a parallel counting sort. Cells can be assigned an index using a classic strided indexing: For an $m \times n$ grid and $\texttt{cell}(\vec x) = (i, j)$,

\begin{equation*} \texttt{index}(\vec x) = nj + i \end{equation*}

Sorting the points based on cell indices nicely arranges points in a cell contiguously in memory, effectively repurposing the particle buffer as several dynamically sized arrays for cells.

Let's work through an example. To start, we compute for each particle its corresponding chunk index:

Then, we perform a sort of the particles by cell index, and compute a lookup table to store the starts of the cells.

With this structure, you can find all the neighbor candidates of a point by computing its cell, offsetting the index to find adjacent cells, look up the slices in the lookup table, and iterate just over these slices. For example, to find the set of neighbors of $C$, the indices of adjacent cells are $\{ 1,2,3, 5,6,7, 9,10,11\}$, and by following the lookup table, we narrow down the set of neighboring points to the following set of slices:

This approach is simple and fast, but assigning a unique index in the lookup table to every cell means we are are restricted to a finite number of cells!

A common strategy for mapping a large set into a small table is using hashing. You might reasonably assume that we could simply replace the cell indices by a simple spatial hash, but the devil is in the details.

Approaches incorporating spatial hashing have already been implemented with some pretty nice results, but have failed to gain traction on the GPU.

Since you can't preallocate space for the buckets ahead of time, this hash table can not be implemented exactly as described, but sorting particles by a spatial hash is a nice starting point. This way, you could theoretically get the benefit of index sorting, namely repurposing the particle buffer as a set of dynamically sized arrays, and also the benefits of hashing, namely mapping an infinite simulation domain into a finite table. Compared to index sorting, you only pay the extra overhead of calculating a simple hash.

The problem is that hashes can collide within the neighborhood of a particle, which means that the set of neighbors can not be found by naively traversing the slices without some additional bookkeeping to deal with hash collisions, or we risk counting particles more than once.

As an example, consider the following example assignment of hash values to cell indices:

\begin{equation*} \left\{\hspace{.5em} \begin{array}{lccl} \texttt{hash} : & 2 & \longmapsto & 3 \\ & 3 & \longmapsto & 1 \\ & 4 & \longmapsto & 0 \\ & 6 & \longmapsto & 2 \\ & 7 & \longmapsto & 5 \\ & 8 & \longmapsto & 2 \\ & 9 & \longmapsto & 3 \end{array} \right. \end{equation*}

To compute the hash of a point, we just apply $\texttt{hash}$ to its cell index.

Now notice, particles in cells $2$ and $9$ both hash to $3$, so the bucket with index $3$ in the hash table will contain the set $\{ B, H, I \}$. Now, consider searching for neighbors of $C$:

If we naively traverse $C$s adjacent cells $\{ 1,2,3 , 5,6,7 , 9,10,11 \}$ and compute their hash, the set $\{ B, H, I \}$ will correspond to two adjacent grid cells – both $2$ and $9$:

If this is not deduplicated, and we just move on to iterating over the slice corresponding to each adjacent cell, these particles will in turn be counted twice.

Obviously, this is bad.

This imposes deduplication further down the line, or in the worst case, if we don't do anything about it, we have a subtle bug where in rare cases points may contribute twice to a mean, or particles may contribute forces twice to the pressure gradient, and so on.

This deduplication is fairly trivial on the CPU; there are any number of ways to solve it: You could for example store the cell hashes in a hash-set, or do a quick scan of the hashes visited prior and skip repeat visits, to name two simple solutions. Crucially, on the CPU threads don't pay for the bookkeeping done on another, so simple bookkeeping is cheap. Not so on the GPU: Any test to see if a hash leads to a repeat visit, or comparison of hashes to remove duplicates is an opportunity for the program to branch, causing serialized execution. Highly divergent branchy loopy code should be avoided if possible.

This seems like quite the pickle.

It would be best, it seems, if we could avoid hash collisions specifically in the neighborhood of a particle.

For $d = 2$, $M=3^2$, and a hash table of size $K = MN$, a hash function can be defined as follows:

\begin{equation*} \left\{\hspace{.5em} \begin{array}{lccl} \texttt{hash} : & \mathbb R^2 & \longrightarrow & \mathbb Z_K \\ & \vec x & \longmapsto & M \beta + \kappa \end{array} \right. \end{equation*}

where $\kappa$ is the equivalence class, and

\begin{equation*} \beta = (P_1 \, i \;\texttt{xor}\; P_2 \, j) \;\texttt{mod}\; N \end{equation*}

where $P_1$ and $P_2$ are relatively large prime numbers and the $\texttt{xor}$ is taken bitwise, selects a "bucket" in the hash table within the class $\kappa$. The reason this works is that adding a multiple of $M$ to $\kappa$ does not change its value $\texttt{mod}\; M$, and thus $\texttt{hash}(\vec x)$ will preserve the cyclic structure of $\kappa$!

For example, looking at the same neighborhood as before,

we see that indeed, the cells in this neighborhood does end up in distinct "columns" of the hash table:

Thus, there is no need to do any bookkeeping to deduplicate hash collisions between cells. Every cell adjacent to a point will always hash to a distinct class, and thus to a distinct bucket, and thus correspond to a distinct slice of the underlying sorted buffer of particles; there is no longer any risk of scanning the same slice twice when iterating over adjacent grid cells.

The size of the hash table can be changed freely by modifying $N$. This allows you to tune the hash collision rate depending on performance and memory constraints, while retaining locally perfect hashing.

In other words, a hash-based sorting scheme with this kind of hash function actually works.

The hash function determines the order of the particles in memory, and as such the choice would have a pretty high impact on performance. The hash function should ideally be chosen with the GPU kernels' memory access patterns in mind, to minimize uncoalesced reads and writes.

The hash function I presented here is almost certainly not optimal with respect to memory access patterns, it is simply chosen because it demonstrates the point in a simple manner and lends itself nicely to illustrations.

I reckon optimizing the interaction between how you structure the kernels implementing the table construction and how the hash function scatters particles in memory is an interesting topic to investigate in its own right, and maybe a good starting point for a future post.