Chapter 37: Testing, Documentation, and Software Engineering for Scientific Fortran

"Program testing can be used to show the presence of bugs, but never to show their absence." — Edsger W. Dijkstra

Overview

In Chapter 36 you learned to navigate a large scientific code and reshaped your own solver into a real package — a src/ + app/ + test/ tree with a build system. You can now find your way around a hundred thousand lines, and your solver looks like software an engineer would respect. But looking like software is not the same as being trustworthy software, and a navigable code is not yet a correct one. Two hard questions remain, and this chapter answers both.

The first is how do you know it still works? You refactored the solver in Chapter 36; you will optimize it, and parallelize it, and someone after you will change it again. Every one of those edits is a chance to silently break the physics — to compute a plausible-looking wrong answer that no compiler error will ever flag. In ordinary software a bug crashes or returns nonsense you can see. In numerical software the most dangerous bugs return smooth, believable numbers that are simply incorrect, and they can sit undetected in a published result for years. The defense is automated testing: a suite of checks that runs on every change and screams the moment the answer moves.

The second question is how does anyone else — including you, in three years — reproduce what you did? A scientific result is not just a number; it is a number plus everything needed to regenerate it: the exact code, the compiler and its flags, the inputs, the library versions, the random seeds. Leave any of those unrecorded and you get the most quietly corrosive failure in all of computational science — the result that "worked on my machine" and works nowhere else, not even on your own machine after an upgrade. This chapter is about closing that gap: testing so the code is right, documenting and version-controlling so it is understandable, continuous integration so it stays right automatically, and reproducibility so a result is a result and not an anecdote.

In this chapter, you will learn to:

• Write unit tests for numerical code — and know what is worth testing when the "right answer" is itself an approximation: exact solutions, invariants, symmetry, and edge cases.
• Build regression tests that pin a solver's output against a known-good result, and decide between comparing bit-for-bit and comparing within a tolerance.
• Verify the solver against an analytical solution, turning the convergence experiment of Chapter 22 into a test.
• Set up continuous integration with GitHub Actions so every push builds and tests the code across several compilers, and catches a break before it reaches anyone else.
• Document the code with FORD doc comments, version-control it with git the way scientific software demands, and record enough provenance that a run is reproducible years later.

Learning Paths

How to read this chapter by track. - 🔬 Scientist ("I need to trust my own results") — this whole chapter is for you, but §37.1 (what to test in numerical code), §37.2 (regression against the analytical solution), and §37.5 (reproducibility) are the ones that will change how you work. The Project Checkpoint is the habit to build. - 📖 Standard — read straight through. Testing, docs, and CI are engineering topics the language reference cannot teach you, and they are what turns Chapter 36's package into the capstone's research artifact. - 🔧 Legacy ("I inherited old code I'm afraid to touch") — §37.2 is your safety net: a regression test that pins the old code's output is exactly what lets you modernize without fear, the discipline Chapter 18 leaned on. Read it first. - ⚡ HPC ("I run big parallel jobs") — §37.2's bit-for-bit-vs-tolerance discussion is your problem, because a parallel reduction reorders floating-point sums; §37.5 (record the build and the run) is how you make a thousand-core result reproducible. Read the ⚠️ and ⚡ notes closely.

37.1 Unit Testing: What to Test in Numerical Code

Start with the smallest unit of trust. Your solver is built from procedures — laplacian, step, stable_dt — and if each one is individually, provably correct, then a bug in the whole is a bug in how they are wired together, a much smaller place to look. Testing each procedure in isolation is unit testing.

Definition (unit test). A unit test is a small, automated check that exercises one unit of code — a single procedure, or a small cluster of them — in isolation, with known inputs, and verifies that it produces the expected output. "Automated" is the load-bearing word: a unit test is a program that returns pass or fail with no human reading the numbers, so it can be run a thousand times, by anyone, on every change, for free. A good unit test is fast (milliseconds), independent of the others, and specific enough that when it fails you know roughly where the bug is, not merely that there is one.

The mechanism underneath a unit test you already have: the assertion from Chapter 13 — a check that a condition you believe is true actually is, failing loudly if not. A unit test is an assertion with a name, collected into a suite. You can build a serviceable test harness in a couple of dozen lines of standard Fortran, which is exactly the right thing to do when you want to understand what a test framework is before adopting one, or when you are on a machine where installing one is more trouble than it is worth:

module assert_mod           ! a minimal test harness: assert, count, report
  use kinds, only: dp
  implicit none
  private
  public :: assert_close, assert_true, report, n_fail
  integer :: n_run = 0, n_fail = 0
contains
  subroutine assert_close(name, got, want, tol)   ! for floating-point: compare within a tolerance
    character(*), intent(in) :: name
    real(dp),     intent(in) :: got, want, tol
    n_run = n_run + 1
    if (abs(got - want) <= tol) then
      print '(a,a)', 'PASS  ', name
    else
      n_fail = n_fail + 1
      print '(a,a,es13.6,a,es13.6)', 'FAIL  ', name // '  got=', got, ' want=', want
    end if
  end subroutine assert_close

  subroutine assert_true(name, cond)              ! for logical invariants
    character(*), intent(in) :: name
    logical,      intent(in) :: cond
    n_run = n_run + 1
    if (cond) then
      print '(a,a)', 'PASS  ', name
    else
      n_fail = n_fail + 1
      print '(a,a)', 'FAIL  ', name
    end if
  end subroutine assert_true

  subroutine report()
    print '(a,i0,a,i0,a)', '--- ', n_run - n_fail, ' / ', n_run, ' checks passed'
  end subroutine report
end module assert_mod

Notice the single most important design choice in that whole module, and the one beginners get wrong: assert_close compares with a tolerance, never with ==. Two real(dp) values that should be equal frequently differ in their last bit or two, because floating-point arithmetic rounds (Chapter 20); an equality test on reals is a test that will fail for reasons that have nothing to do with correctness. We return to exactly how tolerant to be in §37.2.

So what do you test? Here the numerical setting makes life genuinely harder than ordinary software, and it is worth naming the difficulty. To write a test you need to know the right answer in advance — you need a source of truth, sometimes called a test oracle. For a function that sorts a list, the oracle is obvious. For a function that approximates the solution of a partial differential equation, the "right answer" is itself something the code is only estimating — so where does the oracle come from? There are four reliable sources, and every good numerical test suite draws on them:

