Case Study 2: Building a Distributed Heat Rod

"A parallel program is correct when it gives the serial answer — and useful when it decides, together, that it is done."

Executive Summary

The Project Checkpoint distributed the 2D plate; here you build, from scratch, the one-dimensional sibling — a heated rod split across images — and take it two steps further than the checkpoint did. First you make the images agree on a global stopping criterion: instead of a fixed number of steps, they iterate until the largest temperature change anywhere on the rod drops below a tolerance, which requires combining a per-image quantity into a global one with the co_max collective. Second you analyze the design: you prove the parallel answer equals the serial one, you work out the communication-to-computation ratio and what it says about scaling, and you identify the one optimization — overlapping communication with computation — that the next chapter will make real. The result is a small but complete piece of distributed-memory science: decomposition, halo exchange, a collective reduction for convergence, a validation, and a scaling argument, all in one program you can hand-trace.

Skills applied: SPMD decomposition and per-image ownership (§32.1); allocatable coarrays and coindexed halo reads (§32.2); segment ordering and the two-barrier step (§32.3); the co_max collective for a global reduction (§32.4); building and validating with -fcoarray=single then caf/cafrun (§32.5); the communication-cost reasoning that connects to strong/weak scaling (Chapter 31).

Background

The 1D heat equation on a rod is the plate's stencil with the two $y$-neighbours dropped (Chapter 24): $u_i^{n+1} = u_i^n + r\,(u_{i-1}^n - 2u_i^n + u_{i+1}^n)$, stable for $r = \alpha\Delta t/\Delta x^2 \le \tfrac12$. Take a rod of six nodes, the two ends held Dirichlet — the left cold at $0$, the right hot at $100$ — the four interior nodes starting cold, and $r = 0.25$. Marched serially, it warms from the hot end leftward, and after three steps the interior reads $0,\ 0,\ 1.5625,\ 12.5,\ 45.3125$ (with the ends $0$ and $100$), eventually settling to the linear steady state $0, 20, 40, 60, 80, 100$. Those exact dyadic fractions are our validation target: the distributed version must reproduce them to the bit.

Phase 1 — Design: One Rod, Many Images

The decomposition mirrors the plate's, one dimension simpler. The four interior nodes are split into equal contiguous segments, one per image. Each image's local array carries its owned nodes plus one halo node on each side, holding a copy of its neighbour's edge node:

   global rod:   [1]  2   3   4   5  [6]          ends [1]=0 (cold), [6]=100 (hot)
                  |   \___/   \___/   |
   2 images:     img1 owns 2,3     img2 owns 4,5
   img1 local:  ( 1 )  2   3  ( 4 )   local 1 = left halo (=node 1, cold end)
                              \__ local 4 = right halo  <-  img2's node 4
   img2 local:  ( 1 )  4   5  ( 6 )   local 4 = right halo (=node 6, hot end)
                \__ local 1 = left halo  <-  img1's node 3

Owned nodes are local indices 2 : nloc+1; halos are local 1 and nloc+2. The physical ends are simply the halos of the end images, set once and never exchanged (image 1's left halo stays $0$; the last image's right halo stays $100$). The interior halos — image 1's right, image 2's left — are refreshed every step by coindexed reads. The invariant that makes it all work: the coarray's local shape is identical on every image (nloc+2 nodes), as coarrays require, so u(2)[me+1] means the same slot on every image.

Phase 2 — Build: Halo, Stencil, and a Collective Convergence Test

The new idea beyond the checkpoint is the global stopping test. Each image can compute the largest change over its own nodes this step, but "are we converged?" is a question about the whole rod — so we combine the per-image maxima into one global maximum with co_max, and every image learns the same residual. Here is the complete program:

