Case Study 2: A Stencil in Nine Lines

DataField.Dev

Case Study 2: A Stencil in Nine Lines

"Every neighbor computation is a shifted copy of the same grid. Once you see that, the loops disappear."

Executive Summary

Where the first case study ported existing code, this one asks you to build something from scratch — and to build it the array way from the first keystroke. The target is Conway's Game of Life: a grid of cells, each alive or dead, that evolves by a rule depending on how many of its eight neighbors are alive. The naive implementation is a quadruple-nested loop (two loops over cells, two over the neighborhood). We will instead compute the entire generation with array sections — eight shifted copies of the grid summed into a neighbor count, then one masked assignment for the rule — and verify it on the classic "blinker." The payoff is not just elegance: this is a stencil computation, structurally identical to the heat-solver Laplacian you built in the Project Checkpoint, so the technique you practice here is the technique that will carry your simulation to the capstone.

Skills applied: representing a 2D field as an array (§5.1–5.2); the shifted-section idiom for stencils (§5.2); summing sections as a whole-array operation (§5.3); masked assignment with where and compound logical masks (§5.7); the connection between neighbor stencils and column-major performance (§5.6).

Background

The Game of Life runs on an infinite grid of cells, each alive (we store 1) or dead (0). Time advances in generations; every cell's next state depends only on its current state and its eight immediate neighbors, by three rules:

A live cell with two or three live neighbors survives; otherwise it dies (loneliness or overcrowding).
A dead cell with exactly three live neighbors becomes alive (birth).
Every other cell is dead next generation.

We approximate the infinite grid with a finite one and a dead border: we compute only the interior (2:n-1, 2:n-1) and let the outer ring stay dead, which is exact as long as the pattern does not reach the edge. Our test pattern is the blinker — three live cells in a row — which is famous for oscillating between horizontal and vertical with period two. If our code turns a horizontal blinker into a vertical one in a single step, it is right.

Phase 1 — Represent the Board

The board is the most natural array there is: a rank-2 integer grid, 1 for alive and 0 for dead. We need three like it — the current grid, the neighbor counts, and the next grid — all the same shape:

integer, parameter :: n = 5
integer :: grid(n, n), neigh(n, n), next(n, n)

grid = 0;  neigh = 0;  next = 0
grid(3, 2) = 1;  grid(3, 3) = 1;  grid(3, 4) = 1     ! a horizontal blinker

Storing alive/dead as integer 1/0 (rather than logical) is a deliberate choice: it lets us add neighbors arithmetically. Counting is just summing, and summing is what arrays do best.

Phase 2 — Count Neighbors with Eight Shifted Sections

Here is the idea that dissolves the loops. The count of live neighbors of every interior cell is the sum, over the eight compass directions, of the grid shifted one step in that direction. Each shift is an array section, and each section lines up so that its (a,b) entry is one particular neighbor of interior cell (a+1, b+1). Sum the eight sections and every interior cell receives its neighbor total at once:

neigh(2:n-1, 2:n-1) =                                              &
      grid(1:n-2, 1:n-2) + grid(1:n-2, 2:n-1) + grid(1:n-2, 3:n)   &   ! up-left,  up,   up-right
    + grid(2:n-1, 1:n-2)                      + grid(2:n-1, 3:n)   &   ! left,           right
    + grid(3:n,   1:n-2) + grid(3:n,   2:n-1) + grid(3:n,   3:n)       ! down-left, down, down-right

Read the three "row" groups as the row above, the same row, and the row below, and the three columns in each as left, center, right. The center of the middle group is deliberately absent — a cell is not its own neighbor. Every section is (n-2) × (n-2), so they are conformable and the sum is a single whole-array operation. This is the exact same shifted-section trick as the heat Laplacian; there we summed four neighbors and subtracted the center, here we sum eight and keep the center out.

Phase 3 — Apply the Rule with One Masked Assignment

With neighbor counts in hand, the birth-and-survival rule is a single logical mask. A cell is alive next generation if it is alive now with 2 or 3 neighbors, or dead now with exactly 3. Express that compound condition as the mask of a where, assigning 1 to the interior of next where it holds (the rest of next was zeroed already):

where ( (grid(2:n-1,2:n-1) == 1 .and.                             &
         (neigh(2:n-1,2:n-1) == 2 .or. neigh(2:n-1,2:n-1) == 3))  &   ! survives
   .or. (grid(2:n-1,2:n-1) == 0 .and. neigh(2:n-1,2:n-1) == 3) )      ! born
  next(2:n-1, 2:n-1) = 1
end where

No loop over cells, no if per cell. The rule is written once, as an array condition, and applied to the whole interior simultaneously. The complete program:

