Chapter 38: Capstone — Your Complete Parallel Scientific Simulation, From Physics to Publication

"Essentially, all models are wrong, but some are useful." — George E. P. Box, statistician

Overview

Thirty-seven chapters ago you made a decision, not a program: a square plate, hot on one edge, cold on the others, and a promise to simulate how the heat spreads. Since then, every chapter has handed you one more piece. You made the temperature a 2D array, lifted the update into a step procedure, split the code into modules, bundled its state into a field_t, and then — in Chapter 24 — replaced the placeholder physics with a real five-point stencil, an explicit scheme, and the CFL condition that keeps it from exploding. You taught it to write VTK so you could watch it; you made it fast; you made it run across threads and across a cluster; you organized it into a real package and tested it.

This is the chapter where all of that becomes one thing. Not a longer program — a finished one: a modular, validated, optimized, parallel, visualized scientific simulation, and, just as important, a document that lets someone else believe your results. Because in computational science the code is only half the deliverable. The other half is the argument that the code is correct — that the numbers it prints are a faithful solution of the equation you meant to solve, and not an artifact of a bug, an unstable timestep, or a race condition. That argument has a standard form, and it is the form of a paper.

So we are going to do what a working computational scientist does at the end of a project: assemble the solver, verify it against a case with a known answer, measure how it scales, produce the figures, and write it up. By the end you will not just have a program that runs. You will have a result — one you could defend to a skeptical reviewer, hand to a colleague, and put in a portfolio. That is the difference between a student exercise and science, and your solver is about to cross it.

In this chapter, you will learn to:

• Assemble the whole solver — the six canonical modules, the parallel back-ends, the I/O — behind the single step interface that has been frozen since Chapter 6, and see why a stable interface is what let forty chapters compose.
• State the problem and method the way a paper does: the heat equation, the five-point stencil, the explicit FTCS scheme, and the stability limit, written for a reader who was not in the room.
• Verify the solver against an exact analytical solution, run a convergence study, and confirm it recovers the second-order accuracy the Chapter 22 theory predicts — the single most important thing you can do to earn trust.
• Analyze performance without fooling yourself: a serial baseline, Amdahl's ceiling, a scaling table, and the roofline reasoning that explains why this kernel is limited by memory bandwidth, not arithmetic.
• Present it — the figures, the paper structure, and what a reviewer is actually looking for.

Learning Paths

How to read this chapter by track. - 🔬 Scientist — this is the chapter you have been working toward. Read §38.2–§38.3 with a pencil; the verification is the part of your future job that most separates a trusted result from a retracted one. - 📖 Standard — read straight through. It is a synthesis: almost every earlier chapter reappears here, doing the job it was built for. Watching the pieces click together is the reward for the whole book. - 🔧 Legacy — §38.6 (presentation and reproducibility) is the discipline that keeps a validated old code trusted for forty years; §38.3's regression-against-analytical is exactly what protects the science in a code you inherit when you modernize it. - ⚡ HPC — §38.4 is yours: Amdahl, strong vs weak scaling, and the roofline. Read it beside Chapter 31; the honest framing of scaling numbers is a professional skill, not a formality.

38.1 Assembling the Pieces

For most of this book the solver has grown one file at a time, and you may not have stepped back to see the whole. Here it is. Five roles, six modules, one driver — the layered package you built in Chapter 36, now carrying real physics, tests, and parallel back-ends:

heat-solver/
├── fpm.toml                     the build manifest (records the build configuration)
├── README.md                    what it solves; how to build, run, and reproduce
├── src/                         THE LIBRARY
│   ├── kinds.f90                    UTILITY: the dp precision kind            (Ch. 3)
│   ├── timers.f90                   UTILITY: tic/toc around system_clock      (Ch. 28)
│   ├── heat_types.f90               CORE:    the field_t derived type         (Ch. 9)
│   ├── heat_solver.f90              SOLVER:  laplacian + step + stable_dt      (Ch. 24)
│   └── heat_io.f90                  I/O:     read_config, write_field, write_vtk (Ch. 7, 26)
├── app/
│   └── main.f90                     DRIVER:  program heat — setup, time loop, output
└── test/
    └── test_solver.f90              regression vs the analytical solution     (Ch. 37)

Nothing in that tree is new to you. What is new is the view: a program that a maintainer, a reviewer, or the you-of-next-year could open cold, read the README, glance at the module map, and understand in minutes. That readability is not decoration — it is what makes the code reviewable, and a result you cannot review is not yet science.

The one idea that made this possible. Look back at how much changed under the hood between Chapters 24 and 34. The body of step went from a serial array-section update, to a cache-tuned loop, to an OpenMP parallel do, to an MPI halo exchange. And through every one of those rewrites, its signature never moved:

