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.