Case Study 2: A Parallelization Plan for the Heat Solver

"Weeks of coding can save you hours of planning." — a programmer's rueful proverb, inverted on purpose

Executive Summary

You are about to spend the next four chapters making the heat solver parallel — coarrays, OpenMP, MPI, GPU. Before you write one directive, you will do what a professional does: build a model of the payoff and a plan for the work, on paper and in a few lines of Fortran, so that every later chapter has a target to hit and a way to know whether it hit it. This study turns the abstract method of §31.5 into a concrete artifact: a profile-driven serial-fraction estimate, an explicit map of what is data-parallel and what is dependency-bound, a small program that predicts both the strong-scaling (Amdahl) and weak-scaling (Gustafson) outcomes across core counts, and a reasoned choice of memory model. The deliverable is a one-page plan you could hand to a colleague — the thing that separates engineering from hopeful directive- sprinkling. Where Case Study 1 analyzed someone else's finished measurements, this one designs the experiment before it exists.

Skills applied: estimating a serial fraction from a profile (§31.5, and Chapter 28); Amdahl and Gustafson modelling (§31.2); data vs task parallelism and the time-loop dependency (§31.4, §31.5); the shared/distributed/GPU taxonomy and tool choice (§31.3); building a compilable planning tool.

Background

By the end of Part VII the solver is a capable serial program: a heat_solver module with a laplacian and a step, a field_t derived type, namelist configuration, VTK output, timers, all built with -O3 -march=native -flto. It runs one time-stepping loop; inside each step it sweeps the five-point stencil over the interior of the plate. It is fast — and it uses exactly one core of whatever machine it runs on. The task now is to plan its escape onto many cores without guessing. We follow the five steps of §31.5, and we produce, at Step 4, a reusable modelling tool.

Phase 1 — From Profile to Serial Fraction

Step 1 (make it fast serially) and Step 2 (profile) are already done — that was Part VII. We take the profile's verdict as our starting data. Suppose a representative run reports the following breakdown of a $1000 \times 1000$, 10,000-step simulation (round numbers, for the model — your real profile will differ):

The parallel fraction is the stencil work: $p = 196/200 = 0.98$. The serial fraction — setup plus I/O — is $1 - p = 4/200 = 0.02$. That single number, 0.98, will drive the entire model. It is worth pausing on how small the serial part is and how much it still matters: 2% of the run is going to cap the whole parallelization, as we are about to see.

Phase 2 — Map the Parallelism and the Dependencies

Step 3 is structural, not numerical, and it is where most of the design lives. Two questions: what is data-parallel, and what dependency forbids parallelism?

• The stencil sweep is data-parallel. Within one time step, every interior cell's new temperature is computed from the old values of its four neighbours. No interior update reads another interior update's new value, so all of them are independent and could, in principle, happen at once — the purest data parallelism (§31.4). This is the 98%, and it is what every parallel model in Part VIII will target. That it reads only old neighbours is the two-array ("current" and "next") structure from Chapter 24; if the update read new neighbours mid-sweep, the cells would depend on each other and the easy parallelism would vanish.
• The time loop is sequential. Step $n+1$ needs the finished field of step $n$ — you cannot compute the future before the present. This dependency is inherent and cannot be parallelized away; you parallelize the work inside each step and let the steps run in order. Anyone who proposes running consecutive time steps concurrently has misread the algorithm.

There is essentially no task parallelism here worth chasing — the solver does one kind of work, not several independent kinds — which is entirely normal for a numerical simulation and, in fact, good news: data parallelism is the kind that scales.

Phase 3 — Build the Payoff Model

Step 4 is the artifact. Rather than estimate in our heads, we write a small program that, from the profiled $p = 0.98$, prints the strong-scaling (Amdahl, fixed problem) speedup and efficiency and the weak-scaling (Gustafson, growing problem) speedup at each core count. Having both in one table forces the strong-vs-weak decision into the open:

program solver_plan_model
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  real(dp), parameter :: p = 0.98_dp          ! parallel fraction from the profile
  integer,  parameter :: cores(7) = [1, 2, 4, 8, 16, 32, 64]
  integer  :: k, n
  real(dp) :: strong, weak, eff

  print '(a, f4.2)', 'Solver parallelization model, parallel fraction p = ', p
  print '(a)', ' cores      strong  efficiency        weak'
  do k = 1, size(cores)
    n      = cores(k)
    strong = 1.0_dp / ((1.0_dp - p) + p / real(n, dp))   ! Amdahl: fixed problem
    eff    = strong / real(n, dp)
    weak   = (1.0_dp - p) + p * real(n, dp)              ! Gustafson: grow problem
    print '(i6, 3f12.4)', n, strong, eff, weak
  end do
end program solver_plan_model

$ gfortran -std=f2018 -Wall solver_plan_model.f90 -o plan && ./plan
Solver parallelization model, parallel fraction p = 0.98
 cores      strong  efficiency        weak
     1      1.0000      1.0000      1.0000
     2      1.9608      0.9804      1.9800
     4      3.7736      0.9434      3.9400
     8      7.0175      0.8772      7.8600
    16     12.3077      0.7692     15.7000
    32     19.7531      0.6173     31.3800
    64     28.3186      0.4425     62.7400

