Case Study 1: The Overnight Run That Exploded

"The first principle is that you must not fool yourself — and you are the easiest person to fool." — Richard Feynman

Executive Summary

A colleague hands you a 2D explicit heat-diffusion code that "worked perfectly last week" and now fills the terminal with NaN a few hundred steps into every run. The only change: they refined the grid from $26\times26$ to $51\times51$ points "to get a sharper picture," and left the timestep alone. Nothing looks wrong in the code — the stencil is correct, the boundaries are correct, the arithmetic is correct — yet it detonates. This case study is a guided diagnosis: you will reproduce the failure, instrument the code to expose the diffusion number $r = \alpha\Delta t/h^2$, watch it cross the $r = 1/4$ cliff exactly when the grid was refined, and apply the one durable fix — deriving the timestep from the grid. The bug is invisible in any single line and obvious the moment you compute one number.

Skills applied

• Reading a finite-difference stencil and confirming it is correct (§24.2).
• Recognising an explicit FTCS scheme and its stability limit (§24.3, §24.4).
• Computing the diffusion number and applying the CFL condition $r \le 1/4$ (§24.4).
• Instrumenting numerical code to expose a hidden invariant (§24.4, and -ffpe-trap from Chapter 13).
• Distinguishing a coding bug from a numerical bug — the code is right; the timestep is wrong.

Background

The code solves the 2D heat equation on a unit square with a hot top edge and cold sides (Dirichlet everywhere), using the explicit five-point stencil. It is, in fact, structurally identical to your own project's step. Here is the relevant core as your colleague wrote it:

integer,  parameter :: n = 51                     ! was 26 last week
real(dp), parameter :: alpha = 0.1_dp, dt = 0.002_dp
real(dp) :: u(n,n), u_new(n,n), h
integer  :: i, j, step

h = 1.0_dp / real(n - 1, dp)                      ! grid spacing on the unit square
! ... initialise: u = 0, u(1,:) = 100 (hot top edge) ...

do step = 1, nsteps
  u_new = u
  do j = 2, n-1
    do i = 2, n-1
      u_new(i,j) = u(i,j) + alpha*dt/h**2 *                     &
                   (u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1) - 4.0_dp*u(i,j))
    end do
  end do
  u = u_new
end do

The stencil is right. The snapshot update (into u_new) is right. The Dirichlet edges are held. And yet.

Phase 1 — Reproduce and confirm it is not a coding bug

First, establish that the code is correct by running the old configuration. Set n = 26, and the run is well-behaved: the hot edge warms the interior, the field marches smoothly toward a steady temperature gradient, no NaN. Flip n to 51 and it explodes. Same code, same physics, same boundaries — only the grid resolution changed. That single fact is the tell: a coding bug would fail at both resolutions. A bug that appears only when you refine the mesh is almost always a stability bug, not a logic bug.

Sanity check. Before hunting further, rule out the boring explanations. Compile with -fcheck=all -ffpe-trap=invalid,zero,overflow (Chapter 13). The trap fires on an overflow, not an out-of-bounds access — confirming the numbers are growing without limit, not that an index ran off the array. Overflow-to-Inf-then-NaN is the fingerprint of instability.

Phase 2 — Instrument the hidden invariant

The FTCS scheme has one number that decides its fate, and the code never prints it. Add two diagnostics: the diffusion number $r$, computed once, and the maximum interior magnitude, printed every so often.

real(dp) :: r
r = alpha*dt/h**2
print '(a, f8.4)', 'diffusion number r = ', r
print '(a, f8.4)', 'CFL limit (2D)     = ', 0.25_dp
! ... inside the time loop, occasionally:
!     print '(a,i6,a,es12.3)', 'step ', step, '  max|u| = ', maxval(abs(u))

Now compute $r$ by hand for both grids, because this is the whole case. On the unit square, $h = 1/(n-1)$:

There it is. Halving $h$ quadrupled $r$ (because $r \propto 1/h^2$), carrying it from a comfortable $0.125$ straight past the cliff to $0.5$ — exactly double the limit. The instrumented run confirms it: at $n=26$, max|u| sits at $100$ (the hot edge, as it should); at $n=51$, max|u| climbs geometrically — hundreds, then thousands, then Infinity, then NaN.

Phase 3 — Locate the cliff precisely

The colleague's instinct is to "use a smaller timestep until it works." That works, but blindly. The precise statement is: with $\alpha = 0.1$ and $h = 0.02$, stability requires

$$ r = \frac{\alpha\,\Delta t}{h^2} \le \frac14 \quad\Longrightarrow\quad \Delta t \le \frac{h^2}{4\alpha} = \frac{0.0004}{0.4} = 0.001. $$

