Approach

The CUDA implementation parallelizes the same simulation model by keeping the full simulation state on the GPU and decomposing each timestep into kernels that operate over particles, cells, or reductions. The final goal is to preserve the semantics of the CPU pipeline while replacing the serial loops with SIMT execution. Most kernels use one thread per particle. Some shared global quantities, such as the adaptive timestep and die reductions, are handled with reductions or atomic accumulation. Importantly, the only CPU <-> GPU data transfers after initialization (which is done entirely on the CPU) are to send control inputs to the GPU, and send particle information (minimum required) back to the CPU for rendering.

Initialization

This project is mostly concerned with asymptotic speedup for long-running simulations, and so we chose to exempt initialization from our GPU-acceleration efforts. The most important portion of initialization when understanding the GPU fork (which leads to longer initialization times on the GPU) is that the CPU, upon computing the full initial state, uploads it all to the GPU.

GPU State and Memory Layout

The GPU implementation mirrors the CPU structure of the array design with persistent device buffers for each of the arrays present in the CPU state. Initialization uploads the arrays once, allocates CUDA-side runtime buffers, uploads simulation parameters, and stores the die into the device's memory. After initialization, steady state behavior avoids copying the full particle state back to the host. Host-Device transfers are only used for diagnostics and control values.

Parallel Ghost maintenance

Ghost particles are rebuilt on the GPU to avoid round-tripping particle positions. The CUDA path compacts real fluid particles into the prefix of arrays, removing old ghosts. It then marks all fluid particles within one smoothing radius of the sphere wall, uses an exclusive scan to compute each new ghost's destination offset, scatters new ghost particles into the tail of the arrays, and writes the new particle count to device memory.

Parallel Spatial indexing and neighbor construction

We use spatial hashing to avoid checking every particle against every other particle every step. Since our particles are inside a small fixed sphere, we can use a fixed-length grid instead of a general hash table. Each particle is mapped to a grid cell, the cell coordinate is biased so negative positions become valid indices, and that cell index becomes the particle’s key.

On the GPU, we compute these keys in parallel, sort particles by key, and store the start/end range for each grid cell. Then each particle only searches particles in nearby cells contiguously when building its neighbor list.

We are allowed to do this because our simulation domain is bounded and known ahead of time. A normal spatial hash needs to handle large or unbounded space and hash collisions. Ours does not, so we can trade generality for speed.

Parallel Die Coupling and Reductions

We also parallelize the die coupling with one thread per particle. Each thread checks whether its particle is a real fluid particle, evaluates the die-signed distance contact condition, computes the force and drag forces when applicable, and adds the resulting fluid force to that particle. The die reaction is a shared quantity, and so we can't just use simple independent stores. So, we use block-local accumulation and global atomic additions to combine the particle reaction forces and torques into device-side die accumulators.

Adaptive time is done by having a kernel compute per block maxima for speed and acceleration, then we use a scalar reduction to combine block maxima, and then use a finalization kernel to evaluate each time step limit to get our real time step limit.

Parallel integration

Fluid integration is embarrassingly parallel once forces and timestep are known. A kernel updates velocity and position for each fluid particle using the given time step and force information. Then we have a small kernel to swap the velocity buffers. Then, finally, we have a boundary kernel check each particle against the sphere radius and apply the clamp, restitution, and damping logic as with the CPU implementation.

The rigid die update is done after the accumulation has already been loaded into device buffers. We therefore just have to use an integration kernel to apply caps, gravity, buoyancy, damping, linear and angular updates, and quaternion normalization. Then, a kernel keeps the die within the spherical container by clamping it. Although there isn't as much parallelism in these kernels as in the particle kernels, we keep them on the device to avoid unnecessary synchronization and host-device transfers during time steps.

GPU Pipeline

Iterative Optimization Journey

Our first GPU-accelerated fork was done in two stages. First, we converted the compute_density, compute_pressure, compute_forces, integrate, handle_boundary functions into kernels. The speedup is attached below.

Round 1a

The speedup was incredibly modest because at every step, there was a lot of CPU-GPU data transfer. This was fixed when we moved the rest of our functions to CUDA and sought to eliminate all data transfers (more on this later). This led to a much improved speedup. We also realized while combining our code that we were having our ping ponged buffers set up to be vel[i][ping/pong] and so we flipped this so that each warp gets a contiguous block of memory to update (this was not done on the CPU as it slowed our CPU implementation).

All GPU

To note, we recorded the speedup below on the exact same code.

Bad bad bad bad bad so bad