1. Exact solutions on special inputs. Most numerical methods are exact, not approximate, on some class of inputs, and those inputs make perfect tests. Your five-point Laplacian is exact for polynomials up to cubic (Chapter 24), so feeding it $u = x^2 + y^2$, whose true Laplacian is $2 + 2 = 4$ everywhere, must return exactly $4$ — a value you know without running anything. Feeding it a linear field must return exactly $0$. These are oracles you can compute in your head.
2. Invariants and conservation laws. Even when you cannot predict the exact output, you often know a property the output must have. Pure diffusion obeys a maximum principle: with a stable timestep, no interior point can become hotter than the hottest boundary or colder than the coldest — heat spreads out, it does not concentrate. So all(u <= 100.0) after any number of steps of the hot-plate problem is a test that needs no oracle for the individual values, only the physics.
3. Symmetry. A symmetric problem must produce a symmetric answer. The hot-top plate is left–right symmetric, so after stepping, column $j$ and column $n+1-j$ must match to round-off. Break the stencil and the symmetry usually breaks first and most visibly.
4. Convergence order. When the method is only approximate, its error still has a predictable shape — it shrinks at a known rate as you refine the grid. Verifying that rate is a test, and it is the subject of §37.2 and the reason Chapter 22 built you the halving experiment.

Here are the first two oracles as real unit tests, driving the actual laplacian and stable_dt from your solver. The complete self-contained program is code/example-01-hand-rolled-assert.f90; this is its heart:

program test_solver
  use kinds,       only: dp
  use heat_solver, only: laplacian, stable_dt
  use assert_mod,  only: assert_close, assert_true, report, n_fail
  implicit none
  real(dp) :: u(5,5), lap(5,5), dt, r
  integer  :: i, j

  ! Oracle 1 — exactness on a quadratic: laplacian(x^2 + y^2) == 4 everywhere.
  do j = 1, 5
    do i = 1, 5
      u(i,j) = (real(i-1,dp)*0.5_dp)**2 + (real(j-1,dp)*0.5_dp)**2   ! h = 0.5
    end do
  end do
  lap = laplacian(u, 0.5_dp, 0.5_dp)
  call assert_close('laplacian(x^2+y^2) at (2,2) == 4', lap(2,2), 4.0_dp, 1.0e-12_dp)
  call assert_close('laplacian(x^2+y^2) at (3,3) == 4', lap(3,3), 4.0_dp, 1.0e-12_dp)

  ! Oracle 1 again — exactness on a linear field: laplacian == 0 (a steady state).
  do j = 1, 5
    do i = 1, 5
      u(i,j) = 25.0_dp * real(j-1, dp)          ! linear ramp 0,25,50,75,100 across columns
    end do
  end do
  lap = laplacian(u, 1.0_dp, 1.0_dp)
  call assert_close('laplacian(linear) at (3,3) == 0', lap(3,3), 0.0_dp, 1.0e-12_dp)

  ! Oracle 2 — an invariant: stable_dt must return a CFL-safe step, r = alpha*dt/h^2 <= 1/4.
  dt = stable_dt(1.0_dp, 1.0_dp, 1.0_dp, safety=0.8_dp)     ! dx=dy=1, alpha=1 -> dt = 0.2
  call assert_close('stable_dt(1,1,1,0.8) == 0.2', dt, 0.2_dp, 1.0e-12_dp)
  r = 1.0_dp * dt / 1.0_dp**2
  call assert_true ('stable_dt gives r <= 1/4 (stable)', r <= 0.25_dp)

  call report()
  if (n_fail > 0) error stop 1        ! nonzero exit code: the run FAILED (Chapter 13)
end program test_solver

$ gfortran -std=f2018 -Wall example-01-hand-rolled-assert.f90 -o test01 && ./test01
PASS  laplacian(x^2+y^2) at (2,2) == 4
PASS  laplacian(x^2+y^2) at (3,3) == 4
PASS  laplacian(linear) at (3,3) == 0
PASS  stable_dt(1,1,1,0.8) == 0.2
PASS  stable_dt gives r <= 1/4 (stable)
--- 5 / 5 checks passed

Every one of those expected values was computed by hand, not by running the code. At $(2,2)$ with $h = 0.5$, the field $x^2 + y^2$ has centre $0.5$, $x$-neighbours $0.25$ and $1.25$, $y$-neighbours $0.25$ and $1.25$, so the stencil is $(0.25 - 1.0 + 1.25)/0.25 + (0.25 - 1.0 + 1.25)/0.25 = 2 + 2 = 4$. The linear ramp has second difference $(25 - 50 + 75) = 0$ in one direction and $(50 - 100 + 50) = 0$ in the other, so its Laplacian is $0$ exactly. And stable_dt with safety = 0.8 returns $0.8 / (2\cdot 1 \cdot(1 + 1)) = 0.8/4 = 0.2$, giving $r = 0.2 \le 0.25$. The suite confirms what you already proved on paper — which is the point: a test encodes a proof so a machine can re-check it on every future edit.

Look at the last line of the program, because it is what makes a test automated rather than informational. On any failure the program calls error stop 1, ending with a nonzero exit code. That exit code is the entire contract between your tests and the outside world: the shell, fpm test, and — in §37.3 — a continuous-integration server all learn "the tests failed" from that single nonzero integer, exactly the mechanism Chapter 13 introduced when it distinguished stop from error stop. A test that prints "FAIL" but exits zero is invisible to automation; the exit code is what gives it teeth.

💡 Intuition — the oracle problem is the hard part of numerical testing. In most software the hard part of a test is arranging the inputs; the expected output is obvious. In numerical software it is the reverse: arranging inputs is easy, but knowing the correct output — when the code exists precisely because the answer is hard to compute — is the whole challenge. That is why the four oracle sources above are worth memorizing. You almost never test a solver by comparing to "the right answer" (you do not have it); you test it against an exact special case, an invariant it must obey, a symmetry it must preserve, or a convergence rate it must hit. Master those four and you can test code whose output you could never predict directly.

pFUnit: a real unit-test framework

Hand-rolled asserts are perfect for learning and for small codes, but a mature project wants a framework — one that discovers tests automatically, runs them all even when one fails, produces structured output a CI server can parse, and gives you rich assertions (compare arrays, compare within a tolerance, check for NaN) without rewriting them each time. For Fortran the standard choice is pFUnit, which you met named in Chapter 16's tour of the ecosystem, now put to work.

