Particle–Cell Operations Optimization: Algorithm, GPU & CUDA

Years ago, I was working on the Uintah software (a large-scale multiphysics simulation framework developed at the University of Utah). My focus was on optimizing the performance of the software by code profiling, locating the bottlenecks and improving the performance. One of the main bottlenecks I identified was the particle–cell coupling, which is the essential part of the two-way interaction between the Lagrangian particle phase and the Eulerian gas phase.

We presented this work at multiple conferences (SIAM CSE 2017 & NC 2017), however, never published a full article on it. This is a draft on how we approached the Particle-Cell operation problem and addressed it.

Background: Uintah Software

The Uintah Computational Framework is an open-source, massively parallel software suite designed for simulating complex fluid-structure interaction, combustion, and fracture mechanics on supercomputers. Originally funded by the Department of Energy, it manages automated load balancing and asynchronous task execution across hundreds of thousands of CPU and GPU cores.

We designed our own plugin named "Wasatch" for gas and particle simulation. For the gas phase, we used Eulerian methods, while the particle phase was treated in a Lagrangian manner. We used two-way coupling, meaning gas phase and particle phase influence each other reciprocally. But, particles do not interact with each other, as in our simulation the particle-particle interaction is negligible.

After heavy code profiling, I identified multiple places that were slowing down the simulation, and a couple of them were easy targets for optimization, and particle-cell interaction was one of them!

Simplifying Particle-Cell Interaction

In every time step, the simulation must do these calculations:

• Interacting Particles and Cells: The first step is to identify which particles are interacting with which gas-phase cells, meaning determining the spatial location of each particle and mapping it to the corresponding cell in the computational grid.
• Particle-to-Cell (P2C): This means the source terms that particles release (or absorb) into the given gas-phase cells. For example, in our coal particle combustion simulation, particles release heat and combustion products (CO, CO₂) into surrounding gas-phase cells.
• Cell-to-Particle (C2P): This means the impact of gas-phase cells on interacting particles. For example, in our coal particle combustion simulation, gas-phase cells affect the motion (drag) and heating of the particles.

Technically, in every time step, in addition to the conservation equations which we are solving for gas and particles, the code must account for the interactions between particles and gas cells, which is incorporated in the source terms.

Note: In our code, we used two-way coupling, meaning that gas has an impact on particles, particles have an impact on cells, and there is no particle-particle interaction.

To make “every particle” concrete, here is one of the simulations from the original talk. This is what the particle phase actually looks like: a cloud of Lagrangian particles carried through the domain by the gas, their positions evolving every time step:

Particle-Cell Interaction

The original code, and the most scientifically accurate approach, is volume-weighted interpolation: a particle sitting between cells distributes its contribution in proportion to how much of its volume lies in each cell.

However, given the fact that the particle size is much smaller than the cell size (0.1x), it is reasonable to assume that each particle interacts with only the cell that its center is located in. This assumption has no or negligible impact on the accuracy of the simulation while significantly reducing the computational cost and code complexity.

Code complexity is a very important factor, because if we want to move the computation to GPU, as we will do in the next section, the simpler algorithm will reduce the chance of divergence.

for( ; isrc != ise; ++ipx, ++ipy, ++ipz, ++ipsize, ++isrc ){

  const double rp = *ipsize * 0.5;

  // Identify the location of the particle boundary
  // (assuming that it is a cube)
  const double pxlo = px_ ? *ipx - rp : 0;
  const double pylo = py_ ? *ipy - rp : 0;
  const double pzlo = pz_ ? *ipz - rp : 0;

  const double pxhi = px_ ? pxlo + *ipsize : 0;
  const double pyhi = py_ ? pylo + *ipsize : 0;
  const double pzhi = pz_ ? pzlo + *ipsize : 0;

  const int ixlo = ( pxlo - xloface ) / dx_;
  const int iylo = ( pylo - yloface ) / dy_;
  const int izlo = ( pzlo - zloface ) / dz_;

  // hi indices are 1 past the end
  const int ixhi = px_ ? std::min( nmax[0], int(( pxhi - xloface ) / dx_) + 1 ) : ixlo+1;
  const int iyhi = py_ ? std::min( nmax[1], int(( pyhi - yloface ) / dy_) + 1 ) : iylo+1;
  const int izhi = pz_ ? std::min( nmax[2], int(( pzhi - zloface ) / dz_) + 1 ) : izlo+1;

  const double pvol = (px_ ? 2*rp : 1)
  * (py_ ? 2*rp : 1)
  * (pz_ ? 2*rp : 1);
  for( int k=izlo; k<izhi; ++k ){
    const double zcm = zloface + k*dz_;
    const double zcp = zcm + dz_;
    // determine the z bounding box for the particle in this cell
    const double zcont = pz_ ? std::min(zcp,pzhi) - std::max(zcm,pzlo) : 1;
    for( int j=iylo; j<iyhi; ++j ){
      const double ycm = yloface + j*dy_;
      const double ycp = ycm + dy_;
      // determine the y bounding box for the particle in this cell
      const double ycont = py_ ? std::min(ycp,pyhi) - std::max(ycm,pylo) : 1;
      for( int i=ixlo; i<ixhi; ++i ){
        const double xcm = xloface + i*dx_;
        const double xcp = xcm + dx_;
        const double xcont = px_ ? std::min(xcp,pxhi) - std::max(xcm,pxlo) : 1;
        // contribution is the fraction of the particle volume in this cell.
        const double contribution = xcont*ycont*zcont / pvol;
        dest(i,j,k) += *isrc * contribution;
      }
    }
  }
} // particle loop