subroutine step(field, alpha, dt)     ! frozen since Chapter 6; body rewritten five times
  type(field_t), intent(inout) :: field
  real(dp),      intent(in)    :: alpha, dt

Because the driver only ever knew step by that signature, we could swap a serial body for a parallel one without touching a single line of the driver, the I/O, or the tests. The interface was the contract; the implementation was free to change. This is the payoff of a discipline the book has pressed since Chapter 6's intent and Chapter 8's module boundaries, and it is worth naming as the idea it is.

🚪 Threshold Concept — a stable interface is what lets a program grow. The reason forty chapters could keep bolting capability onto one solver is that each piece talked to the others only through frozen signatures — step(field, alpha, dt), laplacian(u, dx, dy), write_vtk(field, filename, step). Change a body, and everything that calls it is undisturbed; change a signature, and the whole edifice shifts. Real scientific codes live for decades precisely because their internal interfaces are treated as promises. Once you see the interface — not the code behind it — as the unit of design, you can build a system too large to hold in your head, because you only ever have to hold one interface at a time.

So "assembling the pieces" is, gratifyingly, almost no work. The modules already compose; fpm already derives the compile order from the use statements; the parallel step already produces the identical answer the serial one did (you proved that in Chapter 33, digit for digit). What remains is not more code. It is the three things that turn a working program into a result: verify it, measure it, and present it. The rest of the chapter is those three things.

🔗 Connection: the "same interface, many bodies" pattern is exactly how production codes ship a serial path for a laptop and a GPU path for a supercomputer from one source tree. WRF, for instance, selects among serial, OpenMP, and MPI builds of the same dynamics routines at configure time. Your step with a serial, an OpenMP, and an MPI body (Chapters 33–34) is a small, honest rehearsal of that architecture — and a demonstration that modern Fortran is a modern language, one whose module and interface facilities scale from a class project to a climate model.

38.2 Problem Statement and Numerical Method

Here begins the write-up. From this section on, read the prose as if it were the methods of a short paper — because that is exactly what we are drafting. A methods section has one job: to let a competent stranger reproduce your computation. That means stating the equation, the discretization, the scheme, and the stability constraint precisely enough that no guessing is required.

The governing equation. We solve the two-dimensional heat (diffusion) equation for a temperature field $u(x, y, t)$ on a square domain $\Omega = [0, L]\times[0, L]$:

$$ \frac{\partial u}{\partial t} = \alpha \nabla^2 u = \alpha\left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right), $$

where $\alpha > 0$ is the thermal diffusivity. The edges are held at fixed temperatures (Dirichlet boundary conditions), and an initial field $u(x, y, 0) = u_0(x, y)$ is prescribed. This is the problem you chose in Chapter 1 and made real in Chapter 24; we restate it here because a paper never assumes its reader has read the previous thirty-seven chapters.

Spatial discretization. We lay a uniform structured grid of spacing $h = \Delta x = \Delta y$ over $\Omega$ and store the temperature at the grid points as field%u(i,j). The continuous Laplacian is replaced by the five-point stencil, the sum of two central second differences:

$$ \nabla^2 u\big|_{i,j} \approx \frac{u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} - 4\,u_{i,j}}{h^2}. $$

This approximation is second-order accurate in space: its truncation error is $O(h^2)$, a fact Chapter 22 derived from a Taylor expansion and one we will measure in §38.3. The $1/h^2$ scaling is what turns a comparison of neighbours into an actual second derivative — drop it and the result changes when you refine the mesh, the signature of a broken discretization.

Time integration. We march in time with the explicit forward-Euler scheme (Forward-Time, Centred-Space, or FTCS), evaluating the stencil on the current field and stepping directly to the next:

$$ u^{n+1}_{i,j} = u^n_{i,j} + r\left(u^n_{i+1,j} + u^n_{i-1,j} + u^n_{i,j+1} + u^n_{i,j-1} - 4\,u^n_{i,j}\right), \qquad r \equiv \frac{\alpha\,\Delta t}{h^2}. $$

The dimensionless group $r$ (the diffusion number) collects the physics, the timestep, and the grid into one number, and it governs everything that follows.

Stability. The explicit scheme is only conditionally stable. A von Neumann analysis (Chapter 24, §24.4) shows the worst-case (checkerboard) mode is amplified each step by $G = 1 - 8r$, so stability — $|G| \le 1$ — requires

$$ r = \frac{\alpha\,\Delta t}{h^2} \le \frac{1}{4} \qquad\Longleftrightarrow\qquad \Delta t \le \frac{h^2}{4\alpha}\quad(\text{2D}). $$

Our code never guesses the timestep; it computes one safely under this limit with stable_dt, so refining the grid automatically shrinks $\Delta t$ (as $\Delta t \sim h^2$) and the run stays stable by construction. That single design choice — derive the timestep from the grid — is the difference between a solver a reviewer trusts and one that blows up the first time someone changes the resolution.