Definition (pFUnit). pFUnit (the parallel Fortran unit-testing framework, developed at NASA) is a unit-testing framework for modern Fortran. You write tests in .pf files — ordinary Fortran decorated with annotations such as @test (this procedure is a test) and assertions such as @assertEqual — which a pFUnit preprocessor turns into standard .F90 source, compiles, links against the pFUnit library, and runs through a generated driver that reports how many tests passed. Its "parallel" heritage matters for scientific code: pFUnit understands MPI and can run tests across ranks, so you can unit-test the parallel routines of Chapter 34, not just serial ones.

The same two oracle tests, written as pFUnit, look like this. Because a .pf file is preprocessed — not compiled directly by gfortran — this is a snippet you read and adapt, not one you feed to the plain compiler:

! test_heat.pf  --  processed by the pFUnit preprocessor, then compiled and run
@test
subroutine test_laplacian_quadratic()
  use funit                                  ! pFUnit's assertions and machinery
  use kinds,       only: dp
  use heat_solver, only: laplacian
  implicit none
  real(dp) :: u(5,5), lap(5,5)
  integer  :: i, j
  do j = 1, 5
    do i = 1, 5
      u(i,j) = (real(i-1,dp)*0.5_dp)**2 + (real(j-1,dp)*0.5_dp)**2
    end do
  end do
  lap = laplacian(u, 0.5_dp, 0.5_dp)
  @assertEqual(4.0_dp, lap(2,2), tolerance=1.0e-12_dp)   ! same oracle, framework syntax
  @assertEqual(4.0_dp, lap(3,3), tolerance=1.0e-12_dp)
end subroutine test_laplacian_quadratic

@test
subroutine test_stable_dt_is_cfl_safe()
  use funit
  use kinds,       only: dp
  use heat_solver, only: stable_dt
  implicit none
  real(dp) :: dt
  dt = stable_dt(1.0_dp, 1.0_dp, 1.0_dp, safety=0.8_dp)
  @assertEqual(0.2_dp, dt, tolerance=1.0e-12_dp)
  @assertLessThanOrEqual(dt / 1.0_dp**2, 0.25_dp)        ! r = alpha*dt/h^2 <= 1/4
end subroutine test_stable_dt_is_cfl_safe

You run the suite with your build tool — fpm test if pFUnit is wired in as a test dependency — and it prints a summary such as Tests run: 2, Failures: 0, returning nonzero if any assertion fails. The virtues over the hand-rolled version are real: @assertEqual already knows how to compare with a tolerance and how to report a mismatch legibly; the framework runs all tests and tallies them rather than stopping at the first failure; and the output is in a format CI tools understand. The concept, though, is identical to the module you wrote by hand — which is exactly why we wrote that module first. A framework is a convenience, not a mystery.

🐍 Python Comparison. If you have used pytest, pFUnit will feel familiar: @test is pytest's test discovery, @assertEqual(a, b, tolerance=t) is numpy.testing.assert_allclose(a, b, atol=t), and both frameworks exist to turn a pile of assert statements into a runnable, self-tallying suite. The lesson transfers in both directions: whether you test the Python orchestration layer with pytest or the Fortran kernel with pFUnit, the discipline — exact special cases, invariants, tolerances not == — is the same. When you wrap the kernel for Python with f2py (Chapter 15), you can even test the Fortran through pytest, comparing its output to NumPy on the same inputs — a reference implementation is a fifth kind of oracle.

🔄 Check Your Understanding. 1. Why does assert_close compare abs(got - want) <= tol instead of got == want? 2. You want to test laplacian but you cannot predict its output on a general temperature field. Name two oracles — sources of a known-correct answer — that let you test it anyway. 3. What does error stop 1 at the end of the test program accomplish that a print '(a)', 'some tests failed' would not?

Answers

(1) Floating-point arithmetic rounds, so two reals that should be equal often differ in the last bit or two; an == test fails for reasons unrelated to correctness. A tolerance accepts the round-off while still catching a real error. (2) Any two of: an exact special case (the Laplacian is exact for polynomials up to cubic, so $x^2+y^2 \to 4$ and a linear field $\to 0$); an invariant (a linear field is a steady state, so its Laplacian is $0$); symmetry (a symmetric field gives a symmetric Laplacian). (3) It sets a nonzero exit code, which is how the shell, fpm test, and a CI server detect failure automatically. A print is invisible to automation; the exit code is the contract that lets a machine act on the result.

37.2 Regression Testing: Comparing Against a Known-Good Result

Unit tests check the pieces. But you also want to guard the whole — to catch the day someone changes the solver and the final field quietly shifts, even though every unit still passes, because the bug is in how they combine or in a subtle change of behaviour no single unit test anticipated. The guard for that is a regression test.

🔗 Connection. You have met regression testing before, in a specific guise: Chapter 18 used it as the safety net for modernizing legacy code — capture the old program's output, then insist the modernized version reproduces it, so you can refactor fearlessly. Here we generalize that idea from a one-time migration to a permanent, automated guardrail on a living code, and we get precise about how to compare two outputs.

Definition (regression test). A regression test captures the output a program is known to produce on a fixed input — a known-good result, often called a golden file or reference — and then, on every later version, reruns that input and checks the output still matches. Its job is to catch a regression: a change that breaks something that used to work. Where a unit test asks "is this piece correct in isolation?", a regression test asks "did anything I did not intend to change, change?" It is the test that lets you optimize, refactor, and parallelize a solver while sleeping at night, because the moment the answer moves, the suite tells you.

Your solver has a perfect golden result already: the $5\times 5$ plate, hot top edge, two steps, that you hand-computed in Chapter 24 and carried unchanged through Chapter 36's reorganization. Freeze it as the reference and compare:

real(dp), parameter :: golden(5,5) = reshape([ &
    100.0_dp, 100.0_dp, 100.0_dp, 100.0_dp, 100.0_dp, &
      0.0_dp,  28.0_dp,  32.0_dp,  28.0_dp,   0.0_dp, &
      0.0_dp,   4.0_dp,   4.0_dp,   4.0_dp,   0.0_dp, &
      0.0_dp,   0.0_dp,   0.0_dp,   0.0_dp,   0.0_dp, &
      0.0_dp,   0.0_dp,   0.0_dp,   0.0_dp,   0.0_dp], [5,5], order=[2,1])
