Case Study 2: Writing the Paper — Assembling Your Solver into a Reproducible Result

"If it disagrees with experiment, it is wrong. In that simple statement is the key to science." — Richard P. Feynman

Executive Summary

Case Study 1 taught you to judge someone else's result. Now you build your own. In this case study you take your finished heat solver and produce the whole deliverable: the assembled program, the verification study that proves it correct, the figures, the paper, and the reproducibility package that lets a stranger regenerate every number. This is the capstone of the capstone — the moment your forty-chapter project stops being code you wrote and becomes a result you can defend. You will confirm the assembled solver, generate the convergence study with a compilable harness (recovering second-order accuracy), lay out the figures, draft the seven-section paper, and assemble the README + build configuration + test that make it reproducible. Where Case Study 1 refereed, this one authors — the harder and more lasting skill.

Skills applied: assembling the modular solver behind one interface (§38.1, Chapter 36); building a convergence study against the analytical solution (§38.3); structuring results and figures (§38.5, Chapter 26); writing the paper and its reviewer-facing sections (§38.6); packaging for reproducibility (Chapter 37).

Background

Your solver is, by now, six modules behind the frozen step interface: kinds, timers, heat_types, heat_solver, heat_io, and the heat driver, plus a test/ directory. It runs, it parallelizes, it writes VTK. What it does not yet have is the argument that its numbers are right and the package that lets someone else reproduce them. We build both, in the order a paper is actually assembled: confirm the code, verify it, show the results, write it up, package it.

Phase 1 — Confirm the Assembled Solver

Before writing a word of the paper, re-run the one result you have hand-verified since Chapter 24 — the $5\times5$ plate, two steps — because a paper's implementation section stands on a solver you know is assembled correctly. The complete code/project-checkpoint.f90 produces it:

$ gfortran -std=f2018 -Wall -O2 project-checkpoint.f90 -o heat && ./heat
[1] Demonstration: the Chapter 1 plate (hot top edge)
    field after 2 steps:
      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
    max temperature:  100.00

That 28/32/28 and 4/4/4 is the same field you first computed by hand in Chapter 24, reproduced now by the fully assembled, parallel-ready solver. It is the anchor the whole paper rests on: the physics is right on the case you can check by pencil, so the machinery underneath is sound.

Phase 2 — Build the Verification Study

Now the paper's most important artifact. We verify against the analytical solution $u = \sin(\pi x)\sin(\pi y)e^{-2\alpha\pi^2 t}$, using the fact (§38.3) that the discrete mode is an exact eigenvector of the stencil, so the numerical field is exactly $G^K$ times the initial field — which lets us compute the whole-solution max error without a full march:

! cs02-convergence.f90 -- the full-solution convergence study for the write-up.
! The analytical mode is a stencil eigenvector, so the numerical field is exactly
! G^K times the initial field; the max error at time T is |G^K - exp(-2 alpha pi^2 T)|,
! with G = 1 - 8 r sin^2(pi h/2) and K = T/dt. Refining h -> error ~ h^2 (order 2).
! Compile: gfortran -std=f2018 -Wall cs02-convergence.f90 -o conv && ./conv
program cs02_convergence
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  real(dp), parameter :: pi = 3.141592653589793_dp
  real(dp), parameter :: alpha = 1.0_dp, r = 0.2_dp, tfinal = 0.05_dp
  real(dp) :: exact, h, dt, g, gk, err, prev
  integer  :: j, nsteps

  exact = exp(-2.0_dp*alpha*pi**2 * tfinal)
  print '(a,f8.6)', 'exact amplitude at T = 0.05 : ', exact
  print '(a)',      '     h        dt      steps   max error     ratio'
  prev = 0.0_dp
  do j = 1, 4
    h      = 0.5_dp**j                                  ! 1/2, 1/4, 1/8, 1/16
    dt     = r*h**2/alpha
    nsteps = nint(tfinal/dt)
    g      = 1.0_dp - 8.0_dp*r*sin(pi*h/2.0_dp)**2      ! amplification factor
    gk     = g**nsteps
    err    = abs(gk - exact)
    if (j == 1) then
      print '(f8.4,f11.6,i7,f13.6,a)',   h, dt, nsteps, err, '      --'
    else
      print '(f8.4,f11.6,i7,f13.6,f9.3)', h, dt, nsteps, err, prev/err
    end if
    prev = err
  end do
  print '(a)', 'error ratio -> 4 => second-order accuracy (observed order -> 2)'
end program cs02_convergence

$ gfortran -std=f2018 -Wall cs02-convergence.f90 -o conv && ./conv
exact amplitude at T = 0.05 : 0.372708
     h        dt      steps   max error     ratio
  0.5000   0.050000      1     0.172708      --
  0.2500   0.012500      4     0.028990    5.957
  0.1250   0.003125     16     0.006762    4.287
  0.0625   0.000781     64     0.001663    4.066
error ratio -> 4 => second-order accuracy (observed order -> 2)

There it is — the single most persuasive table in the paper. The error falls by a factor approaching four each time the grid is halved, so the observed order converges to two, exactly what the five-point stencil's theory (Chapter 22) demands. (The coarsest grid is pre-asymptotic — one interior point — so its ratio overshoots; the trend is what counts.)