💡 Intuition: a good methods section reads like a recipe a stranger could follow to get your dish. Equation, grid, stencil, scheme, timestep rule — five sentences, and someone in another lab could rebuild your solver from scratch and expect the same numbers. If a reader would have to guess any of those five, the section is not yet finished. Vagueness in a methods section is not modesty; it is the enemy of reproducibility.

🔄 Check Your Understanding. 1. Why does a paper restate the governing equation and boundary conditions even though "everyone knows the heat equation"? 2. In one sentence, what does the diffusion number $r$ combine, and what is its stability limit in 2D? 3. Why does deriving $\Delta t$ from $h$ (rather than hard-coding it) make the solver safe under grid refinement?

Answers

1. Because reproducibility demands it — the specific $\alpha$, domain, boundary values, and initial condition define your problem, and "the heat equation" alone is not enough to reproduce a number. 2. $r = \alpha\Delta t/h^2$ bundles diffusivity, timestep, and spacing; the explicit 2D scheme is stable only for $r \le 1/4$. 3. Refining the grid halves $h$, which quarters the stability-limited $\Delta t$; computing $\Delta t$ from $h$ enforces $r \le 1/4$ automatically, so a resolution change can never silently push the run across the CFL cliff.

38.3 Verification and Validation

Now the heart of the matter, and the section a good reviewer reads first. You have a program that produces plausible-looking heat maps. So what? Plausible is not correct. A subtly wrong sign, a factor of two in the scaling, an off-by-one at the boundary — any of these produces output that looks like diffusion and is wrong. The question that separates science from wishful thinking is: how do you know the numbers are right? The professional answer has a name and a structure.

Definition (verification and validation, V&V). Verification asks "are we solving the equations right?" — does the code correctly and accurately solve the mathematical model we intended, free of bugs, converging to the true solution of the equations at the expected rate? Validation asks "are we solving the right equations?" — does the mathematical model itself match physical reality, as measured by experiment or observation? The two are different and both are necessary: a perfectly verified code can faithfully solve the wrong physics, and a physically apt model can be wrecked by a coding bug. Together they are abbreviated V&V, and in computational science a result without them is not yet a result.

The distinction matters for what we are about to do, so be precise about it. Comparing our solver to an exact solution of the same heat equation is verification: it confirms the code solves the equation correctly. It is not validation — validation would compare the simulation to a real heated plate in a laboratory, testing whether the heat equation with our chosen $\alpha$ actually describes that plate. This chapter does verification, thoroughly, and is honest that full validation needs experimental data we do not have. Saying so is not a weakness of the write-up; it is the kind of precision that makes the rest of it credible.

📜 From History: the "solving the equations right / solving the right equations" phrasing was popularized in the computational-fluid-dynamics community (notably by Patrick Roache) and is now enshrined in the verification-and-validation standards of professional engineering societies. It exists because the field learned, sometimes expensively, that a simulation nobody verified is a very confident way to be wrong — and that the two failure modes (a bug versus a bad model) need different cures.

An analytical solution to check against

Verification needs a case whose exact answer you know in closed form. The heat equation obliges. On the unit square $[0,1]\times[0,1]$ with all four edges held at zero, the function

$$ u(x, y, t) = \sin(\pi x)\,\sin(\pi y)\; e^{-2\alpha\pi^2 t} $$

is an exact solution. You can verify it by substitution in three lines: each spatial second derivative brings down a factor $-\pi^2$, so $\nabla^2 u = -2\pi^2 u$ and $\alpha\nabla^2 u = -2\alpha\pi^2 u$; the time derivative brings down exactly $-2\alpha\pi^2$; the two sides match. It satisfies the zero-Dirichlet edges (each sine vanishes at $0$ and $1$) and starts from $u_0 = \sin(\pi x)\sin(\pi y)$. Physically, it is the slowest-decaying mode of the plate: a single smooth hump that keeps its shape and simply fades, its amplitude decaying by the fixed rate $2\alpha\pi^2$. That is the reference truth we will hold the code to.

Note that this is a different problem from the Chapter 1 plate (one hot edge, three cold), which has no simple closed form. This is standard practice: you verify on a problem chosen for its known answer, then present results on the problem you actually care about. Reaching for a solvable case to test the machinery is not cheating — it is the entire method.

💡 Intuition: the sine mode is the perfect test case because it is an eigenfunction of the Laplacian — the operator returns the same shape, scaled. So the exact solution never changes shape; only its height shrinks, at a rate you can write down. Any discrepancy your code shows against it is therefore pure numerical error, cleanly separated from the physics.

The convergence study