program game_of_life
  implicit none
  integer, parameter :: n = 5
  integer :: grid(n, n), neigh(n, n), next(n, n)
  integer :: i

  grid = 0;  neigh = 0;  next = 0
  grid(3, 2) = 1;  grid(3, 3) = 1;  grid(3, 4) = 1     ! horizontal blinker

  neigh(2:n-1, 2:n-1) =                                              &
        grid(1:n-2, 1:n-2) + grid(1:n-2, 2:n-1) + grid(1:n-2, 3:n)   &
      + grid(2:n-1, 1:n-2)                      + grid(2:n-1, 3:n)   &
      + grid(3:n,   1:n-2) + grid(3:n,   2:n-1) + grid(3:n,   3:n)

  where ( (grid(2:n-1,2:n-1) == 1 .and.                             &
           (neigh(2:n-1,2:n-1) == 2 .or. neigh(2:n-1,2:n-1) == 3))  &
     .or. (grid(2:n-1,2:n-1) == 0 .and. neigh(2:n-1,2:n-1) == 3) )
    next(2:n-1, 2:n-1) = 1
  end where

  print '(a)', 'generation 0 (1 = alive):'
  do i = 1, n
    print '(5i2)', grid(i, :)
  end do
  print '(a)', 'generation 1:'
  do i = 1, n
    print '(5i2)', next(i, :)
  end do
end program game_of_life

$ gfortran -std=f2018 -Wall game_of_life.f90 -o life && ./life
generation 0 (1 = alive):
 0 0 0 0 0
 0 0 0 0 0
 0 1 1 1 0
 0 0 0 0 0
 0 0 0 0 0
generation 1:
 0 0 0 0 0
 0 0 1 0 0
 0 0 1 0 0
 0 0 1 0 0
 0 0 0 0 0

Phase 4 — Verify on the Blinker

The horizontal bar became a vertical bar — the blinker's signature. Confirm it cell by cell against the rules, at the three cells that change:

Cell (2,3) is dead, and its neighbors include the whole bar (3,2),(3,3),(3,4) — exactly 3 live neighbors — so it is born. Likewise (4,3) sees the same three and is born. Those are the new top and bottom of the vertical bar.
Cell (3,3), the bar's center, is alive with neighbors (3,2) and (3,4) — 2 live — so it survives.
The bar's ends, (3,2) and (3,4), are alive but each has only 1 live neighbor (the center), so they die of loneliness.

Three born-or-surviving cells at (2,3),(3,3),(4,3); the two ends gone: a vertical blinker, exactly as printed. (Every other interior cell works out to fewer than three neighbors and stays dead — check (2,2), which sees (3,2) and (3,3), only 2, and is not born.)

Phase 5 — Performance, and the Bridge to Your Solver

Count what we did: eight section reads, one array sum, one masked assignment — a fixed handful of whole-array operations, regardless of grid size. The naive version would be four nested loops with an if in the innermost, executed n² times. Ours does the same arithmetic, but expressed so the compiler sees whole arrays of known shape and can vectorize the neighbor sum. And because every section is read in column-major-friendly order (§5.6) — each is a contiguous-ish sweep with the memory grain — the access pattern is exactly the one Fortran is fast at. This is the general lesson of stencil computing: a neighbor update is a sum of shifted copies of the field, and shifted copies are array sections.

That is not a coincidence you will use once. Your heat solver's Laplacian is the same pattern with four neighbors instead of eight; a blur filter, an edge detector, a wave equation, a fluid solver — all are stencils, all are shifted-section sums. The technique you just practiced on a toy is the technique Chapter 24 formalizes for the real heat equation and Chapter 29 tunes for speed. You built a cellular automaton; you also built a rehearsal for the capstone.

Discussion Questions

Our dead-border trick is exact only while the pattern avoids the edge. Describe two other boundary policies (for example, a wrap-around periodic grid), and sketch how each would change the section expressions. (Periodic boundaries return in Chapter 24.)
We used integer 1/0 so we could add neighbors. What would break, or become awkward, if the board were logical? What intrinsic could still count .true. values?
To run many generations, you must feed next back in as grid. Why is it a bug to update grid in place, cell by cell, rather than computing a whole next and swapping? Relate this to the no-aliasing idea from Chapter 1.

Your Turn: Extensions

Option A. Wrap the update in a do gen = 1, 4 loop that swaps next into grid each generation (grid = next then re-zero next), and confirm the blinker returns to horizontal after two steps. Print each generation.
Option B. Replace the blinker with a glider (a five-cell pattern that walks diagonally) on a larger grid, and watch it move. Does your dead border eventually corrupt it? At what generation, and why?
Option C. Swap the rule to compute the heat-style four-neighbor Laplacian instead of the eight-neighbor count, on a real(dp) field, reusing the exact same shifted-section skeleton. You have now converted a Game-of-Life engine into a diffusion kernel by changing two lines — the clearest possible evidence that a stencil is a stencil.

Key Takeaways

A neighbor/stencil update is a sum of shifted array sections; seeing that turns nested loops into a handful of whole-array statements.
Storing states as integer 1/0 lets you count by summing, the operation arrays are built for.
A compound rule becomes a single where with a logical mask; no per-cell if, no per-cell loop.
The Game of Life, the heat Laplacian, blur filters, and wave solvers are the same computational shape — master the shifted-section idiom once and you have the core of every structured-grid simulation, straight through to the capstone.