parallel stencil computation
**Parallel Stencil Computation** is the **numerical method where each grid point is updated based on a fixed pattern of neighboring values (the stencil) — ubiquitous in computational fluid dynamics, weather simulation, image processing, and PDE solvers — and one of the most important parallel computing patterns because the regular, local data access pattern enables highly efficient parallelization through domain decomposition with halo exchange, achieving near-linear scaling to millions of cores when communication is properly overlapped with computation**.
**Stencil Pattern**
A 2D 5-point stencil:
```
new[i][j] = w0*old[i][j] + w1*old[i-1][j] + w2*old[i+1][j]
+ w3*old[i][j-1] + w4*old[i][j+1]
```
Each point depends only on its immediate neighbors. Applied to every point in a 2D/3D grid for each timestep. Examples: Jacobi iteration, Gauss-Seidel (with dependency ordering), heat equation, wave equation, weather prediction.
**Domain Decomposition**
The grid is divided into subdomains, one per processor. Each processor updates its local subdomain independently — except at subdomain boundaries, where stencil calculations need values from adjacent processors' domains.
**Halo Exchange (Ghost Cells)**
- **Ghost/Halo Region**: Each subdomain is padded with an extra layer of cells (1-3 layers depending on stencil radius) copied from neighboring processors.
- **Exchange Protocol**: Before each timestep, each processor sends its boundary cells to neighbors and receives neighbors' boundary cells into its ghost region. For a 2D decomposition with 4 neighbors, 4 send/receive pairs per timestep.
- **Communication Volume**: For an N×N local subdomain with a 1-cell halo, communication per timestep = 4N elements (surface) while computation = N² elements (volume). The surface-to-volume ratio decreases as N increases → larger subdomains have better computation-to-communication ratio.
**Optimization Techniques**
- **Communication-Computation Overlap**: Start halo exchange (non-blocking MPI_Isend/Irecv), compute interior points (which don't need ghost cells), then wait for halo exchange completion and compute boundary points. Hides communication latency behind useful computation.
- **Temporal Blocking (Tiling)**: Instead of exchanging halos every timestep, expand the halo by k cells and compute k timesteps before exchanging. Reduces communication frequency by k× at the cost of computing redundant cells in the expanded halo.
- **Cache-Oblivious Tiling**: Tile both spatial and temporal dimensions to maximize data reuse within the cache hierarchy. Achieved through recursive decomposition (space-time wavefront tiling).
- **Vectorization (SIMD)**: Stencil operations on contiguous grid rows vectorize naturally — adjacent grid points are processed by adjacent SIMD lanes. Array padding to cache-line boundaries maximizes vectorization efficiency.
**GPU Stencil Implementation**
Load a tile of the grid (plus halo) into shared memory. Each thread computes one grid point using shared memory reads (fast, no bank conflicts for stencil patterns). Thread blocks process tiles; the grid is tiled across the entire GPU grid of blocks.
Parallel Stencil Computation is **the poster child of structured parallel computing** — combining regular data access, predictable communication, and natural domain decomposition into a pattern that scales to the largest supercomputers on Earth, underpinning the simulations that predict weather, design aircraft, and model physical phenomena.