Here is the beautiful part, and the reason this particular test is worth its weight. Sample the exact initial mode onto the grid, $u^0_{i,j} = \sin(\pi x_i)\sin(\pi y_j)$. That discrete field is an exact eigenvector of the five-point stencil: applying the discrete Laplacian returns the same grid pattern, scaled by the discrete eigenvalue

$$ \lambda_h = -\frac{8}{h^2}\sin^2\!\left(\frac{\pi h}{2}\right). $$

Because the mode is an eigenvector, the FTCS update multiplies the entire field by a single number each step — the amplification factor $G = 1 + \alpha\Delta t\,\lambda_h = 1 - 8r\sin^2(\pi h/2)$ — so after $K$ steps the numerical solution is exactly $u^K_{i,j} = G^K\,u^0_{i,j}$. No other modes are excited, no round-off analysis is needed: the error at time $T = K\Delta t$ is the clean difference between the discrete and exact decay of one number,

$$ \text{error}(h) = \big|\,G^K - e^{-2\alpha\pi^2 T}\,\big|, $$

which we can compute by hand and which the code must reproduce. Fix $\alpha = 1$, hold $r = 0.2$ (safely under $1/4$), march to $T = 0.05$, and refine the grid. Each halving of $h$ quarters $\Delta t$ (to keep $r$ fixed) and so quadruples the step count — the $\Delta t \sim h^2$ tax of explicit diffusion — and the maximum error over the plate falls like $h^2$:

Read the last two columns. As the grid refines, the error drops by a factor approaching four each time $h$ is halved, and the observed order — $\log_2(\text{ratio})$ — settles on 2. That is the whole game: the code recovers the second-order spatial accuracy the Chapter 22 theory predicts for the five-point stencil. (The coarsest grid is outside the asymptotic regime — one interior point is barely a simulation — which is why its ratio overshoots; the trend is what matters, and it converges cleanly to 2.) The numbers above are exact, computed from the closed-form $G^K$; your program will reproduce them, and that agreement is the verification.

🚪 Threshold Concept — a convergence study, not a single run, is what proves correctness. A solver that gives a "reasonable" answer on one grid has proven nothing: a bug and a correct code can both look reasonable. What a bug almost never survives is a convergence study. If refining the mesh drives the error down at exactly the theoretical rate — order 2 here — then the code is discretizing the right operator correctly, because a wrong stencil converges at the wrong order (or not at all). The observed order of accuracy is the single most diagnostic number in computational science. When a reviewer asks "how do you know it's right?", the convergence table is the answer, and almost nothing else is.

⚠️ Common Pitfall — the wrong order is a specific bug signature. A convergence study does more than pass or fail; it diagnoses. Forget the $1/h^2$ scaling and the error will not converge at all. Drop the stencil to a lopsided three-point form and you will measure order 1, not 2. Overwrite the field in place (mixing time levels, the §24.3 trap) and the order degrades. Measuring order 2 when you expected order 2 is strong evidence the discretization is correct; measuring order 1 is not "close enough" — it is a bug announcing its own address.