We realized that our benchmarking metrics were introducing a lot of noise in our speedup calculations. We mainly care about how quickly we can produce new sim frames at steady state, not really how long it takes to initialize the data structures. By including initialization time in the speedup calculation, we were allowing for differences in initialization time to appear as large speedup gains, even though they were not. We saw this issue a lot when we started optimizing to the point where each step was very fast, because then any changes in initialization time represented a much larger portion of overall runtime, making the output noisier. To address this, we excluded init time from future speedup calculations, since it takes a nondeterministic amount of time (Poisson disk sampling is random), and we care more about frames per second than overall time.

Kernel Execution Profiling Results

The profiling information told us that we should focus on optimizing kernel_build_neighbor_lists, since it took up 87% of the runtime. We realized that the hash-based method led to poor spatial locality as it sent contiguous memory accesses far apart from each other in location. The reason people do spatial hashing is to enforce load-balanced bins, but the nature of our sim leads to natural load balancing throughout the sphere, and so we used a simple geometric binning approach. We also discovered a lot more allocation of memory than we anticipated. Eventually, we found out that it was the Thrust sorting function we were using that was allocating buffers under the hood every time we called it. To fix this we switched over to CUB for sorting, since we were able to pre-allocate a single buffer and reuse it for every step, which significantly improved performance. These optimizations led to our speedup imrpoving from around 7× to 13×.

After profiling again, we noticed that there was a ton of data being transferred with cudaMemcpyToSymbol. Confused, we looked through our code to find any data transfers being done in the main step loop. We realized that we were passing the number of ghost particles back and forth between the CPU and GPU every step to calculate the total number of particles. We were also reducing particle forces onto the die using a partial reduction on the GPU, and then doing the final reduction on the CPU. This required that all the partial reductions be sent back to the CPU, just for a simple sum to be computed, and then to send all the data back, many times per step. This was extremely inefficient, and to fix it, we decided to keep all the sim state only on the GPU. We wanted to prevent the CPU from doing any computations if possible. We did this by replacing some of the CPU operations with 1x1 kernels that did the same thing without having to move the data. This resulted in a massive speedup, bringing us to 65×.

No more data transfer!

This was the first version of the simulator that was able to simulate particle motion at 60 Hz, which was our target for the real-time demo. However, we were not able to actually render the frames in real time, since we were still using a matplotlib-based renderer, which was too slow. So we switched to a C++ based renderer and were able to stream the sim in real time.

We then tried a few more smaller rounds of optimizations, such as reducing the number of kernel launches by consolidating neighboring kernels into larger ones, and folding the 1x1 kernels to be within other kernels, but these produced modest speedups. We realized that any significant speedup would be found in kernel_build_neighbor_lists, since that still took up around 45% of the runtime

With this in mind, we went back to build_neighbor_lists and looked for anything else we could do to squeeze out more performance. We ended up finding a small debug if statement that verified that the cell key was the same as the sorting key when building the neighbor list. This was a holdover from when we were first implementing build_neighbor_lists, and would never actually trigger under normal operation. Just by removing this redundant check, we were able to get a speedup (+20x), even though it didn't affect the algorithmic complexity at all. This brought our overall speedup to 96×.

Final speedup

Another idea we had to reduce time spent building neighbor lists was to only rebuild every k steps. We know that particles don't move much between steps (we guarantee this by reducing the $\Delta t$ based on the max particle speed. So we tested different semi-static building intervals to 2, 4, and 8 to see how much it affected performance. Unsurprisingly it had a big impact, but we were more curious how it affected correctness. If a neighboring particle moves outside of a neighboring cell without it rebuilding, then the particle interaction will be missed. We experimented with dynamic rebuilding intervals based on when particles actually leave the cell, but we found that the overhead of tracking particle position erased any gains that we would get by rebuilding less frequently.

Afterwards, we ran benchmarks on our final implementation to see the effects of each optimization.

Effect of each optimization sorted by effect

Results

In our results we observe that running the sim for more timesteps leads to linear growth. This is to be expected as the work per time step is constant, no matter how many time steps we take.

We measured our performance with computation time (time per sim step) against a serial CPU baseline. We did not modify the CPU implementation across our various improvements so that it would remain a fair baseline. We did not consider initialization time very much since the goal of the project was to have a sim that can run in real time, and the user experience is much more affected by the smoothness and frame rate of the simulation in stead-state rather than how long it takes for the sim to start up. By virtue of maintaining a static baseline, increases in speedup are directly proportional to decreases in clock time and vice versa.

The size of our inputs were built to simulate a real 8-ball. We sized the 8 ball accordingly, and the die was also sized to approximate a real 8 ball (albeit with a cubic die). From here, we tuned our physics parameters to make the sim stable and have realistic behavior, and simply maintained these parameters while we optimized. We've attached relevant parameters below.