! ... build the 5x5 field, hot top edge, step twice ...
maxdev = maxval(abs(field%u - golden))
call assert_true('field matches golden within 1e-9', maxdev < 1.0e-9_dp)

The complete program is code/example-02-regression-golden.f90, and it prints:

$ gfortran -std=f2018 -Wall example-02-regression-golden.f90 -o test02 && ./test02
PASS  field matches golden within 1e-9
max deviation from golden:  0.00000E+00
bit-for-bit identical to golden: T   [true on a standard IEEE build; see below]
--- regression: 1 / 1 checks passed

Two lines of that output carry the deepest idea in the section, so slow down on them. The maximum deviation from the golden field is exactly zero, and the field is bit-for-bit identical to the reference — on this build. That raises the question at the heart of every regression test: how equal is equal enough?

Definition (bit-for-bit; tolerance). Two floating-point results agree bit-for-bit when they are identical to the last bit — the same IEEE 754 pattern in memory, so that a == b is exactly true. They agree within a tolerance when they differ by no more than an allowed amount, abs(a - b) <= tol (an absolute tolerance) or abs(a - b) <= tol*abs(b) (a relative one). Bit-for-bit is the strictest possible standard and the most fragile; a tolerance is looser and far more robust. Choosing between them is the central judgement call of a numerical regression test.

When can you demand bit-for-bit, and when must you settle for a tolerance? The rule follows from a fact you learned in Chapter 20: floating-point addition is not associative, so the order in which operations happen can change the last bits of the result. Anything that changes that order changes the bits:

• A different compiler, or a different version of the same compiler, may reorder or fuse arithmetic differently.
• Different optimization flags — -O3, and especially -Ofast/-ffast-math, which explicitly permit the compiler to reassociate floating-point math (Chapter 30) — change results in the low bits.
• A different math or linear-algebra library (a different BLAS/LAPACK) computes the same quantity by a different sequence of operations.
• Parallel execution reorders a reduction: an OpenMP reduction(+:s) (Chapter 33) or an MPI mpi_allreduce (Chapter 34) sums partial results in an order that depends on how many processors you used — so the same code gives different last bits on 4 cores versus 16.

So bit-for-bit reproducibility is achievable only when you pin the entire build and execution: same compiler, same version, same flags, same libraries, same number of processors. Within that frozen configuration it is a wonderfully strong check — it catches a change of a single bit, which is exactly what you want when verifying that a refactor changed nothing. But demand it across compilers, or across processor counts, and it will fail constantly for reasons that are not bugs, drowning you in false alarms. That is why the portable regression test in the code above uses a tolerance as its real assertion, and treats the bit-for-bit line as informational. Our golden values happen to be bit-identical here because $28$, $32$, and $4$ are small integers that land exactly on representable doubles — but you should not rely on that, and a robust suite never does.

⚠️ Common Pitfall — the bit-for-bit test that cries wolf. A tempting regression test stores a full-precision golden field and asserts exact equality. It passes beautifully on your laptop and then fails on every colleague's machine, in CI on a different compiler, and the first time anyone runs it in parallel — none of which are bugs. The team learns to ignore the red X, and a real regression sails through unnoticed because the test that would have caught it is the boy who cried wolf. Reserve bit-for-bit for a pinned configuration where reproducing exact bits is itself the thing you are testing; for everything portable, compare within a physically meaningful tolerance, and set that tolerance from the problem, not from wishful thinking.

Regression against an analytical solution

The golden field above is a regression against a previous run — it guards against change, but if that first run was itself wrong, the test faithfully preserves the bug. The stronger test compares against an answer you know is correct independent of any run: an analytical solution of the equation you are solving. This is the difference between "the code still does what it did" and "the code does the right thing," and it is worth its own name.

Definition (verification). Verification is confirming that a code correctly solves the equations it claims to solve — that the numerics are right. The gold standard is comparison against an analytical solution: a case where the true solution of the governing equation is known in closed form, so the code's error can be measured exactly and its convergence rate checked. (Verification is distinct from validation — confirming the equations themselves describe reality — a distinction the Chapter 38 capstone draws sharply. Testing lives mostly in verification.)

The heat equation hands you a clean analytical case. At steady state the field stops changing, so $\partial u/\partial t = 0$ and the equation collapses to $\nabla^2 u = 0$ — Laplace's equation. For linear Dirichlet boundary data, the exact steady solution is simply the linear field itself, and because the five-point stencil is exact for linear functions, the discrete solver reproduces it with zero discretization error. That makes the analytical steady state a fixed point of step: apply the solver to the exact linear field and it comes back unchanged, to round-off. It is a verification test with an oracle you can compute by hand:

! The analytical steady state for linear edges is the linear field; ∇²u = 0 exactly,
! so `step` must leave it unchanged. Residual should be ~ 0 (a discrete fixed point).
do j = 1, 5
  do i = 1, 5
    field%u(i,j) = 25.0_dp * real(j-1, dp)     ! exact steady state: 0,25,50,75,100
  end do
end do
before = field%u
call step(field, alpha=1.0_dp, dt=0.2_dp)
residual = maxval(abs(field%u - before))
call assert_true('linear steady state is a fixed point', residual < 1.0e-12_dp)   ! PASS: residual = 0

For the transient — the plate actually warming over time — the heat equation also has an exact separable solution. On the unit square with all edges held at zero, the fundamental mode $u(x,y,t) = \sin(\pi x)\sin(\pi y)\,e^{-2\alpha\pi^2 t}$ satisfies the equation exactly: it keeps its spatial shape and decays in time at the rate $2\alpha\pi^2$. Seed the solver with $\sin(\pi x)\sin(\pi y)$ and it should track that exponential decay, with an error that shrinks as you refine the grid and the timestep — and how fast it shrinks is the verification. This is precisely the convergence experiment Chapter 22 taught you: halve the grid spacing, and a second-order-accurate solver's error should fall by a factor of four. Turning that experiment into an automated test — refine, measure the error against the analytical solution, assert the ratio is near four — is the most rigorous regression test a PDE solver can have, because it verifies not just an answer but the whole order of accuracy of the method. The full convergence study is the capstone's validation section (Chapter 38); the Project Checkpoint below wires in the exactly checkable pieces.

🐛 Find the Bug. A colleague adds this "regression test" to the suite and is proud that it always passes:

fortran call step(field, alpha, dt) call assert_true('solver ran', allocated(field%u))