This is also exactly what the Chapter 37 regression test guards: it pins one point of this convergence curve (the code's error on a fixed grid) so that the day a careless refactor breaks the physics, the test fails loudly instead of the science drifting silently. The analytical solution is both the verification and the regression oracle — one idea doing two jobs.

🔄 Check Your Understanding. 1. What is the difference between verification and validation, and which one does this section perform? 2. Why is the sine mode $\sin(\pi x)\sin(\pi y)$ an especially convenient function to verify against? 3. Your convergence study measures order 1 where you expected order 2. Is that "close enough"? What does it most likely mean?

Answers

1. Verification asks whether the code solves the chosen equations correctly (compare to an exact solution of the same PDE); validation asks whether those equations describe reality (compare to experiment). This section does verification. 2. It is an eigenfunction of the Laplacian, and its grid samples are an exact eigenvector of the five-point stencil — so the numerical solution is exactly $G^K$ times the initial field, isolating pure numerical error with no other modes to muddy it. 3. No — it is a bug. A first-order result where the method is second-order almost always means a broken stencil (a dropped scaling, a lopsided difference) or mixed time levels. The wrong order is a diagnostic, not a rounding issue.

38.4 Performance Analysis

A verified solver is correct. A useful one is also fast enough to run the problem you care about before the deadline — and reporting its performance honestly is its own discipline, with its own ways to fool yourself. This section establishes a baseline, applies the theory, and — crucially — states plainly which numbers are measurements and which are expectations. The frame is the theme all of Part VII pressed: performance is not accidental — it is Amdahl's ceiling, the roofline's walls, and the memory layout of Chapter 5, reasoned about deliberately rather than hoped for.

⚠️ A word on the numbers in this section. This book never runs code, so every timing and speedup below is illustrative — an order-of-magnitude expectation, not a measurement. They are shaped to be realistic and to teach the right reasoning, but the only honest performance number is one you measure on your machine with the profiler of Chapter 28. Run it yourself; treat the tables as a map, not the territory. What is exact here is the arithmetic — the Amdahl formula, the flop counts — and those you can trust.

The serial baseline. Everything starts from one honest measurement: how long the tuned single-core solver takes. On a moderate grid — say $1000\times1000$, marched a few thousand steps — a representative single-core run might spend, per step, essentially all of its time in the interior stencil sweep, with setup, boundary handling, and (infrequent) VTK output a rounding error by comparison. That profile is the whole reason the solver parallelizes well, and it is the number every speedup is measured against. A speedup with no stated baseline is meaningless; always report what you compared to.

Amdahl's ceiling. Chapter 31 gave us the law that bounds any parallel effort: if a fraction $f$ of the work is parallelizable and $1 - f$ is irreducibly serial, then $p$ processors can speed the whole up by at most

$$ S(p) = \frac{1}{(1 - f) + f/p}. $$

For our solver the parallel fraction is the stencil sweep and the serial remainder is the time loop's bookkeeping (the buffer swap, the step counter — step $n+1$ genuinely needs step $n$, so the loop itself cannot be parallelized). Take $f = 0.98$ as a representative estimate for a real run. The exact ceilings follow from the formula alone:

Two things to read off it. First, the ceiling: with $f = 0.98$, no number of cores beats $1/(1-f) = 50\times$ — the serial 2% is the wall. Second, the erosion: efficiency falls as you add cores, because the fixed serial cost becomes a larger share of a shrinking total. These are exact consequences of the formula, and they are the honest expectation your measured strong-scaling curve should be compared against.

Definition (roofline model; arithmetic intensity). The arithmetic intensity of a kernel is the number of floating-point operations it performs per byte of data it moves to and from memory, in flops/byte. The roofline model plots achievable performance against arithmetic intensity: a kernel with low intensity is capped by a sloped memory-bandwidth ceiling (it starves waiting for data), while a kernel with high intensity is capped by the flat peak-compute ceiling (it is limited by the arithmetic units). The name comes from the plot's shape — a slanted line meeting a horizontal one, like a roofline. Which ceiling you are under tells you what to optimize: bandwidth-bound kernels want better memory access, compute-bound kernels want better vectorization.

Where does our stencil sit? Count the arithmetic: each interior cell costs about ten flops (a handful of adds and multiplies for the two scaled second differences and the update). Count the traffic: to update one cell the machine must read its neighbours and write the result — and even with perfect cache reuse, the kernel moves on the order of tens of bytes per cell. That puts its arithmetic intensity well below one flop/byte — squarely in memory-bound territory, far under the compute ceiling. The consequence is the honest correction to the Amdahl table above:

⚡ Performance Note — Amdahl is the optimist; memory bandwidth is the realist. Amdahl's $50\times$ assumes the only limit is the serial fraction. But a memory-bound kernel hits a second wall: once enough cores are sweeping the grid to saturate the shared memory bandwidth, adding more cores buys almost nothing, because they all queue for the same starved memory bus. So a real strong-scaling curve for this stencil typically tracks Amdahl for the first few cores and then plateaus below it — not because the parallelism is wrong, but because the hardware runs out of bandwidth before it runs out of cores. This is precisely why Chapter 29 tuned the loop order and cache blocking: for a bandwidth-bound kernel, using the memory you fetch well matters more than adding arithmetic units. The roofline tells you which wall you are about to hit; the profiler tells you when.

Strong versus weak scaling. Chapter 31 drew the distinction and it is the right lens for the two parallel back-ends. Strong scaling — fixed total problem, more processors — is what OpenMP does on one node, and it is bounded by Amdahl and then by bandwidth, as above. Weak scaling — the problem grows with the processor count, so each keeps the same local work — is the natural regime for the MPI solver: give every rank its own $500\times500$ strip, and a run on 64 ranks solves a $4000\times4000$ plate in roughly the time one rank solves $500\times500$, because each does constant work and the only overhead is the thin halo exchange at the strip boundaries. Distributed solvers are usually sold on weak scaling for exactly this reason: you do not use a thousand nodes to solve today's problem faster; you use them to solve a problem a thousand times larger. Which curve you report depends on which question you are answering — and stating which is part of reporting honestly.

🐍 Python Comparison: you could prototype and plot all of this scaling analysis in Python — read the timings your Fortran run writes, and let matplotlib draw the speedup and roofline curves. That is the sixth theme of the book, Fortran and Python are better together, in its natural habitat: Fortran runs the ten-thousand-step kernel at full speed, and Python does the analysis and the figures. The one thing you must not do is time the kernel from Python — the interpreter's overhead would swamp exactly the microseconds you are trying to measure. Time in Fortran with system_clock (Chapter 28); analyze in Python.

38.5 Results and Visualization

With the solver verified and its performance characterized, we can finally show what it does — and this is where the Chapter 26 VTK output pays off. We present results on the problem we actually care about: the Chapter 1 plate, hot on its top edge at $100°$, cold ($0°$) on the other three, released from a cold interior and left to reach steady state.

The transient. Early in the run, heat has only just begun to leave the hot edge, and the interior is a sharp front creeping downward and sideways. You have already computed the first two steps of the small-grid version by hand, and it is worth seeing again as the concrete anchor beneath the pictures — the $5\times5$ plate after two steps, the warmth spreading symmetrically inward from the top row:

after step 2 (5x5 plate, hot top edge, dt = 0.2):
  100.00  100.00  100.00  100.00  100.00
    0.00   28.00   32.00   28.00    0.00
    0.00    4.00    4.00    4.00    0.00
    0.00    0.00    0.00    0.00    0.00
    0.00    0.00    0.00    0.00    0.00

On a real $500\times500$ grid this same physics produces a smooth thermal front. Opening the VTK time series in ParaView and pressing play, you watch the hot edge bleed into the cold interior: a bright band along the top that softens and reaches downward, frame by frame, until the whole plate settles.

Figure 38.1 (described). The transient, as ParaView renders it. Four frames of the $500\times500$ plate, colored by temperature with a perceptually uniform colormap (viridis — never jet, per §26.5). Frame 1 (early): a thin brilliant strip along the top edge, the rest deep blue-black — heat has barely entered. Frame 2: the strip has thickened into a downward-bulging tongue, warm in the upper third. Frame 3: the warm region reaches past the middle, its contours smooth arcs bowing away from the hot edge. Frame 4 (near steady state): a stable, unchanging gradient from the hot top to the three cold edges, brightest at the top and darkening smoothly downward — the plate has stopped evolving. A play button at the top of the ParaView window scrubs through all of them; a .pvd collection labels each frame with its physical time in seconds.

The steady state. Left to run, the plate reaches the state where nothing changes anywhere — $\partial u/\partial t = 0$, so $\nabla^2 u = 0$, Laplace's equation. The result is a smooth, unchanging temperature field: hottest against the top edge, falling away toward the three cold edges, with gently curved isotherms. It is the picture the whole project was aiming at from Chapter 1 — and the moment the number in your array and the picture in your head finally become the same thing.

The verification figure. A paper's results section shows not only the pretty picture but the evidence. The convergence table of §38.3 is that evidence, and it has a canonical visual form: a log–log plot of error versus grid spacing $h$, on which second-order accuracy appears as a straight line of slope 2. Drawing that line — your four error points falling on it — is the single most persuasive figure in the whole write-up, because it shows, at a glance, that the code does what the mathematics promised.

Figure 38.2 (described). The convergence plot. Error (max norm, log scale) on the vertical axis against grid spacing $h$ (log scale) on the horizontal. The four data points from the §38.3 table fall very nearly on a straight line; a reference triangle of slope 2 is drawn beside them for comparison, and the points parallel it. Caption: "Maximum error against the analytical solution versus grid spacing at $T = 0.05$; the slope-2 trend confirms the expected second-order spatial accuracy." This one figure is worth more to a skeptical reviewer than every heat map combined.

🔗 Connection: the division of labor here is the one real codes use. Fortran computes the field and writes VTK; ParaView renders the heat maps and the animation; matplotlib (driven by a few lines of Python reading the error table) draws the convergence plot for the paper. Each tool does the half it is best at — arrays are Fortran's superpower for the computation, and the visualization ecosystem is Python's. You do not choose between them; you use both.

38.6 Presenting It: The Structure of a Computational-Science Paper

You have the science. The last skill — the one that turns a result into something the field can use — is presenting it. A computational-science paper has a conventional structure, and it is conventional for a good reason: it is the order in which a reader needs the information to decide whether to believe you. Learn the skeleton once and you can read any paper in the field faster and write your own without staring at a blank page.

What reviewers actually look for. Having refereed and been refereed, a computational scientist reads a paper with a short mental checklist, and it is worth knowing what it is, because writing to it is how you get accepted:

• Is it correct? This is first and it is non-negotiable. A reviewer looks straight for the V&V: is there a convergence study, does it hit the expected order, is the comparison to a genuine reference? A paper with a beautiful method and no verification is a paper a reviewer cannot trust, and will not.
• Is it reproducible? Are the equations, discretization, parameters, compiler, flags, and inputs stated precisely enough to rebuild the result? This is the Chapter 37 discipline — record the build configuration, pin the inputs, publish the code — read from the other side of the table. "Works on my machine" is not a result; a result travels.
• Are the claims supported by the evidence? If you claim second-order accuracy, is there a slope-2 plot? If you claim it scales, is there a scaling curve with a stated baseline? Reviewers are allergic to claims that outrun their evidence — "blazingly fast" with no number, "scales well" with no curve.
• Is it honest about limitations? Every method has a regime where it fails. Ours is explicit, so it pays the $\Delta t \sim h^2$ tax on fine grids and would want an implicit scheme (the Chapter 21 dgesv route) for stiff, high-resolution runs. Naming that strengthens the paper: it shows you understand your own tool. The honesty about limits is not a confession; it is a credential.
• Is it significant? Does it teach the reader something — a method, a result, a reusable code? Significance is the hardest bar and the one furthest outside your control, but it is much easier to clear on a foundation of correctness, reproducibility, and honest scope.

💡 Intuition: notice that four of the five reviewer questions are about trust, not cleverness. Correct, reproducible, supported, honest — only the fifth, significance, is about the idea. That ratio is the real lesson of this chapter: in computational science, the credibility of a result is built, deliberately and visibly, and the building is the work. The cleverest method nobody can trust is worth nothing; a modest one, verified and reproducible, is worth citing.

🔄 Check Your Understanding. 1. In what order does a computational-science paper present its material, and why is verification placed before results? 2. Name three of the five things a reviewer checks, and say which single one is non-negotiable. 3. Why does openly stating a method's limitations tend to strengthen a paper rather than weaken it?

Answers

1. Problem → method → V&V → implementation/performance → results → conclusion → reproducibility. Verification comes before results because a reader must be convinced the code is correct before the results mean anything. 2. Any three of: correctness (V&V), reproducibility, claims-supported-by-evidence, honesty about limitations, significance. Correctness is non-negotiable — without it nothing else matters. 3. It shows the author understands the tool's regime of validity, which builds trust; a paper that pretends its method has no limits invites the reviewer to find the one it ignored.

Project Checkpoint

This is the one every checkpoint has pointed to. Your solver is finished — and "finished," for a scientific code, means presented. So the capstone deliverable is not more code; it is the whole thing written up as a short computational-science paper, with the code as its reproducible core. Here is the shape of that write-up, assembled from the six sections above, in miniature.

A Verified, Parallel Finite-Difference Solver for the 2D Heat Equation in Modern Fortran

Abstract. We present a modular Fortran 2018 solver for the two-dimensional heat equation on a structured grid, using a five-point spatial stencil and an explicit FTCS time scheme with a CFL-safe timestep. The code is verified against an analytical solution; a convergence study confirms second-order spatial accuracy (observed order → 2). The solver parallelizes behind a single step interface via OpenMP (shared memory) and MPI (domain decomposition), producing bit-identical results to the serial path. We report Amdahl-bounded strong scaling and note the kernel is memory-bandwidth-bound. Output is written as VTK for ParaView.

Method. The equation $\partial u/\partial t = \alpha\nabla^2 u$ is discretized with the five-point Laplacian ($O(h^2)$) and marched explicitly, $u^{n+1} = u^n + r(\text{stencil})$, with $r = \alpha\Delta t/h^2 \le 1/4$ for stability; the timestep is derived from the grid by stable_dt. Dirichlet boundaries are held fixed; the interior is swept each step.

Verification. Against $u = \sin(\pi x)\sin(\pi y)e^{-2\alpha\pi^2 t}$ on the unit square with zero edges, the max error at $T = 0.05$ ($\alpha=1$, $r=0.2$) falls by ~$4\times$ per grid halving — observed order 2, as the five-point stencil's truncation error requires.

Results & performance. The Chapter 1 plate reaches a smooth steady state (Figure 38.1); the convergence plot is a slope-2 line (Figure 38.2). With a parallel fraction $f \approx 0.98$, Amdahl bounds the shared-memory speedup near $50\times$; the memory-bound stencil plateaus below that as bandwidth saturates. The MPI solver weak-scales — constant work per rank as the plate grows.

Reproducibility. Build with fpm build --profile release (gfortran -std=f2018 -O3 -march=native, or add -fopenmp; MPI via mpif90). The test/ suite regresses against the analytical solution. Complete source is in Appendix I.

Your deliverable. Produce this document for your solver: the six sections, the two figures (steady state and convergence plot), the results table, and a README that lets a stranger reproduce it. The assembled program — the six canonical modules behind the frozen step interface, with the analytical-solution verification wired into test/ — is code/project-checkpoint.f90, and the full annotated source is reproduced in Appendix I. The self-contained checkpoint program runs the canonical $5\times5$ demonstration (reproducing the exact 28/32/28, 4/4/4 field you have hand-computed since Chapter 24) and the one-step analytical check on the $5\times5$ sine mode, printing the numerical value ($0.7657$) beside the exact ($0.7813$) so the verification is visible in the output itself. Physics that is correct, code that is organized, results that are presented — that is the whole of it, and it is now yours.

This is what you built: not a toy, but a small, complete, honest piece of computational science. The next two chapters send it into the world — Chapter 39 refreshes one of its snippets with the newest standard, and Chapter 40 turns it into a portfolio piece you can show an employer.

Summary

This chapter assembled the running solver into a finished scientific result and taught the discipline of presenting it — the half of computational science that is not code.

The two things to remember. First, a convergence study against a known solution is how you prove a solver is correct — the observed order of accuracy is the most diagnostic number you can report, and no heat map, however pretty, substitutes for it. Second, in computational science the result is not the code; it is the code plus the argument that it is right — verified, reproducible, honestly scoped — and building that argument, visibly, is the work.

Spaced Review

A broad synthesis, fittingly for the capstone: this pulls on the three chapters it most directly assembles — the PDE core (Chapter 24), OpenMP (Chapter 33), and testing and reproducibility (Chapter 37). Answer before peeking.

1. (Ch. 24) Our convergence study holds $r = \alpha\Delta t/h^2 = 0.2$ fixed and halves $h$. Why does the number of timesteps to reach a fixed final time $T$ quadruple each time, and what property of the explicit scheme forces that?
   

   Answer

   Holding $r$ fixed means $\Delta t \propto h^2$, so halving $h$ quarters $\Delta t$ and quadruples the step count $K = T/\Delta t$. It is forced by the CFL stability limit $r \le 1/4$: the explicit scheme's timestep must shrink as $h^2$, the notorious $\Delta t \sim h^2$ tax that makes explicit diffusion expensive on fine grids and motivates implicit methods.

2. (Ch. 24) The verification relies on the discrete sine mode being an exact eigenvector of the five-point stencil, so the numerical field is exactly $G^K u^0$. What is $G$ in terms of $r$ and $h$, and why does that make the error trivial to compute?
   

   Answer

   $G = 1 - 8r\sin^2(\pi h/2)$ (from the discrete eigenvalue $\lambda_h = -\frac{8}{h^2}\sin^2(\pi h/2)$). Because the mode never changes shape — only its amplitude scales by $G$ each step — the whole-field error at time $T = K\Delta t$ collapses to the one-number difference $|G^K - e^{-2\alpha\pi^2 T}|$, no round-off or multi-mode analysis needed.

3. (Ch. 33) The OpenMP step produces bit-identical results to the serial step. Why is that determinism required, not merely nice — and what one scoping mistake in the parallel loop would break it, silently?
   

   Answer

   The physics is deterministic, so a correct parallel solver *must* match the serial one exactly — if it does not, you have a bug, not a faster answer. The classic silent breaker is failing to make the inner loop index `private` (or writing back into the field in place instead of a separate buffer): both introduce a data race, giving different, run-dependent numbers. `default(none)` forces you to scope every variable and catches it at compile time.

4. (Ch. 37) The analytical solution serves as both the verification reference and the regression-test oracle. Explain the difference in how each uses it, and why pinning one point of the convergence curve as a test protects the science.
   

   Answer

   Verification uses the analytical solution *once*, across many grids, to establish the order of accuracy (a study). The regression test uses it *repeatedly*, on one fixed grid, to assert the code still reproduces a known error to tolerance on every commit. Pinning that one point means the day a refactor breaks the physics, the test fails loudly — so correctness, once earned, cannot drift away unnoticed.

5. (Ch. 37) A colleague reports "the solver is 8× faster on 8 cores" with no other detail. Name two pieces of information required before that number means anything, and connect your answer to reproducibility.
   

   Answer

   At minimum: (1) the **baseline** it is measured against (8× versus *what* — an untuned serial run, or the best serial code?), and (2) the **build configuration and problem size** (compiler, flags, grid, step count, machine). Without them the speedup is unreproducible and possibly meaningless — a fast-but-unoptimized baseline inflates any speedup. Reproducibility demands recording the baseline, the build, and the inputs, exactly the Chapter 37 discipline.

What's Next

Your solver is complete: correct, organized, parallel, verified, and written up. There are two things left to do with it, and they are the subject of Part X. First, a glance forward: Chapter 39 takes one snippet of your finished solver and refreshes it with Fortran 2023 — conditional expressions, cleaner enumerations — and shows you the living, still-evolving language your code is written in, so you can see that Fortran is not dead in the most direct way possible: by writing tomorrow's Fortran today. Then the closing chapter, Chapter 40, turns this project into a portfolio piece — the README, the license, the figure, the one-paragraph abstract you just drafted — and shows you where the people who write Fortran work, and why they cannot find enough of you. You started, forty chapters ago, with a decision about a square plate. You are ending with a piece of computational science you built, understand completely, and can defend. Let's take it into the world.