Case Study 2: Building a Validated Steady-State Heat Solver

"Testing shows the presence, not the absence, of bugs — but a solver you have validated against a known answer is one you can build on." — after Edsger W. Dijkstra

Executive Summary

Case Study 1 fixed a broken run. This one builds — and, crucially, validates — a complete 2D steady-state heat solver on top of your project's field_t, step, and stable_dt. A solver that produces a plausible-looking heat map is not the same as a solver that is correct, and the difference is a validation you can defend. You will add a run_to_steady driver that marches until the field stops changing, then verify it three ways: against a 1D exact solution (a linear temperature profile), against a 2D manufactured exact solution (boundaries chosen so the answer is $u = x$ everywhere), and against two cheap physical invariants (the maximum principle and left–right symmetry). By the end you have not just a working solver but a trustworthy one — the foundation the Chapter 38 capstone's validation study is built on.

Skills applied

• Assembling the canonical solver pieces — field_t (Ch. 9), laplacian/step/stable_dt (§24.2–24.4).
• Writing a convergence-monitored time loop to steady state (§24.3).
• Imposing Dirichlet boundary conditions, including a linear ramp (§24.5).
• Validating a PDE solver against exact and manufactured solutions — the core discipline of Chapter 38.
• Reasoning about the optimization and parallelization path (Chapters 27, 29, 33–34).

Background

At steady state, temperatures stop changing: $\partial u/\partial t = 0$, so the heat equation collapses to Laplace's equation $\nabla^2 u = 0$. Physically, every interior point has reached the average of its neighbours. Our explicit time-marcher reaches this state naturally — run it long enough and the field stops moving — which gives us a simple, robust way to compute steady states and, more importantly, cases whose answers we know exactly and can check against.

The solver pieces are already canonical from the Project Checkpoint. The new piece is the driver that knows when to stop.

Phase 1 — Build the convergence-monitored driver

