Parallel Stencil Computation

Parallel Stencil Computation is the structured-grid numerical technique where each grid point's value is updated based on a fixed pattern of neighboring values — fundamental to finite-difference methods in CFD, weather simulation, seismic imaging, and image processing — where the regular access pattern enables highly efficient GPU and multi-node parallelization through domain decomposition with halo exchange, achieving 50-80% of peak memory bandwidth on modern hardware when properly optimized with tiling, vectorization, and temporal blocking.

Stencil Pattern

A stencil operation updates point (i,j,k) from its neighbors:

u_new[i][j][k] = c0*u[i][j][k] +
  c1*(u[i-1][j][k] + u[i+1][j][k]) +
  c2*(u[i][j-1][k] + u[i][j+1][k]) +
  c3*(u[i][j][k-1] + u[i][j][k+1]);

This 7-point 3D stencil (Jacobi/Laplacian) reads 7 values and writes 1. Arithmetic intensity: 7 FLOPS / 8 memory accesses × 4 bytes = 0.22 FLOPS/byte — severely memory-bandwidth-bound.

Parallelization Strategy

Domain Decomposition: Divide the 3D grid into sub-domains, assign one to each processor/GPU. Each sub-domain is updated independently for interior points. Boundary points require neighbor data from adjacent sub-domains → halo exchange.

Halo Exchange: Before each time step, each processor sends its boundary layer to neighbors and receives their boundary layers:
- 3D domain with P processors: each processor exchanges 6 faces (±x, ±y, ±z).
- Communication volume: proportional to surface area of sub-domain.
- Computation: proportional to volume of sub-domain.
- Surface-to-volume ratio decreases with larger sub-domains → strong scaling limited by communication.

GPU Stencil Optimization

• Thread Mapping: One thread per grid point. 3D thread blocks (e.g., 32×8×4) map to 3D grid regions. Adjacent threads access adjacent memory → coalesced global memory reads.
• Shared Memory Tiling: Load a tile (including halo) into shared memory. Compute stencil from shared memory (fast, reusable). For a 7-point stencil with tile 32×32: load 34×34 into shared memory (32+2 halo in each dimension). Interior reads hit shared memory instead of L2/global.
• Register Tiling (2.5D Blocking): Load one z-plane into registers, compute stencil using current + cached previous/next planes. Walk through z-dimension, sliding the register window. Reduces shared memory pressure.

Temporal Blocking

Execute multiple time steps on a tile before exchanging halos:
- Standard approach: compute 1 step → exchange halos → compute 1 step → ...
- Temporal blocking: compute T steps locally (tile shrinks by halo width per step) → exchange once. Communication reduced by factor T.
- Overlapped tiling: extend tile by T×halo_width in each direction. Compute T steps, then trim the overlap region (which may be incorrect due to missing neighbor data). Interior results are correct. Trades redundant computation for reduced communication.

Performance Metrics

A 7-point 3D stencil on H100 GPU achieves ~2.5 TB/s effective bandwidth using FP32 — approaching the 3.35 TB/s HBM3 peak. On a 1000-GPU cluster with NVLink/IB interconnect, weak scaling efficiency of 85-95% is achievable for large domains.

Parallel Stencil Computation is the canonical example of memory-bandwidth-bound parallel computing — the regular, predictable access pattern that serves as the benchmark for memory system optimization and whose performance directly determines the time-to-solution for the fluid dynamics, weather, and geophysics simulations that model the physical world.