Case Study 2: A Convergence-Driven, Non-Blocking Distributed Solver

"A good parallel program does two things at once: it computes, and it hides the fact that it is also communicating."

Executive Summary

The Project Checkpoint gave you a distributed heat solver that marches a fixed number of steps and exchanges halos with the simple, blocking mpi_sendrecv. That is the right place to start and the wrong place to stop. Real solvers do not march a guessed number of steps; they run until the physics settles — until the plate reaches steady state — and they decide collectively, because a distributed simulation has no single process that can see the whole field. And real solvers do not sit idle inside a blocking exchange while the network works; they overlap the communication with computation that does not need it. This study builds both. You will design a global convergence test using mpi_allreduce — so every process stops on the same step, together — and then restructure the halo exchange to be non-blocking, splitting the update so the strip's deep interior is computed while its edge rows are still in flight. The result is a solver that is both self-terminating and latency-hiding: the shape of code that runs on real clusters. Where Case Study 1 repaired a broken solver, this one designs a better one, one tier up.

Skills applied: mpi_allreduce for a global reduction every process needs (§34.3); why all-reduce, not reduce, for a stopping decision (§34.3); non-blocking mpi_irecv/mpi_isend/mpi_waitall and computation/communication overlap (§34.4); the surface-to-volume argument for when overlap pays (§34.4, 💡); verifying determinism against process count (Project Checkpoint honesty note).

Background

We start from the Checkpoint solver: each process owns a strip u(nx, 0:nloc+1), exchanges halos, and applies the Chapter 24 stencil. Two design goals drive the redesign:

1. Stop at steady state, not at a guessed step count. The heat equation relaxes toward a steady field where nothing changes; the natural stopping rule is "halt when the largest temperature change anywhere on the plate, in one step, falls below a tolerance." The catch: anywhere on the plate spans every process, so no single process can evaluate the rule alone.
2. Hide the halo behind computation. Inside a blocking mpi_sendrecv, the process waits on the network. If instead it starts the exchange, computes the rows that do not depend on the incoming ghosts, and only then waits, the network latency disappears under useful work — at least when there is enough interior work to hide it.

Phase 1 — A Global Convergence Test with mpi_allreduce

Each process can compute its own largest cell change in a step, local_max = maxval(abs(u_new − u)) over its owned rows. The stopping rule needs the global maximum over all strips — and, crucially, every process must learn it, because every process runs the time loop and every process must decide, identically, whether to take another step. If only one process knew, the others would march on and the program would diverge into ranks doing different numbers of steps — a deadlock waiting to happen at the next collective.

This is precisely the case mpi_allreduce exists for: combine a value across all processes and give the result to all of them.

local_max = maxval(abs(u_new(:,1:nloc) - u(:,1:nloc)))          ! my biggest change this step
call mpi_allreduce(local_max, global_max, 1, MPI_DOUBLE_PRECISION, &
                   MPI_MAX, MPI_COMM_WORLD, ierr)               ! every rank now holds the global max
if (global_max < tol) exit                                     ! ...so every rank exits TOGETHER

Using mpi_reduce here would be a bug: only the root would hold global_max, so only the root would see the exit condition, and the others would loop forever waiting at the next mpi_allreduce for a root that has already left. The stopping decision must be unanimous, so the reduction must be an all-reduce. This is the distributed-memory face of a simple truth: a collective decision needs collective information.

Phase 2 — The Convergence Loop, and What It Prints

Here is the redesigned time loop (the full program is a direct extension of code/project-checkpoint.f90, adding the tolerance test and a step cap for safety):