Marching to steady state means stepping until the field stops changing. We measure "change" as the maximum absolute difference between successive fields, and stop when it drops below a tolerance (or a step cap, so a mistake cannot loop forever — Chapter 13's defensive habit):

subroutine run_to_steady(field, alpha, dt, tol, max_steps, nsteps, final_change)
  type(field_t), intent(inout) :: field
  real(dp), intent(in)  :: alpha, dt, tol
  integer,  intent(in)  :: max_steps
  integer,  intent(out) :: nsteps
  real(dp), intent(out) :: final_change
  real(dp), allocatable :: prev(:,:)
  integer :: k
  final_change = huge(1.0_dp)
  do k = 1, max_steps
    prev = field%u                        ! snapshot before the step
    call step(field, alpha, dt)           ! one FTCS step (from heat_solver)
    final_change = maxval(abs(field%u - prev))
    if (final_change < tol) exit
  end do
  nsteps = k
end subroutine run_to_steady

The whole-array maxval(abs(field%u - prev)) is one line and reads like its intent — the array thinking of Chapter 5 paying off again. The routine returns both how many steps it took and the final change, so the caller can report convergence honestly rather than assume it.

Phase 2 — Validate against a 1D exact solution

The cleanest validation uses a case whose answer you know on paper. A 1D rod with ends held at $0$ and $100$ has the exact steady state of a straight line: $u(x) = 100x$ on the unit interval, i.e. the profile $0, 25, 50, 75, 100$ on five equally spaced nodes. This is exact because a linear function has zero second derivative, so $\nabla^2 u = 0$ everywhere — Laplace's equation is satisfied identically.

You already watched example-02-ftcs-1d.f90 march toward this: after three steps the interior read $1.5625,\ 12.5,\ 45.3125$, climbing toward $25, 50, 75$. Run it to convergence and it lands on the line. The validation is quantitative and hand-anchored: the solver's steady state must match $[0, 25, 50, 75, 100]$ to within your tolerance, and the discrete Laplacian of that profile must be zero (each interior node equals the average of its two neighbours: $25 = (0+50)/2$, $50 = (25+75)/2$, $75 = (50+100)/2$ — check them). If your solver converges to anything else, it is wrong, and you know it without a plot.

Sanity check (do it by hand). The steady interior node $u_3 = 50$ must satisfy the stencil exactly: $\nabla^2 u|_3 \propto u_2 - 2u_3 + u_4 = 25 - 100 + 75 = 0$. Zero, as Laplace demands. A single line of mental arithmetic certifies the answer.

Phase 3 — Validate against a 2D manufactured solution

One dimension is reassuring; two is the real test, and here we use a powerful trick: manufacture a problem whose exact 2D answer you know. Choose the boundary values so that the function $u(x, y) = x$ satisfies them — that is, set every boundary point equal to its own $x$-coordinate (left edge $x=0 \to 0$, right edge $x=1 \to 100$ if we scale by 100, and the top and bottom edges to a linear ramp across $x$). Then the exact steady solution is $u(x, y) = 100x$ everywhere, independent of $y$, because $\nabla^2 x = 0$.

Why this is a genuine test: $u = 100x$ is an exact fixed point of the discrete stencil. If the field is linear in the first index, its discrete second difference in that direction is zero, and it does not depend on the second index at all, so the five-point Laplacian is exactly zero at every interior point — step leaves it untouched. You can verify the fixed-point property on a tiny grid by hand: for a row of interior values $\dots, 25, 50, 75, \dots$ (spacing $25$), the stencil gives $25 - 2(50) + 75 = 0$ in $x$ and $50 - 2(50) + 50 = 0$ in $y$. Start the solver from a cold interior with these ramped boundaries, march to steady state, and it must converge to $u(i,j) = 100\,x_i$ — a plane tilted along $x$, flat along $y$.

! Manufactured exact solution u = 100*x : set boundaries to the linear ramp, interior cold.
call f%init(nx=n, ny=n, dx=1.0_dp/real(n-1,dp), dy=1.0_dp/real(n-1,dp))
do j = 1, n
  f%u(1, j) = 0.0_dp                                ! left edge  x=0
  f%u(n, j) = 100.0_dp                              ! right edge x=1
end do
do i = 1, n
  f%u(i, 1) = 100.0_dp * real(i-1,dp)/real(n-1,dp)  ! bottom edge: ramp in x
  f%u(i, n) = 100.0_dp * real(i-1,dp)/real(n-1,dp)  ! top edge:    ramp in x
end do
! interior stays 0 from init; march to steady state, then compare to 100*x_i

After convergence, compute the error against the exact solution, maxval(abs(f%u(i,j) - 100*x_i)). It should be at the level of your convergence tolerance, not the discretisation error — because for this solution the discretisation is exact. That is the strongest kind of validation: a case where the numerical answer should match the analytical one to (nearly) machine precision, so any real discrepancy is a genuine bug.

Why "manufactured"? The method of manufactured solutions is a standard verification technique in computational science: pick an answer, work out what boundary data (and source terms) make it exact, then check the code reproduces it. It decouples "is my code correct?" from "is my grid fine enough?" — here the linear solution removes discretisation error entirely, isolating coding errors. The Chapter 38 capstone uses exactly this idea, with a time-dependent exact solution, to certify the full solver.

Phase 4 — Cheap physical invariants as always-on checks

Exact solutions are gold but specific. Two invariants hold for every Dirichlet steady state and cost almost nothing to check, so leave them on as assertions:

The maximum principle is deep physics — heat cannot spontaneously make an interior point hotter than every wall — and a violated maximum principle is a red flag for a boundary bug or an instability that slipped past your CFL check. Symmetry catches indexing errors that a plausible-looking heat map hides: if you hold the top edge hot and the left/right edges equal, the steady field must be mirror-symmetric about the centre column, and the tiniest asymmetry means a stencil or boundary index is off. In the Project Checkpoint's two-step trace, you saw this symmetry appear immediately — row two came out $28, 32, 28$, symmetric about the centre. That was not luck; it was the physics, and it is a check you can automate.

Phase 5 — The path from correct to fast

You now have a correct solver, which is the only kind worth optimizing. The road ahead is the rest of the book, and it never touches the physics you just validated:

• See it: Chapter 26 writes VTK per step so you can watch the plate warm in ParaView.
• Measure it: Chapter 27 shows the stencil's loop order is a several-fold speed difference (column-major!), and Chapter 28 profiles where the time goes.
• Tune it: Chapter 29 blocks and vectorises the update.
• Scale it: Chapters 33–34 run the update across cores and cluster nodes — behind the same step interface, so your validation still holds.

Every one of those steps must reproduce the validated answer to within tolerance; validation is what makes optimization safe. That is the discipline the Chapter 38 capstone formalises.

Discussion Questions

1. Phase 3's manufactured solution is linear, so the discretisation error is exactly zero. Why is that a feature for finding coding bugs, but a limitation if you want to measure the solver's order of accuracy? What kind of exact solution would you manufacture to test convergence order instead?
2. run_to_steady stops when the max change drops below tol. Name two ways this can report "converged" when the solver has not reached the true steady state, and how you would guard against each.
3. The maximum principle is a physical law the continuous equation obeys. Does the explicit discrete scheme obey it for all stable $r$, or only some? (Hint: consider the update coefficients as weights, and when they are all non-negative.)
4. Steady state via time-marching can be slow (many steps). What faster route to $\nabla^2 u = 0$ does Chapter 21 offer, and what do you trade for the speed?

Your Turn: Extensions

• Option A (build). Add a convergence history: record final_change every $k$ steps and write it to a file. Plotting it (log scale) reveals the geometric decay rate — and a curve that flattens above tol warns you the solver has stalled, not converged.
• Option B (validate deeper). Manufacture a non-linear exact solution, e.g. $u = \sin(\pi x)\sinh(\pi y)/\sinh(\pi)$, which exactly satisfies $\nabla^2 u = 0$ with sinusoidal boundary data. Measure the solver's error versus $h$ and confirm the expected $O(h^2)$ convergence — the order-of-accuracy study Phase 3's linear case could not provide.
• Option C (design). Generalise run_to_steady to accept a boundary-condition procedure argument so the same driver runs Dirichlet, Neumann, or periodic problems by swapping one callback — the abstraction that turns a one-off solver into a small framework (Chapter 10's class machinery, if you want it typed).

Key Takeaways

• A solver is not correct until it is validated. Build the driver, then prove it against answers you know.
• Exact solutions (the 1D linear profile) and manufactured solutions (2D $u = 100x$) give you cases where the numerical answer should match the analytical one to near machine precision — the strongest tests, because any real discrepancy is a genuine bug, not discretisation error.
• Cheap invariants — the maximum principle and symmetry — cost nothing and catch boundary and indexing bugs a plausible heat map would hide; leave them on.
• Validation is what makes the later optimization and parallelization safe: every faster version must reproduce the validated answer, behind the frozen step interface.