for( ; isrc != ise; ++ipx, ++ipy, ++ipz, ++isrc ){
  const int i = ( *ipx - xloface ) / dx_;
  const int j = ( *ipy - yloface ) / dy_;
  const int k = ( *ipz - zloface ) / dz_;
  dest(i,j,k) += *isrc;
} // particle loop

Simplifying the particle–cell operators to this nearest-cell approximation yields roughly a 2× speedup on CPU.

Two patterns in this data are worth noting. The gather starts near 2.6–2.8× on small meshes but drifts down toward 1.7× at 128³: as the mesh grows the cost is dominated less by the arithmetic the approximation eliminates and more by memory access, which it doesn’t change. And for the scatter, more particles per cell means more speedup — the exact method’s overlap computation is per-particle work that the approximation removes entirely. That growing role of memory access foreshadowed the real fight on the GPU.

GPU Implementation

With the simplified algorithm implemented and tested, the next step was to move both operators (P2C and C2P) onto the GPU. The benchmarks below are measured on cubic meshes from 16³ to 128³ cells, with particles randomly distributed in the domain, as they are in a real simulation, at average loadings of 0.1, 1, and 10 particles per cell (ppc). Speedup is GPU time versus a single CPU core, for the approximate method on both sides.

Three things stand out in the data:

• Small problems don’t pay for the GPU. At 16³ with 0.1 ppc the GPU is actually slower than the CPU (0.3–0.4×) — there simply isn’t enough work to amortize kernel launch and keep the device occupied. The crossover comes quickly: by 64³ every configuration is 10× or better.
• Speedup grows with work. More cells and more particles per cell both push the numbers up, peaking around 33–41× for the largest workloads.
• Scatter consistently trails gather. At 128³ C2P reaches ≈ 41× while P2C peaks around 15–22×. That gap is structural, and worth unpacking.

The challenges

Using CUDA, I have transferred both the particle and gas-phase equation solvers to the GPU. Following the Eulerian approach for the gas-phase, each cell is mapped to a GPU thread, whereas for the particles (Lagrangian approach), each particle is mapped to a GPU thread. After the code implementation, I used nvprof and NVIDIA Visual Profiler (now being replaced by Nsight products) to profile the performance of the GPU kernels and found a couple of places where it was slowing the code down.

Atomic Operations

Using the simplified algorithm, it guarantees that each particle interacts with only one cell at a time. However, we have multiple particles per cell. This means that in the P2C operation multiple particles can try to update the same cell simultaneously. In order to resolve this issue, atomic operations are required to ensure that updates to the same cell are performed safely and correctly. This may introduce a performance reduction, as threads must wait for one particle to finish updating the cell information, for the next particle (that is in the same cell) to apply the changes. The orange lines in Fig. 5 show the instances where we have multiple particles in a cell simultaneously, where atomic operations are necessary to ensure correct updates.

Coalesced Memory Access

The second problem, especially for the P2C operation, is memory access patterns. Particles are randomly distributed across the grid, and are interacting with the cells. Imagine P₁ is interacting with cell C₃; however, P₂, which is adjacent to P₁ in memory, might be interacting with a completely different cell, say C₅₀. This forces the GPU to do a separate memory transaction (read the cell information) for each particle.

Suggestions

Coalesced memory access was really slowing our code down, as it had a significant impact on the performance of both P2C and C2P operations. One of the suggestions we came up with was to sort the particles by their cell indices to improve memory coalescing and thus enhance the performance of both P2C and C2P operations. This approach introduces an additional sorting step, which itself has a computational cost. Also, our main simulations, which were coal particle combustion, were mostly very particle-dispersed simulations, which made it less likely for sorting by cell indices to yield significant memory coalescing benefits.

However, if you are dealing with a particle-laden simulation, sorting the particles by cells can lead to significant improvements in memory coalescing and overall performance. The key is that you do not want to do the sorting at every timestep, but rather at intervals that balance sorting overhead with performance gains. A static equation for the sorting interval, determined at the initialization of the simulation, can be suggested as the following:

where u_x,char, u_y,char, and u_z,char are representative initial particle velocity magnitudes in each direction. A high percentile, such as the 95th percentile of the initial particle velocities, can be used to avoid making the interval overly sensitive to individual outliers.

The fixed number of timesteps between sorting operations is then:

where α is a safety factor, typically between 0.25 and 0.50. The value of N_sort is calculated once at initialization and remains constant throughout the simulation.

Particles are therefore sorted at:

I didn't have a chance/time to test this idea, as I moved to the next bottleneck in the simulation, which was radiation properties calculation (need another blog post).