do step = 1, max_steps
  call exchange_halos(u, nx, nloc, up, down)                   ! ghosts first (Case Study 1's lesson)
  call apply_stencil(u, u_new, nx, nloc, r)                    ! Chapter 24 update over owned rows
  local_max = maxval(abs(u_new(:,1:nloc) - u(:,1:nloc)))
  call mpi_allreduce(local_max, global_max, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ierr)
  u(:,1:nloc) = u_new(:,1:nloc)
  if (rank == 0) print '(a,i0,a,f6.1)', 'step ', step, '  global max change = ', global_max
  if (global_max < tol) then
    if (rank == 0) print '(a,i0,a)', 'converged at step ', step, ' (max change below tolerance)'
    exit
  end if
end do

Run it on the Checkpoint's $5 \times 6$ plate (top edge 100, $r = 0.2$), on 2 processes, with a tolerance of 1.0 and a generous step cap. The global max change is the same sequence you hand-computed for the Checkpoint — 20.0, then 12.0, then 6.4 — decreasing monotonically as the plate approaches steady state:

$ mpif90 -std=f2018 -Wall -O2 heat_converge.f90 -o hc && mpirun -np 2 ./hc
step 1  global max change =   20.0
step 2  global max change =   12.0
step 3  global max change =    6.4
   ... (continues; each step's change is the global MAX across both strips, via allreduce)
converged at step NN (max change below tolerance)

The first three numbers are exact and hand-checkable (they are the Checkpoint's); the step NN at which the change finally drops below 1.0 depends on running the full relaxation, which we do not tabulate here because the point is the mechanism, not the count. What matters is that the change is a single global number every process agrees on, computed by mpi_allreduce, so all processes stop on the same step. Change to -np 1 or -np 4 and the printed sequence of global max changes is identical — the reduction combines the same per-cell changes regardless of how the strips are divided.

Honesty note. We present the per-step max change (a computed quantity we can verify) but never a per-step time. How long each step takes, and how the two processes interleave, depends on the machine and is not something this study measured. The physics is deterministic; the schedule is not.

Phase 3 — Non-Blocking Halo Exchange

Now the latency-hiding redesign. The blocking exchange_halos makes the process wait on the network before it does any update. But observe which rows actually need the ghosts: only the edge rows (local rows 1 and nloc) read u(:,0) and u(:,nloc+1). The deep interior rows (2 through nloc-1) read only owned data — they can be updated the instant the step begins, without waiting for anything. So we split the work:

! post the receives first (landing spots ready), then the sends -- none of these block
call mpi_irecv(u(:,0),      nx, MPI_DOUBLE_PRECISION, up,   1, MPI_COMM_WORLD, reqs(1), ierr)
call mpi_irecv(u(:,nloc+1), nx, MPI_DOUBLE_PRECISION, down, 0, MPI_COMM_WORLD, reqs(2), ierr)
call mpi_isend(u(:,1),      nx, MPI_DOUBLE_PRECISION, up,   0, MPI_COMM_WORLD, reqs(3), ierr)
call mpi_isend(u(:,nloc),   nx, MPI_DOUBLE_PRECISION, down, 1, MPI_COMM_WORLD, reqs(4), ierr)

call update_rows(u, u_new, 2, nloc-1, r)     ! DEEP interior: needs no ghosts -- compute while halo flies

call mpi_waitall(4, reqs, stats, ierr)       ! now the ghosts have landed...
call update_rows(u, u_new, 1, 1,       r)    ! ...top edge row, which needed the upper ghost
call update_rows(u, u_new, nloc, nloc, r)    ! ...bottom edge row, which needed the lower ghost

The tags follow §34.4's direction convention (up-going messages tag 0, down-going tag 1), so each isend matches the neighbour's irecv. Between the isends and the mpi_waitall, the process computes its deep interior; on a fat strip that is most of the work, and it runs while the network moves the edge rows. When mpi_waitall returns, the ghosts are current and the two edge rows finish the step. The communication has been hidden behind computation.

⚠️ The rule you must not break: do not read or write the halo buffers u(:,0) / u(:,nloc+1), or the send buffers u(:,1) / u(:,nloc), between the i-calls and the mpi_waitall. The transfer is in flight; touching those rows early is a race. The deep-interior update is safe precisely because it touches none of them.

Phase 4 — When Does the Overlap Actually Pay?

Be honest about the payoff, because it is not free and it is not always worth it. The deep interior has nloc − 2 rows; the edges are always 2 rows. On the Checkpoint's tiny strips, nloc = 2, so the deep interior is empty — there is nothing to overlap, and the non-blocking version buys nothing but complexity. The overlap pays only when nloc is large enough that the deep-interior update takes longer than the halo transfer, so the exchange finishes "for free" underneath it. This is the surface-to-volume argument of §34.4 wearing work clothes: computation scales with the strip's area (nloc × nx), communication with its perimeter (2 × nx), so the fatter the strip, the more thoroughly its interior work hides its boundary exchange.

The design lesson is to keep both exchange routines and choose by scale: blocking mpi_sendrecv for its simplicity when strips are thin or the code is young, non-blocking overlap when the strips are fat and the run is large enough that hiding the halo matters. Reaching for the complex version too early is a premature optimization; never reaching for it at scale leaves performance on the table.

Phase 5 — Verify, Then Reason About Scaling

The correctness check is the same discipline as always: the convergence-driven, non-blocking solver must print the identical field and the identical global-max-change sequence as the simple Checkpoint solver, on any process count, because none of these changes touches the per-cell arithmetic. Run all three — Checkpoint (blocking, fixed steps), convergence (blocking, allreduce stop), and overlap (non-blocking) — on the $5 \times 6$ plate and confirm the step-1–3 max changes are 20.0, 12.0, 6.4 in every case. If they diverge, a communication bug has crept in, and the deterministic reference has caught it.

With correctness pinned, the redesign is what lets the solver scale honestly: it stops when the science says to (not when a guess says to), and it hides its communication when the strips are fat enough to do so. Both are prerequisites for the weak-scaling studies of the Chapter 38 capstone, where the plate grows with the process count and the per-step allreduce and halo exchange must stay a small fraction of the step. You have built the two mechanisms that keep them small.

Discussion Questions

1. Why must the convergence test use mpi_allreduce rather than mpi_reduce followed by a broadcast of the result? (Both would give every rank the value — so is it only about convenience, or is there a correctness difference if a rank misses the broadcast?)
2. The non-blocking version computes the deep interior between the isends and the waitall. Sketch what goes wrong if a programmer, trying to "do more work," also updates an edge row in that window.
3. The overlap buys nothing at nloc = 2 and a great deal at nloc = 1000. Given a fixed plate and a target process count, how would you decide at run time which exchange routine to call? What single quantity would you branch on?

Your Turn: Extensions

• Option A. Add a second reduced quantity to the convergence report: alongside the global max change, compute the global sum of the field (total heat) each step with a second mpi_allreduce using MPI_SUM, and watch it approach the steady-state total. What does a non-monotone total tell you about a bug?
• Option B. Implement both exchange_halos_blocking and exchange_halos_nonblocking behind one exchange_halos that dispatches on nloc (blocking below a threshold, non-blocking above). Justify your threshold in terms of the surface-to-volume table.
• Option C. Extend the convergence loop to also stop if global_max rises between steps (a sign of CFL instability — Chapter 24), printing a clear diagnostic. Why is a rising global max change a reliable distributed detector of an unstable timestep, and why does allreduce make it cheap to check?

Key Takeaways

• A distributed stopping decision needs a global reduction every process sees: compute each strip's local max change, combine with mpi_allreduce/MPI_MAX, and let every rank exit on the same step. mpi_reduce would strand the non-root ranks.
• Non-blocking exchange hides the halo behind computation: post irecv/isend, update the ghost-independent deep interior, mpi_waitall, then finish the edge rows — but never touch the in-flight buffers before the wait.
• Overlap pays by surface-to-volume: it helps only when the strip is fat enough that its interior work outlasts its boundary transfer. Keep the simple blocking version for thin strips.
• Determinism is the check: every redesign must reproduce the reference field and max-change sequence on any process count, because none of it alters the per-cell arithmetic.