Sanity check. Every number is exact, computed from the closed-form $G^K$, so your solver marched step-by-step on these grids must reproduce this table to rounding. That agreement is the verification; the day it fails, you have a bug, and the Chapter 37 regression test — which pins one row of this table — will tell you the same day you introduce it.

Phase 3 — Lay Out the Results and Figures

A results section needs the evidence figure and the result figure. You produce two:

The convergence plot is the one a reviewer trusts: four points falling on a straight line of slope 2, with a reference triangle beside them. The steady-state heat map is the one a reader enjoys: the hot top edge, the smooth gradient falling to the three cold edges. Both come from the division of labor the book has taught — Fortran and Python are better together: Fortran computes and writes the data, ParaView and matplotlib turn it into the pictures. You never plot from Fortran, and you never compute the kernel in Python.

Phase 4 — Draft the Paper

With the solver confirmed, verified, and figured, the paper writes quickly, because you assemble it from parts you already have:

Title:        A Verified Parallel Finite-Difference Heat-Equation Solver in Modern Fortran
Abstract:     what + method + key verification result + parallelism + honest performance  (5 sentences)
1. Intro:     the plate problem; why diffusion solvers matter; the contribution
2. Method:    heat eq, five-point stencil (O(h^2)), FTCS, CFL-safe dt          -> from Sec 38.2
3. V&V:       analytical solution + the Phase-2 convergence table (order 2)    -> from Sec 38.3
4. Impl/perf: the src/ architecture; Amdahl ceiling; memory-bound roofline     -> from Sec 38.1, 38.4
5. Results:   transient + steady state; the two figures                        -> from Phase 3
6. Conclusion:a verified, scalable solver; the explicit dt~h^2 limitation
7. Repro:     build config, inputs, the test suite, code availability          -> from Phase 5

The discipline of the abstract is the discipline of the whole paper: claim only what the evidence supports. "Verified against the analytical solution with observed second-order accuracy" is a sentence you can stand behind; "physically correct and blazingly fast" is not.

Phase 5 — Package for Reproducibility

The last artifact is the one that makes the result travel. A reproducibility package is a README, a pinned build configuration, the inputs, and the test suite:

README.md      : what it solves; `fpm build && fpm run`; expected first lines of output
fpm.toml       : the manifest (records name, version, dependencies)
build config   : gfortran -std=f2018 -O3 -march=native  (release);  add -fopenmp for threads;
                 mpif90 for the MPI build. RECORD THE FLAGS -- they can change the numbers.
inputs         : the namelist (grid, alpha, steps, save_every) that defines the run in the figures
test/          : test_solver.f90 -- regresses one row of the Phase-2 table against the analytical value
LICENSE        : so others may legally build on it (Ch. 40)

This is the Chapter 37 reproducibility discipline as a deliverable: someone clones the repository, reads the README, runs fpm build && fpm run, and gets your numbers — including the convergence table, because the test suite regenerates it. A result that reproduces is a result that can be built upon; a result that does not is an anecdote.

Discussion Questions

1. Phase 2 computes the convergence table from the closed-form $G^K$ rather than by marching the solver. Why is that legitimate for the write-up's purpose, and why must you also run the actual solver on these grids before claiming your code passes the study?
2. The reproducibility package records the compiler flags. Give a concrete example (recall Chapter 30) where the same source built with two different flag sets produces different numbers, and say why that makes flags part of the result.
3. The abstract is drafted last, though it is read first. Why is "write the abstract last" good advice, and what goes wrong when a paper's abstract is written before its verification?

Your Turn: Extensions

• Option A. Write the actual test/test_solver.f90 that marches the real solver on the $9\times9$ grid to $T = 0.05$ and asserts the max error is within $10\%$ of the Phase-2 value ($6.76\times10^{-3}$), printing PASS/FAIL. You have now connected the verification study to a runnable regression test.
• Option B. Add the $33\times33$ ($h = 1/32$) row to cs02-convergence.f90 and confirm the ratio stays near 4 and the order near 2 deeper into the asymptotic regime. Then deliberately break the scheme (use $r = 0.3 > 1/4$) and watch the "convergence" become a blow-up.
• Option C. Write the full one-page README.md for your solver — build commands, the expected first lines of output, and the one-paragraph abstract — then hand it to someone who has never seen the code and have them build and run it from the README alone. If they cannot, the package is not yet reproducible.

Key Takeaways

• A finished result is code plus the argument that it is right: assemble the solver, verify it against a known solution, produce the figures, write the paper, and package it so it reproduces.
• The convergence study is the paper's spine — the error falling like $h^2$ (observed order → 2) is the evidence that turns "trust me" into "here is why," and it doubles as the regression oracle.
• Build the two figures with the right tools: the convergence plot (matplotlib, the reviewer's figure) and the steady-state heat map (ParaView, the reader's figure) — Fortran computes, Python and ParaView draw.
• Claim only what the evidence supports. "Verified, second-order, memory-bound, Amdahl-limited" is defensible; "correct and blazingly fast" is not. Write the abstract last, to the evidence you actually have.
• A reproducibility package (README, build configuration, inputs, tests, license) is what makes the result travel — the difference between science others can build on and a number only you ever saw.