program rod_caf
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  integer,  parameter :: mglob = 6, nint = mglob - 2   ! 6 nodes, 4 interior
  real(dp), parameter :: r = 0.25_dp                   ! <= 1/2 in 1D: stable (Ch.24)
  real(dp), allocatable :: u(:)[:]                      ! local segment + 2 halos (COARRAY)
  real(dp), allocatable :: u_new(:)
  integer  :: me, ni, nloc, it, q, lo
  real(dp) :: change
  real(dp), allocatable :: rod(:)                       ! image 1 assembles the full rod

  me = this_image();  ni = num_images()
  if (mod(nint, ni) /= 0) then
    if (me == 1) print '(a)', 'error: interior nodes must divide evenly across images'
    error stop 1
  end if
  nloc = nint / ni
  allocate(u(nloc+2)[*]);  allocate(u_new(nloc+2))
  u = 0.0_dp
  if (me == ni) u(nloc+2) = 100.0_dp                    ! hot right end (last image's right halo)

  do it = 1, 3                                          ! real code: do while (change > tol)
    sync all
    if (me > 1)  u(1)      = u(nloc+1)[me-1]            ! left halo  <- left neighbour's last node
    if (me < ni) u(nloc+2) = u(2)[me+1]                 ! right halo <- right neighbour's first node
    sync all
    u_new(2:nloc+1) = u(2:nloc+1) + r*( u(1:nloc) - 2.0_dp*u(2:nloc+1) + u(3:nloc+2) )
    change = maxval(abs(u_new(2:nloc+1) - u(2:nloc+1)))  ! this image's largest change
    call co_max(change)                                 ! GLOBAL residual: max over ALL images
    u(2:nloc+1) = u_new(2:nloc+1)                        ! commit
    if (me == 1) print '(a,i0,a,f10.4)', 'step ', it, '  global residual = ', change
  end do

  ! gather the whole rod onto image 1 and print it
  sync all
  if (me == 1) then
    allocate(rod(mglob));  rod = 0.0_dp;  rod(mglob) = 100.0_dp
    do q = 1, ni
      lo = (q-1)*nloc + 2
      rod(lo:lo+nloc-1) = u(2:nloc+1)[q]                ! coindexed read of each segment
    end do
    print '(a, 6f9.4)', 'rod after 3 steps: ', rod
  end if
end program rod_caf

$ gfortran -fcoarray=single -std=f2018 -Wall rod_caf.f90 -o rod && ./rod
step 1  global residual =    25.0000
step 2  global residual =    12.5000
step 3  global residual =     7.8125
rod after 3 steps:    0.0000   0.0000   1.5625  12.5000  45.3125 100.0000

The identical output appears under caf -std=f2018 -Wall rod_caf.f90 -o rod && cafrun -n 2 ./rod — same residuals, same rod — which is the point of Phase 3.

Phase 3 — Validate: Parallel Equals Serial

Correctness first, always. Trace the residuals by hand on two images. Image 1 owns nodes 2–3; image 2 owns 4–5. At step 1 only node 5 sees the hot end, changing $0 \to 0.25 \times 100 = 25$; every other interior node is still surrounded by zeros, so image 1's local maximum change is $0$ and image 2's is $25$. co_max combines them: global residual $= \max(0, 25) = 25$. At step 2, node 4 changes to $6.25$ and node 5 to $37.5$ (a change of $12.5$), while image 1's nodes barely move; global residual $= 12.5$. At step 3, node 3 finally warms to $1.5625$ because image 1's right halo now holds node 4's value $6.25$ from image 2 — the halo doing its job across the strip boundary — and node 5's change is $7.8125$, the new global maximum.

The assembled rod, $0,\ 0,\ 1.5625,\ 12.5,\ 45.3125,\ 100$, is exactly the serial result. Two images, halo exchange, and a collective convergence test reproduce the single-core answer to the bit — and would on 1, 2, or 4 images, since $4$ interior nodes divide all three. That image-count independence is the proof of correctness, the same acceptance test as Case Study 1: the physics cannot depend on how we cut the rod.

Phase 4 — Analyze: Communication Cost and the One Optimization

Now the engineering. Two questions decide whether a decomposition is worth it: how much does each image communicate relative to how much it computes, and can the communication be hidden.

Communication-to-computation ratio. Per step, an interior image of the rod computes $\sim M/P$ node updates (its segment) and communicates exactly $2$ halo nodes — a ratio of $2P/M$, which shrinks as the rod $M$ grows for fixed $P$. That is excellent: a long rod is almost all computation, a sliver of communication. The plate is the more instructive case. A strip image computes $\sim N^2/P$ cell updates and communicates $\sim 2N$ halo cells, a ratio $\sim 2P/N$ — still shrinking as the grid $N$ grows. This is the general and important truth of halo methods: the work is a volume ($N^2$) and the communication is a surface ($N$), so refining the grid improves the ratio. It is precisely why structured-grid codes weak-scale so well (Chapter 31 §31.4): grow the problem with the cores and each core stays compute-bound, its halo a rounding error.

The one optimization: overlap. Look at the two sync all barriers per step. Between them the images do nothing but communicate — the compute waits for the halo. But most of each segment's update (every node except the two touching a halo) needs only local data and could run while the halo is in flight. The optimization is to split the update: fire off the halo exchange, immediately compute the segment's interior interior (the nodes not adjacent to a halo) using data already in hand, and only then, once the halo has landed, compute the two edge nodes. Done well, the communication cost disappears behind computation that had to happen anyway. Coarrays with sync all express the simple version; the overlapped version wants one-sided or non-blocking communication, which is exactly the mpi_isend/mpi_irecv machinery of Chapter 34. We name the optimization here and cash it there.

Phase 5 — The Result

Read it as an engineer would write it up:

Distributed heat rod (v1). Model: 1D FTCS, $r = 0.25$, 6 nodes, ends $0$/$100$ Dirichlet. Decomposition: interior nodes split into equal contiguous segments, one per image; local shape identical across images; two halo nodes per image, exchanged each step by coindexed reads bracketed by sync all. Convergence: per-image maxval of the change, combined by co_max into a global residual; iterate until below tolerance. Validation: reproduces the serial rod ($0,0,1.5625,12.5,45.3125,100$ after 3 steps) exactly, on any image count dividing the 4 interior nodes; residuals $25 \to 12.5 \to 7.8125$. Cost: communication/computation $\sim 2P/M$, shrinking with problem size — weak-scales well. Next optimization: overlap the halo exchange with the interior update (Chapter 34, non-blocking).

That is a complete distributed solver in one page: decomposition, communication, a collective decision, a proof of correctness, and a cost model. The 2D plate of the Project Checkpoint is the same recipe with one more dimension of halo, and the capstone runs it at scale.

Discussion Questions

1. The convergence test uses co_max so every image learns the same global residual and they all stop together. What would go wrong if each image instead decided to stop based only on its own segment's maximum change?
2. Phase 4 argues the communication/computation ratio $\sim 2P/N$ shrinks as the grid grows. Explain why this makes weak scaling (grow the problem with the cores) far easier than strong scaling (fixed problem, more cores) for this solver, connecting to Amdahl and Gustafson from Chapter 31.
3. The overlap optimization computes the "interior interior" while the halo is in transit. What is the risk if you get the split wrong and accidentally use a halo value before it has arrived — and which chapter's tools (and which property of the update) let you do the overlap safely?

Your Turn: Extensions

• Option A. Replace the fixed do it = 1, 3 with do while (change > tol) (tolerance $10^{-3}$) and a guard against runaway iteration, so the rod runs to steady state. Predict (no code) the steady-state rod and confirm it is the straight line $0, 20, 40, 60, 80, 100$. Roughly how does the residual behave as it converges?
• Option B. Extend the rod to $M = 10$ nodes (8 interior) and run on 1, 2, and 4 images. Confirm the assembled rod is identical across all three, and that the first-step global residual is the same regardless of image count. Which line of code enforces the "divides evenly" precondition, and how would you relax it?
• Option C. Add a second collective: after convergence, use co_sum to compute the rod's total thermal energy $\sum_i u_i$ (a conserved-ish diagnostic) across all images, and print it from image 1. Verify by hand against the assembled rod.

Key Takeaways

• A distributed solver is decomposition + halo + a collective decision. Split the domain, exchange edges each step with coindexed reads and two barriers, and combine per-image quantities (a residual, an energy) with a collective like co_max/co_sum.
• Correctness is proven by reproducing the serial answer — bit-for-bit, and independent of the image count. If the parallel rod does not equal the serial rod, the halo or the synchronization is wrong.
• Halo methods weak-scale because work is volume and communication is surface. The ratio $\sim 2P/N$ shrinks as the grid grows, so bigger problems are more efficient — the deep reason structured-grid codes dominate HPC.
• The next win is overlap. Most of each update needs only local data and can run while the halo is in flight; hiding communication behind computation is the optimization Chapter 34's non-blocking messages make real.