Their $\Delta t = 0.002$ is exactly twice the largest stable step. Any $\Delta t \le 0.001$ is stable; anything above detonates. That is not a matter of taste or "enough" — it is a hard threshold you can compute to three digits before running anything.

Phase 4 — The durable fix: derive the timestep from the grid

The lesson is not "use $\Delta t = 0.001$." It is: never hard-code a timestep for an explicit scheme. Compute it from the grid every time, with a safety margin, so the code is correct at any resolution:

real(dp), parameter :: safety = 0.9_dp            ! stay 10% under the CFL limit
real(dp) :: dt
dt = safety * h**2 / (4.0_dp*alpha)               ! 2D limit, recomputed from h

With $n = 51$, $h = 0.02$: $\Delta t = 0.9 \times 0.0004/0.4 = 0.9 \times 0.001 = 0.0009$, giving $r = 0.1\times0.0009/0.0004 = 0.225 \le 0.25$. Stable, with a margin. Refine to $n=101$ and the same line gives $\Delta t = 0.000225$ automatically — no human in the loop, no overnight surprise. The cost, of course, is more steps: the finer grid needs roughly $4\times$ as many to reach the same physical time (the $\Delta t \sim h^2$ tax of §24.4), which is annoying but correct, and infinitely preferable to a screen of NaN.

Verify the fix on a tiny grid. Reduce to the hand-checkable $5\times5$, $h=1$, $\alpha=1$ case of the Project Checkpoint. There dt = 0.9 * 1/(4*1) = 0.225, giving $r = 0.225$; one step of the hot-top-edge plate warms the second row to $0 + 0.225\times100 = 22.5$, not the blow-up of the unstable run. Stable, predictable, and reproducible with a pencil.

Phase 5 — Prevent the recurrence

Fixing this run is easy; preventing the next one is the professional move. Add a guard at setup (Chapter 13's error stop) so the code refuses to run an unstable configuration rather than silently producing garbage:

if (alpha*dt/h**2 > 0.25_dp) then
  print '(a, f8.4, a)', 'FATAL: diffusion number r = ', alpha*dt/h**2, ' exceeds the 2D CFL limit 0.25'
  error stop 'unstable timestep — reduce dt or coarsen the grid'
end if

A run that would have wasted a night now dies in the first millisecond with a message that tells you exactly what to change. This is the difference between a script and a piece of scientific software.

Discussion Questions

1. The stencil, boundaries, and update were all correct, yet the code failed. Argue for why this should be called a numerical bug rather than a coding bug, and what that distinction implies for how you test numerical software.
2. Why is refining the grid — an apparently innocent "make the picture sharper" change — such a common trigger for instability in explicit diffusion codes specifically? Would a wave solver (CFL $\sim h$) have failed as dramatically on the same refinement?
3. The safety factor was $0.9$. What are the risks of setting it to $0.999$? To $0.5$? How would you choose it for a production run versus a quick test?
4. The error stop guard checks $r \le 1/4$. In a code with spatially varying $\alpha$ (a composite material), what is the correct thing to check, and against which $\alpha$?

Your Turn: Extensions

• Option A (analyse). Take the unstable $n=51$ run and, using the amplification factor $G = 1 - 8r$ for the worst mode, predict how many steps it takes for a round-off-seeded perturbation of size $\varepsilon_{\text{mach}} \approx 2\times10^{-16}$ to grow to order $100$. Compare with the observed "few hundred steps."
• Option B (instrument). Add a monitor that prints the diffusion number and aborts if maxval(abs(u)) ever exceeds, say, $10\times$ the maximum boundary temperature — a cheap runtime tripwire that catches instability within a handful of steps instead of overnight.
• Option C (fix forward). Replace the explicit step entirely with an implicit backward-Euler step (Exercise E17, via dgesv), which is unconditionally stable, and show it runs happily at the colleague's original $\Delta t = 0.002$ on the $51\times51$ grid — trading a linear solve per step for freedom from the CFL limit.

Key Takeaways

• A bug that appears only when you refine the grid is a stability bug, not a logic bug — suspect the CFL condition first.
• The diffusion number $r = \alpha\Delta t/h^2$ is the single invariant that decides an explicit heat run's fate; print it, and know it must be $\le 1/4$ in 2D.
• Because $r \propto 1/h^2$, halving the spacing quadruples $r$ — the most common way real explicit diffusion codes explode.
• Never hard-code an explicit timestep. Derive it from the grid (dt = safety * h**2/(4*alpha)) and guard the configuration with an error stop, so an unstable run dies in a millisecond with a clear message rather than an overnight screen of NaN.