Read the table as a decision aid. The strong column (this exact plate, solved faster) is capped: it will never beat $1/(1-p) = 50\times$, it delivers a healthy 7× at 8 cores (88% efficiency) and 12× at 16 (77%), and then efficiency slides to 44% by 64 cores — the classic Amdahl bend. The weak column (a bigger, finer plate in the same wall-clock time) climbs almost linearly to 62.7× at 64 cores, because growing the problem keeps every core busy and shrinks the serial fraction's relative bite. The plan's first decision writes itself from these two columns.

Phase 4 — Choose the Model and the Scale

Step 5 turns the model into commitments. Two decisions:

Strong or weak? If the scientific need is this plate solved faster — a fixed grid, a deadline — you are in the strong column: target 8–16 shared-memory cores, where efficiency is still 77–88%, and do not expect miracles past that on a fixed problem. If the need is resolution — a finer grid, more physics, better accuracy — you are in the weak column, and you can profitably ride many more cores by growing the plate to match them. Most research use of the solver is the second kind, which argues for a design that can scale the grid, not just the core count.

Which memory model? Now the taxonomy of §31.3 pays off, matched to the target:

A sensible plan for most readers: start with OpenMP (Chapter 33) for a quick, high-efficiency 8–12× on one node; then, if the problem outgrows a node, move to coarrays or MPI (Chapters 32, 34) with the plate decomposed into tiles and halos exchanged at the seams; consider the GPU (Chapter 35) if the stencil dominates and the host-device transfer can be tamed. Each is the same 98%-parallel hot spot, approached through a different door.

Phase 5 — Sanity Check and the One-Page Plan

Before trusting the model, sanity-check its anchors by hand. The ceiling: $1/(1 - 0.98) = 1/0.02 = 50$ — matches the strong column's approach. One row: at 8 cores, strong $= 1/(0.02 + 0.98/8) = 1/0.1425 = 7.02$, and weak $= 0.02 + 0.98 \times 8 = 7.86$ — both match the printout. The model is arithmetic we can verify, not a black box. The finished plan reads:

Heat-solver parallelization plan (v1). Profile: $p = 0.98$ (stencil), serial 2% (setup + VTK I/O). Hot spot: the interior stencil update, data-parallel within each step. Hard dependency: the time loop is sequential (step $n+1$ needs step $n$). Strong-scaling ceiling: 50×; efficient range 8–16 cores. First target: OpenMP on the update loop (Chapter 33), expect 7–12× at 88–77% efficiency; measure against this model and read the gap as overhead. Next: coarray/MPI decomposition for beyond-one-node scale (Chapters 32, 34); grow the grid for weak scaling when resolution, not deadline, is the goal. Serial-fraction risk: the 2% I/O caps strong scaling at 50× — revisit with MPI-IO (Chapter 34) if needed.

That is the deliverable: a page that names the target, the tool, the expected number, and the way to check it — written before a single line of parallel code exists.

Discussion Questions

1. The model says weak scaling reaches 62.7× at 64 cores while strong scaling stalls near 28×. Give a concrete scientific scenario for the heat plate in which the strong number is the one you must have, despite being smaller.
2. The plan starts with OpenMP rather than MPI even though MPI scales further. Justify that ordering in terms of effort, the profiled ceiling, and the likely problem size — and say what would change your mind.
3. The 2% serial fraction is entirely setup and I/O. Sketch one change to the solver's output strategy that would lower it, and estimate (with the ceiling formula) the new strong-scaling ceiling if you halved the serial fraction to 1%.

Your Turn: Extensions

• Option A. Change p in the model to 0.90 and to 0.999 and re-run. Tabulate how the strong ceiling and the 64-core efficiency move. Which single input has the largest effect on the plan, and why does that argue for attacking the serial fraction first?
• Option B. Add a fourth column to the model that prints the strong-scaling efficiency as a percentage and flags the largest core count whose efficiency is still above 70% — an automatic "recommended core count" for a fixed problem. What does it recommend at $p = 0.98$? (Hint: check where efficiency crosses 0.70 in the printout.)
• Option C. Extend the model to take a communication-overhead term, $T(N) = (1-p) + p/N + c\,N$ for a small constant $c$, and find the core count that minimizes the modelled time. Relate the result to Code B of Case Study 1 — you are now predicting the sweet spot you there diagnosed.

Key Takeaways

• Model the payoff before you build it. A profiled serial fraction plus Amdahl and Gustafson gives you the ceiling, the efficient range, and the strong-vs-weak decision — in a few lines of verifiable Fortran.
• Map the parallelism first. The solver's stencil update is data-parallel (the 98%); its time loop is a hard sequential dependency (parallelize within a step, never across steps).
• Let the target choose the tool. One node → OpenMP/coarrays; a cluster → MPI/coarrays; a GPU → OpenACC — all attacking the same hot spot through the §31.3 taxonomy.
• The output of this chapter is not code but a plan with a number in it: what you will parallelize, what it should buy, and how you will know if it did. Every chapter of Part VIII is measured against it.