Why is this test worthless as a regression test, and what would make it a real one?

Answer

It asserts only that field%u is still allocated after step — which is true whether the physics is right, wrong, or a no-op. A test that cannot fail when the answer is wrong tests nothing; it is a "cry wolf" in reverse, giving false confidence. A real regression test must check the values: compare field%u against a known-good golden field within a tolerance (as in example-02), or against an analytical solution (the steady-state fixed point, or the separable mode's decay). The rule of thumb: if you cannot describe an input that would make the assertion fail, the assertion is not testing anything.

🔄 Check Your Understanding. 1. Give two reasons the same solver source can produce results that differ in the last bit, making a bit-for-bit regression test fail without any bug. 2. What is the difference between a golden-file regression test and a verification test against an analytical solution — and which one can silently preserve a pre-existing bug? 3. Why is a linear temperature field an especially good analytical case for testing this particular solver?

Answers

(1) Any two of: a different compiler or compiler version; different optimization flags (especially -ffast-math/-Ofast, which reassociate arithmetic); a different math/BLAS library; a different number of processors, which reorders a parallel reduction. All change the order of non-associative floating-point operations. (2) A golden-file test checks the output still matches a previous run (guards against change); a verification test checks it matches a known-true analytical answer (guards against being wrong). The golden-file test faithfully preserves a bug that was present in the run it captured; the analytical test does not, because its oracle is independent of the code. (3) The five-point stencil is exact for linear (indeed cubic) fields, so the discrete solution has zero error — the linear field is both the exact analytical steady state and a discrete fixed point of step, giving an oracle with no tolerance ambiguity at all.

37.3 Continuous Integration: Every Change, Checked Automatically

A test suite that no one runs is worthless, and human beings forget to run tests — especially the ones that take a minute, especially when they are sure their change is trivial, especially right before the deadline when it matters most. The fix is to take the remembering out of human hands.

Definition (continuous integration). Continuous integration (CI) is the practice of automatically building a project and running its test suite on every change — every commit or every proposed merge — on a clean, neutral machine, so that a break is caught within minutes of being introduced rather than weeks later by a colleague or a reviewer. The "continuous" is the point: instead of a big, painful integration-and-testing effort once in a while, the code is built and tested all the time, so it is never far from a known-good state. A CI system watches your repository, and on each change it checks out the code, compiles it, runs the tests, and reports pass or fail — the green check or red X you have seen on open-source projects.

The most widely used CI service for open-source code is GitHub Actions, which runs a workflow you describe in a YAML file under .github/workflows/. Here is a real, complete workflow that builds and tests a Fortran project across three versions of gfortran on every push and pull request:

# .github/workflows/ci.yml  --  build and test the heat solver on every change
name: CI
on: [push, pull_request]

jobs:
  test:
    runs-on: ubuntu-latest
    strategy:
      matrix:
        gcc: [11, 12, 13]          # test across THREE compiler versions
    steps:
      - uses: actions/checkout@v4
      - name: Install gfortran
        run: sudo apt-get update && sudo apt-get install -y gfortran-${{ matrix.gcc }}
      - name: Install fpm
        uses: fortran-lang/setup-fpm@v5
      - name: Build
        run: fpm build --compiler gfortran-${{ matrix.gcc }} --flag "-std=f2018 -Wall"
      - name: Run tests
        run: fpm test --compiler gfortran-${{ matrix.gcc }}

Read it as a recipe. On every push or pull_request, GitHub spins up a fresh Ubuntu machine, and — because of the matrix — it does so three times in parallel, once for each of gfortran 11, 12, and 13. Each run checks out the code, installs that compiler, installs fpm, builds the project, and runs fpm test. If every test suite exits zero, the job is green; if any test calls error stop 1 — the nonzero exit code from §37.1 — fpm test propagates that failure, the job goes red, and (for a pull request) the merge is blocked until it is fixed. That single chain, from an assert_close failure all the way to a blocked merge, is the entire value of CI: a mistake cannot reach the main branch without someone seeing red.

The matrix across compilers is not decoration; it is the most valuable thing CI does for scientific Fortran. Code that compiles cleanly on your gfortran may use an extension the standard does not guarantee, or lean on behaviour that differs between compilers, or trip a bug fixed in a later version. Building on three gfortrans — and, in a serious project, on Intel's ifx and perhaps NVIDIA's nvfortran too — catches portability problems the day they are introduced, while the change is small and fresh in your mind, instead of six months later when a collaborator on a different cluster cannot build your code at all. This is the portability theme of Chapter 30 made into a reflex.

⚠️ Honesty note. The workflows in this chapter are shown as code you can read, adapt, and commit — but, in keeping with this book's rule, none was executed while writing it, so no green check or CI log is reproduced as certified output. The behaviour is what matters and it is reliable: a passing suite yields a green job, and a test that calls error stop with a nonzero code yields a red one. When you commit a workflow like this to a real repository, GitHub does the running; your job is to make the tests meaningful (§§37.1–37.2) and the exit codes honest.

⚡ Performance Note — tests must be fast, or they will be skipped. CI runs on every change, and a suite that takes twenty minutes is a suite people route around with "skip CI" commits until it might as well not exist. Keep unit tests to milliseconds by testing on small grids — the $5\times 5$ plate, not a $1000\times 1000$ one; the oracles of §37.1 are exact regardless of size, so a tiny grid tests the logic just as well and a thousand times faster. Reserve the expensive full-resolution convergence study for a separate, slower job that runs on a schedule or before a release, not on every keystroke. Fast tests are run tests; slow tests are dead tests.

🔄 Check Your Understanding. 1. What triggers a CI run in the workflow above, and what does the matrix: gcc: [11, 12, 13] line cause to happen? 2. Trace the chain: a unit test's assert_close fails. How does that become a blocked pull request? 3. Why is building across several compilers especially valuable for scientific Fortran, more so than for a program that will only ever run on your own laptop?

Answers

(1) Every push and pull_request triggers it; the matrix runs the whole build-and-test job three times in parallel, once each for gfortran 11, 12, and 13. (2) assert_close increments the failure count, so the test program calls error stop 1 and exits nonzero; fpm test sees the nonzero exit and itself fails; the CI step's command fails, so the job goes red; and a red required check blocks the pull request from merging. The exit code is the link at every stage. (3) Because scientific code is meant to run on other people's machines — clusters, collaborators' laptops, national-lab systems — with different compilers and versions. A compiler-matrix build catches non-portable code and compiler-specific bugs immediately, instead of when someone else cannot build your published code.

37.4 Documentation with FORD, and Version Control for Scientific Code

Correct code that no one can understand is only slightly better than incorrect code, and in six months "no one" will include you. Two disciplines keep a scientific code understandable over time: documentation that lives with the source and stays true to it, and version control that records the code's history and ties every result to an exact state of it.

Documentation with FORD. The best documentation is next to the code it describes, because documentation kept in a separate document drifts out of date the moment the code changes and no one updates both. So you write the documentation as specially-marked comments in the source, and a tool extracts them into browsable pages. For Fortran that tool is FORD — named, like pFUnit, in Chapter 16, now put to use.

Definition (FORD). FORD (FORtran Documentator) is a documentation generator for modern Fortran. It reads your source, understands its structure — modules, derived types, procedures, arguments, intents — and produces a cross-linked HTML site from special doc comments you write in the code. By convention FORD reads a comment beginning !> as documentation preceding the thing it describes, and one beginning !! as documentation following it (handy for annotating an argument on its own line). Because those are ordinary Fortran comments, FORD-documented source still compiles unchanged — the documentation and the code are the same file, so they cannot drift apart.

Here is the solver's laplacian, documented for FORD. To gfortran the !> and !! lines are just comments — it compiles exactly as before — but FORD turns them into a page describing what the function does, what each argument means, and how it relates to the rest of the code:

!> Compute the scaled five-point discrete Laplacian of a 2-D field.
!
!  Applies the second-order central stencil in each direction and divides by the
!  grid spacing squared, so the result approximates ∇²u to O(h²). Interior points
!  are filled; boundary rows and columns are returned as zero (they are held fixed
!  by the Dirichlet conditions and never enter the update). Exact for fields that
!  are polynomial of degree ≤ 3 — the basis of the unit tests in Chapter 37.
pure function laplacian(u, dx, dy) result(lap)
  real(dp), intent(in) :: u(:,:)   !! the temperature field, shape (nx, ny)
  real(dp), intent(in) :: dx       !! grid spacing in the first index (x)
  real(dp), intent(in) :: dy       !! grid spacing in the second index (y)
  real(dp) :: lap(size(u,1), size(u,2))  !! the discrete Laplacian, same shape as u
  ! ... body exactly as in Chapter 24 ...
end function laplacian

You generate the site by pointing FORD at a small project file (ford project.md) and it writes an HTML tree you can host or browse locally. The payoff compounds: a newcomer reads the interfaces of your code as a documented API rather than reverse-engineering them from the bodies, and because the docs are the comments, keeping them accurate is part of editing the code, not a separate chore that never happens.

Version control with git. You have surely used git; what changes for scientific code is what you commit and why it matters more. Git records the exact state of your code at every commit, each identified by a hash, and lets you return to any of them — which for science means something specific and powerful: a result can be tied to the exact commit that produced it. When you make a figure for a paper, you record the commit hash of the code that made it (git rev-parse HEAD), and now "the code that produced Figure 3" is not a vague memory but a precise, retrievable state you can check out and rerun years later. A few scientific-git habits are worth stating:

• Commit the code and the inputs. The namelist, the config, the small reference data a test needs — these are part of what makes a result, so they belong in version control beside the source. A run is code plus input.
• Tag releases and paper submissions. When you submit a paper, tag the commit (git tag paper-2026-submission); that tag is the permanent, named pointer to "the code as it was when we made these claims."
• Do not commit large binary outputs or datasets. Gigabytes of simulation output do not belong in git (they bloat the repository and git handles binaries poorly); record instead how to regenerate them — the input and the commit — which is smaller, more useful, and the essence of reproducibility.
• Write commit messages a scientist can use. "Fixed bug" tells a future reader nothing; "Fix sign error in y-Laplacian that broke the convergence test" tells them exactly what changed and why, and turns the git log into a lab notebook.

🔗 Connection — the build configuration belongs under version control too. Chapter 36 warned that the build configuration — compiler, flags, library versions — can change the numbers, and must be recorded. Version control is where you record it: commit the fpm.toml (which pins dependency versions to a git tag), commit the CI workflow (which names the exact compilers and flags), and the "how it was built" travels with the code automatically. A cloned repository at a given commit then carries not just the source but the recipe, which is most of what reproducibility requires — the subject of §37.5.

🔄 Check Your Understanding. 1. Why does writing documentation as !>/!! doc comments in the source keep it more accurate than a separate design document, and why does FORD-annotated code still compile? 2. What does it mean, concretely, to "tie a figure to a commit," and why is that more useful to a future reader than saving the output file? 3. Give two things you should commit for a scientific run and one thing you should not.

Answers

(1) The docs live in the same file as the code, so updating them is part of editing the code and they cannot drift out of sync the way a separate document does; and because !>/!! are ordinary Fortran comments, the compiler ignores them, so the documented source compiles unchanged. (2) It means recording the git commit hash of the code that produced the figure (e.g. git rev-parse HEAD), so anyone can check out that exact state and rerun it. That is more useful than the saved output because it lets a reader regenerate and vary the result — change a parameter, refine the grid — not merely look at a static file whose provenance is otherwise unknown. (3) Commit: the source, and the inputs (namelist/config, small reference data); tag paper submissions. Do not commit: large binary outputs or datasets — record how to regenerate them instead.

37.5 Reproducibility: The "Works on My Machine" Problem

Everything so far — tests, CI, docs, version control — serves one final goal, the one that separates computational science from computational guessing. A result someone else cannot regenerate is not evidence; it is an anecdote.

Definition (reproducibility). A computational result is reproducible when someone else — or you, later, on another machine — can regenerate it from the recorded materials: the same code, inputs, build, and environment producing the same result (bit-for-bit, or within a stated tolerance). Reproducibility is not a nicety bolted on at publication; it is the property that makes a computational result checkable, and therefore scientific. Its failure mode has a name every programmer knows — "it works on my machine" — which, said of a scientific result, is a confession: the result depends on undocumented features of one particular environment, so no one, including its author, can be sure it is real.

To reproduce a run you need everything that could change the answer, written down. Gather the list from the chapter:

The elegant move is to make the program itself record most of this, so the provenance is captured automatically at run time rather than remembered by a fallible human. Fortran helps directly: the intrinsics compiler_version() and compiler_options() (from iso_fortran_env) let a program report how it was built, and a version string and git hash baked in as parameters report what it is. A run then stamps its own provenance into its output:

program reproducibility_stamp
  use, intrinsic :: iso_fortran_env, only: compiler_version, compiler_options
  implicit none
  integer, parameter      :: dp = selected_real_kind(15, 307)
  character(*), parameter :: code_version = '1.0.0'
  character(*), parameter :: git_commit   = 'a1b2c3d'    ! injected by the build in practice
  integer, allocatable    :: seed(:)
  integer  :: nseed
  real(dp) :: x1, x2

  print '(a)', 'code version : ' // code_version         ! WHAT ran
  print '(a)', 'git commit   : ' // git_commit
  print '(a)', 'compiler     : ' // compiler_version()   ! HOW it was built (build-dependent)
  print '(a)', 'options      : ' // compiler_options()

  call random_seed(size=nseed)                            ! a FIXED seed makes the RNG reproducible
  allocate(seed(nseed))
  seed = 20260722                                         ! record THIS integer with the results
  call random_seed(put=seed)
  call random_number(x1)
  call random_seed(put=seed)                              ! re-seed identically...
  call random_number(x2)                                  ! ...and the draw repeats exactly
  print '(a,l1)', 'same seed reproduces the draw : ', (x1 == x2)
end program reproducibility_stamp

$ gfortran -std=f2018 -Wall -O2 example-03-reproducibility-stamp.f90 -o stamp && ./stamp
code version : 1.0.0
git commit   : a1b2c3d
compiler     : GCC version 13.2.0                         [build-dependent]
options      : -std=f2018 -Wall -O2 ...                   [build-dependent]
same seed reproduces the draw : T

Two of those lines are marked build-dependent, and honestly so: the whole purpose of compiler_version() and compiler_options() is to report the local build, which differs from machine to machine — that variability is the thing being recorded, so its exact text cannot be hand-computed and is shown representatively. The last line, though, is exact and it is the reproducibility lesson in miniature: seeding the generator with a fixed integer and drawing gives one value; re-seeding with the same integer and drawing gives the same value, so x1 == x2 is T. Record that seed with your results and the stochastic part of your run repeats.

There is a sharp honesty in that example worth stating, because it is a reproducibility trap. Fortran's standard does not fix which random-number generator random_number uses — so the same seed reproduces the same sequence on the same compiler and runtime, but a different compiler may implement a different generator and produce different numbers from the identical seed. Seeding buys you reproducibility within a build, not necessarily across builds; if you need cross-compiler reproducible randomness, you record and ship the actual random stream, or use a generator whose algorithm you control. This is the whole subject in one microcosm: reproducibility is never automatic, it is a set of things you deliberately record and pin.

🚪 Threshold Concept — a scientific result you cannot reproduce is not a result. This is the idea that reorganizes everything in this chapter, and much of your career. A number your code printed is not, by itself, science. It becomes science only when it is reproducible — when the code, the build, the inputs, and the environment are recorded well enough that another person can regenerate it and check it. Tests make the code right; version control makes it retrievable; recording the compiler, flags, libraries, and seeds makes the run repeatable; and only all of that together turns "here is a number I got" into "here is a result you can verify." Once you internalize that the deliverable of computational science is not the figure but the figure plus everything needed to regenerate it, you stop treating testing and provenance as chores and start treating them as what they are: the difference between doing science and merely running a program. The most important output of your solver is not the temperature field. It is a temperature field anyone can reproduce.

The point is not new, and its clearest statement predates most of the tooling. Computational scientists like to quote the argument, popularized by Jonathan Buckheit and David Donoho from Jon Claerbout's practice, that a published article about computational science is not the scholarship itself but merely advertising of it — the actual scholarship is the complete software and the complete set of instructions that generated the figures. Testing and reproducibility are how you make the scholarship real, and not just the advertisement.

🔗 Connection — this is where the whole book has been heading. Reproducibility ties back to Chapter 1's theme that scientific Fortran is infrastructure: the codes that predict weather and model the climate are trusted precisely because their results are checked, reproduced, and validated by many people over many years. It reaches forward to the Chapter 38 capstone, where you present your solver "as a paper would" — and a computational-science paper's credibility rests on exactly the reproducible artifact this chapter builds. The tested, documented, CI-guarded, provenance-stamped solver is not a tidier version of your code. It is the version a reviewer, a collaborator, or a future you can actually trust.

🔄 Check Your Understanding. 1. A colleague's paper reports a number their code no longer produces after a routine compiler upgrade. Name three things that, had they been recorded, would let you reproduce the original result — and say which single one ties the code itself down. 2. Why does a fixed random seed guarantee reproducibility on the same compiler but not necessarily across compilers? 3. In one sentence, what does "it works on my machine" reveal about a scientific result?

Answers

(1) Any three of: the exact code (a git commit hash — this is the one that pins the code itself), the compiler and version, the flags, the inputs, the library versions, the random seed, and (if parallel) the processor count. The git commit is what ties the code down. (2) Because the Fortran standard does not specify which random-number generator random_number uses, so a given seed reproduces the same sequence only for a given implementation; a different compiler may use a different generator and produce different numbers from the identical seed. (3) That the result depends on undocumented features of one particular environment — so it is not yet reproducible, and therefore not yet a checkable scientific result.

Project Checkpoint

Chapter 36 gave your solver a real source tree with an empty test/ directory waiting. This checkpoint fills it, and in doing so does to your project what a maintainer does to any code worth trusting: it adds tests, docs, and CI so the science inside can be changed without being broken. This is the increment Chapter 24 and Chapter 36 pointed at, and the one Chapter 38 presents as a finished artifact.

The test suite. Add three kinds of test to test/, each drawing on a different oracle:

• Unit tests of laplacian, step, and stable_dt, using the exactness and invariant oracles of §37.1: the Laplacian of $x^2+y^2$ is $4$ and of a linear field is $0$ (both exact); stable_dt returns a step with $r = \alpha\,dt/h^2 \le 1/4$.
• A regression test against the golden $5\times 5$, two-step field — the known-good result carried unchanged since Chapter 24 — compared within a tolerance, with a bit-for-bit check kept informational (§37.2).
• A verification test against the analytical solution: the linear steady state is an exact fixed point of step (residual $\sim 0$), and — cited from Chapter 22's order-of-accuracy experiment — the error against the separable mode falls at second order under grid refinement.

The complete, self-contained suite is code/project-checkpoint.f90: it bundles the solver modules, runs all the checks through a small assert harness, and — the crucial part — calls error stop 1 if any fail, so fpm test and CI can see the result. Here is its spine and its hand-computed output:

! In test/, driven by the same assert harness of section 37.1:
call run_unit_tests()          ! laplacian exactness (4 and 0); stable_dt CFL-safe
call run_regression_test()     ! 5x5 two-step field vs the golden reference (tolerance)
call run_verification_test()   ! linear steady state is a fixed point; max principle holds
call report()
if (n_fail > 0) error stop 1   ! any failure -> nonzero exit -> CI goes red

$ fpm test        # (standalone: gfortran -std=f2018 -Wall project-checkpoint.f90 && ./a.out)
PASS  unit: laplacian(x^2+y^2) == 4
PASS  unit: laplacian(linear)  == 0
PASS  unit: stable_dt gives r <= 1/4
PASS  regression: 5x5 two-step field matches golden (maxdev = 0.0)
PASS  verification: linear steady state is a fixed point (residual = 0.0)
PASS  verification: maximum principle holds (0 <= u <= 100)
--- 6 / 6 checks passed

Every value is hand-derived from the physics of Chapter 24: the golden field's warm cells are $28$, $32$, and $4$; the linear steady state has zero residual because the stencil is exact for linear fields; and the maximum principle holds because a CFL-safe step is a weighted average of a cell and its neighbours, which cannot exceed the boundary maximum of $100$.

The docs and CI. Annotate the public procedures of heat_solver with FORD !>/!! doc comments (§37.4) — they compile unchanged and generate an API page — and add the .github/workflows/ci.yml of §37.3 so every push builds and runs this suite across three gfortran versions. Commit the workflow and the fpm.toml alongside the source, and your solver now proves itself on every change.

That is the whole transformation. After Chapter 36 your solver was navigable; after this checkpoint it is trustworthy — it tells you the instant its answer moves, documents its own interfaces, checks itself on every compiler automatically, and can be reproduced from a commit hash. The capstone inherits software a reviewer would believe.

Summary

This chapter turned a navigable solver into a trustworthy one — tested, documented, continuously checked, and reproducible.

The two things to remember. First, test numerical code against oracles, not against "the right answer" — exact special cases, invariants, symmetry, and convergence order, compared within a tolerance. Second, a scientific result you cannot reproduce is not a result: the deliverable is the number plus everything needed to regenerate it — the tested code, the recorded build, the inputs, and the seeds.

Spaced Review

Retrieval practice on the two chapters this one builds directly on: error handling (Chapter 13), whose assertions and exit codes are the machinery of testing, and the anatomy of a real code (Chapter 36), whose package and build configuration are what we now test and reproduce. Answer before peeking.

1. (Ch. 13) A test program calls error stop 1 when a check fails. Explain the difference between stop and error stop, and why the nonzero exit code — not the printed "FAIL" message — is what makes the test usable by fpm test and CI.
   

   Answer

   `stop` is a normal termination (exit code $0$ by default); `error stop` is an *error* termination that sets a nonzero exit code (and in a parallel run halts every image at once). Automation — the shell, `fpm test`, a CI server — decides pass/fail from the process's *exit code*, not from any text it printed. A test that prints "FAIL" but exits $0$ looks successful to a machine; `error stop 1` makes the failure visible to automation, which is the entire chain from a failed assertion to a blocked merge.

2. (Ch. 13) Section 37.1's assert_close is an assertion with a tolerance. Recall Chapter 13's definition of an assertion and its "fail early, fail loudly" principle: how does a test suite apply that principle, and why is a tolerance (not ==) the right comparison for a floating-point assertion specifically?
   

   Answer

   An assertion checks a condition you believe is true and fails immediately and loudly if it is not, so a bug is caught at its source rather than producing corrupted output far downstream. A test suite applies this by running many assertions on known inputs before any real computation depends on them. For floating-point values a tolerance is required because arithmetic rounds ([Ch. 20](../../part-05-numerical-methods/chapter-20-floating-point/index.md)): two reals that should be equal often differ in the last bit, so `==` would fail on correct code, while `abs(got - want) <= tol` accepts round-off yet still catches a genuine error.

3. (Ch. 36) Chapter 36 warned that the build configuration can change the numbers and must be recorded. Connect that to §37.2: why can the same solver source fail a bit-for-bit regression test, and what part of the build configuration is most often responsible?
   

   Answer

   Because floating-point addition is not associative, so anything that reorders operations changes the last bits. The build configuration controls exactly those things: a different compiler or version, and above all different *optimization flags* — `-O3`, and especially `-Ofast`/`-ffast-math`, explicitly permit the compiler to reassociate arithmetic — as well as a different math/BLAS library. Change any of them and the same source produces bit-different (though equally valid) results, which is why a portable regression test compares within a tolerance and reserves bit-for-bit for a pinned build configuration.

4. (Ch. 36) Chapter 36 reorganized the solver into src/, app/, and an empty test/. Now that test/ is populated, explain how the five architectural roles (driver, solver, physics, I/O, utility) make the physics easy to unit-test — what property of a pure physics function like laplacian is exactly what a unit test wants?
   

   Answer

   A `pure` physics function has no side effects and no hidden state — its output depends only on its arguments — so a unit test can call it with known inputs and check the output in complete isolation, with no files, no global state to set up, and no order-dependence between tests. That is precisely the property the five-role separation was built to preserve: because `laplacian` is `pure` and takes its field and spacing as arguments (rather than reading a module variable or opening a file, which the I/O and driver roles keep away from it), it is trivially testable — the architecture of Chapter 36 is what makes the tests of Chapter 37 short.

What's Next

Your solver is now modular (Chapter 36), and tested, documented, and reproducible (this chapter). Every piece the book promised is in place: it computes real physics, it runs fast and in parallel, it writes visualizations, and it proves itself correct on every change. Chapter 38 is where it all comes together — the capstone, where you assemble the complete solver and present it as a computational-science paper: the problem and method, the validation against the analytical solution (the verification test of this chapter, run in full as a convergence study), the performance and scaling analysis, the ParaView figures, and the reproducible artifact behind them all. The testing and reproducibility discipline you built here is not a detour from the science; it is what makes the capstone science and not merely a program that ran. Let's write the paper.