// Smoothing and discretization
constexpr float H          = 0.05f;
constexpr float PARTICLE_R = H / 2.0f;
constexpr float PARTICLE_M = 0.008f; // kg
constexpr float SPHERE_R   = 0.15f;  // inner radius of sphere
constexpr float DEFAULT_FLUID_DENSITY = 1000.0f; // kg/m^3
constexpr int   N_PARTICLES           = 3500; // target count

From this, our experiements simply ran the sim for some number of timesteps (which we swept across).

Problem Size Sweep

It is worth noting that the number of particles (3500) was chosen because that results in the most stable and realistic looking sim, not because it gave the best speedup (it doesn't). Our speedup would be much greater if we boosted the number of particles in order to get better hardware utilization on the GPU and a higher proportion of parallelizable work. When we tested our GPU implementation with higher particle counts we saw that the speedup is linear with the number of particles.

On the matter of the execution profile, our intiialziation time for the CPU with our default problem size is usually around 2.37 seconds, and for the GPU it is around 2.94 seconds. The extra GPU time is spent copying the buffers onto the GPU. From here, we simply run a series of kernels one after the other, so the rest of the algorithm's execution time can be simply broken down by the nsys profile.

Kernel Execution Profiling Results After Optimization

The final profiling table shows that particle-particle (kernel_build_neighbor_lists and kernel_compute_forces) interaction still takes up the majority of runtime.

Nsight Compute Metrics for Neighbor List Construction

NCU for Force Computation

Profiling with NCU gave us a lot of insight into what was limiting performance. First, we noticed that active warps were only around 20%, which means that a lot of warps were sitting idle. This confirms our hypothesis that the problem size was too small to get full utilization of the GPU. We also notice that load balancing per thread is not ideal: we have an average of 20/32 useful threads per warp. This is probably because the building of the neighbor list and the calculation of the force are not perfectly uniform. Different particles have different numbers of neighbors; some are ghosts, and the distribution of the particles that actually pass the radius check will affect coalescence unpredictably.

We notice that the cache hit rate is very good for building neighbor lists (> 95%). This is because we sort the neighbor lists before traversing, resulting in good spatial locality. The hit rate is a bit lower for force computation, because interacting neighbors are not guaranteed to be adjacent in the neighbor list, and looking up the properties of neighboring particles requires a layer of indirection (indexing into the state array by particle id), which has poor spatial locality. We could address this problem by keeping the particle state parameters nearby in memory, potentially by sorting parameter arrays on the same key as the neighbor list. That way, when we traverse neighbors, the parameters of neighboring particles are more likely to be nearby in memory when we access them, at the cost of another sorting phase.

Soundness of target machine

Our choice of machine was sound. The problem is inherently simulating a massively data parallel simulation, and so a GPU is a very natural choice for this. While our problem size was likely too small to get the maximum speedup from our machine, the alternatives are likely much worse. One that is somewhat interesting would be a larger supercomputer with hundreds of thousands of threads, although we anticipate that for metrics like energy or cost, a GPU-based implementation would far outpace a super computer implementation.

One thing that could have been interesting and maybe led to a different target machine would be optimizing this for training (ML) robot arms to rig magic 8 balls. In which case, the question becomes how do we parallelize many 8 balls and also reset them quicker on machines with multiple GPUs. All this to say, the GPU was definitely the best choice, and I have some interesting projects ahead of me.

References

1. Matthias Müller, David Charypar, and Markus Gross. Particle-based fluid simulation for interactive applications. In Proceedings of the 2003 ACM SIGGRAPH/Eurographics Symposium on Computer Animation (SCA), pages 154--159, 2003.
2. Joe J. Monaghan. Smoothed particle hydrodynamics. Annual Review of Astronomy and Astrophysics, 30(1):543--574, 1992.
3. Matthias Teschner, Bruno Heidelberger, Matthias Müller, Danat Pomerantes und Markus Gross. Optimized spatial hashing for collision detection. In Vision, Modeling and Visualization (VMV), volume 3, pages 47--54, 2003.
4. Hagit Schechter and Robert Bridson. Ghost SPH for animating water. ACM Transactions on Graphics (SIGGRAPH), 31(4):1--8, 2012.
5. Alejandro JC Crespo, Moncho Gómez-Gesteira y Robert A Dalrymple. Boundary conditions based on dynamic particles in SPH method. International Journal for Numerical Methods in Engineering, 71(13):1537--1588, 2007.
6. Markus Becker and Matthias Teschner. Weakly compressible SPH for free surface flows. In Proceedings of the 2007 ACM SIGGRAPH/Eurographics Symposium on Computer Animation (SCA), pages 209--217, 2007.
7. Nadir Akinci, Markus Ihmsen, Gizem Akinci, Barbara Solenthaler and Matthias Teschner. Versatile rigid-fluid coupling for SPH. ACM Transactions on Graphics (SIGGRAPH), 31(4):1--8, 2012.