Case Study 1: The NaN at Step 30,000

"The first step in fixing a bug is making it reproducible; the second is making it loud."

Executive Summary

A solver runs for hours, finishes without complaint, and writes an output field full of NaN. Nothing crashed; the exit code was zero; the program "worked." This is the worst kind of failure — silent, late, and disconnected from its cause — and it is exactly the kind this chapter's tools are built to defeat. In this study we take a heat-solver run that produces NaN output and hunt the bug to its source, using the development flags and debuggers of §§13.3–13.4 in the order a professional actually reaches for them. The lesson is a method: make the failure reproducible, make it loud, read where it becomes loud, then fix the cause and install a guard so it can never return silently.

Skills applied: guarding and diagnosing allocation and arithmetic (§13.1); -fcheck, -ffpe-trap, and -fbacktrace (§13.3); gdb inspection (§13.4); assertions and preconditions as prevention (§13.5); the uninitialized-value and precision suspects (§13.6).

Background

The code is a 2-D heat solver of the kind you are building. Each step it computes a dimensionless coefficient $r = \alpha\,\Delta t / \Delta x^2$ and applies the five-point update across the grid. On this run it read its parameters from a namelist a colleague edited, ran 100,000 steps, and produced a field of NaN. The team first suspects the physics — an unstable timestep, the CFL trouble of Chapter 24 — but the timestep looks conservative. Something else is wrong, and the output alone will not say what. We treat it as a detective problem.

Phase 1 — Reproduce, and Ask -fcheck First

A bug you cannot reproduce you cannot fix, so the first move is to shrink the run until the NaN appears quickly — the same 8×8 grid, only 50 steps — and confirm it still fails. It does: the field is NaN by the end. Reproducibility in hand, we reach for the cheapest instrument, -fcheck=all, on the theory that the bug might be an out-of-bounds write corrupting the field:

$ gfortran -std=f2018 -Wall -fcheck=all -g -O0 heat.f90 -o heat && ./heat

It runs clean — no bounds error — and still produces NaN. That is itself a clue: the bug is not an array-index fault. -fcheck catches out-of-bounds accesses, use of unallocated data, and bad pointers; its silence rules those out. The NaN is being computed, not corrupted into existence. We need the tool that watches arithmetic.

Phase 2 — Make It Loud with -ffpe-trap

By default, a NaN is born quietly and then contaminates everything it touches. §13.3's remedy is to trap the floating-point exception so the program halts at the instant the NaN is created:

$ gfortran -std=f2018 -Wall -ffpe-trap=invalid,zero,overflow -g -fbacktrace heat.f90 -o heat && ./heat

Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.

Backtrace for this error:
#0  ...
#1  ...  in coefficient_  at heat.f90:71
#2  ...  in step_         at heat.f90:52
#3  ...  in MAIN__        at heat.f90:29

That crash text is representative and version- and platform-dependent — we did not run the program to obtain it — but its structure is exactly what you will see: a SIGFPE, and a backtrace naming the call chain. The top frame points at line 71, inside a routine that computes the coefficient. The failure is no longer at "step 30,000, somewhere"; it is at one line. And note which exception fired: we trapped zero among others, which points at a division by zero.

The turn of the investigation: the NaN in the output was a symptom. The trap caught the disease — the first bad operation — thousands of steps earlier. This is the entire value of -ffpe-trap: it collapses the distance between cause and symptom to zero.

Phase 3 — Inspect the State with gdb

Line 71 computes r = alpha * dt / dx**2. A division-by-zero there means the denominator dx**2 is zero, which means dx is zero. But why would the grid spacing be zero? We stop the program just before the fault and look, using gdb (§13.4):

$ gdb ./heat
(gdb) break heat.f90:71
(gdb) run
Breakpoint 1, coefficient () at heat.f90:71
71          r = alpha * dt / dx**2
(gdb) print dx
$1 = 0
(gdb) print alpha
$2 = 9.9999999999999995e-05
(gdb) print dt
$3 = 0.20000000000000001

There it is: alpha and dt are the sensible values from the config, but dx is 0. The transcript is representative — your addresses and float formatting will differ — but the finding is the point: two of the three inputs are right and one is zero. The bug is upstream of the arithmetic; the division merely reported it. We now ask where dx was supposed to be set.

Phase 4 — Root Cause: An Unset, Unvalidated Parameter

Reading the setup code, the story assembles. The grid spacing dx is meant to be computed from the plate size and grid count, but the refactor that introduced the namelist left dx to a default of 0.0 and the line that should have set it was guarded by a branch that this config did not take. Nothing ever gave dx a real value. It was, in effect, an uninitialized parameter (§13.6) that happened to be zero — and because nothing validated it, the zero flowed straight into a division.

The causal chain, start to finish:

Step What happened Which suspect
Setup dx left at its 0.0 default; the setting line skipped uninitialized/unset value (§13.6)
Line 71 alpha*dt / dx**2 divides by zero → Inf no precondition check (§13.5)
Next op Inf combined with a finite update → NaN silent propagation (§13.3)
Step 30,000 NaN has spread across the whole field —
Output a field of NaN, exit code 0 no postcondition check (§13.5)

The physics was never the problem. The problem was a value that was never set and never checked — precisely the failure the Project Checkpoint's validate_config exists to prevent.

Phase 5 — Fix the Cause, Install the Guard

Fixing the immediate bug (compute dx correctly) is necessary but not sufficient; the professional move is to make this class of failure impossible to ship silently again. Two guards do it.

First, a precondition on the configuration, exactly the checkpoint's pattern — reject a non-positive dx before any arithmetic:

if (dx <= 0.0_dp) then
   write(error_unit, '(a, es10.2)') 'config error: dx must be positive, got ', dx
   error stop 2
end if

Second, a postcondition that no step is allowed to produce a NaN, so even a future bug surfaces at its first step instead of in the final file. A NaN is the only value not equal to itself, which gives a one-line test:

program coefficient_check
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  real(dp) :: alpha, dt, dx, r
  alpha = 1.0e-4_dp;  dt = 0.25_dp;  dx = 0.25_dp     ! now dx is a real value

  r = alpha * dt / dx**2
  if (r /= r) error stop 'coefficient is NaN'          ! x /= x  is true only for NaN
  print '(a, es9.2)', 'r = ', r
end program coefficient_check

$ gfortran -std=f2018 -Wall -O2 coefficient_check.f90 -o coeff && ./coeff
r =  4.00E-04

With a real dx = 0.25, the coefficient is $r = (10^{-4} \times 0.25) / 0.25^2 = 2.5\times10^{-5} / 0.0625 = 4.0\times10^{-4}$ — a finite, hand-checkable value, so the NaN guard stays silent and the program prints it. Had dx been zero, the precondition would have stopped the run with a clear message and exit code 2, long before any NaN could form. The bug that cost an afternoon becomes a five-second, self-explaining halt.

Discussion Questions

  1. -fcheck=all ran clean on the buggy program, yet the program was clearly broken. Explain precisely what that silence did and did not tell us, and why it was still useful information.
  2. The team's first hypothesis was a physics problem (an unstable timestep). How did the tools let them falsify that hypothesis quickly, rather than arguing about it? What is the general lesson about debugging by instrument versus debugging by opinion?
  3. The postcondition if (r /= r) error stop catches a NaN after it forms; -ffpe-trap catches it as it forms. When would you want each, given that -ffpe-trap is a build flag and the postcondition is in the source?

Your Turn: Extensions

  • Option A. Take any small numerical program of your own and deliberately introduce a divide-by-zero deep inside it. Find it twice: once by reading the NaN back to its source by hand, and once with -ffpe-trap=zero -fbacktrace. Time both. The ratio is the value of the flag.
  • Option B. Write a reusable elemental logical function is_nan(x) using the x /= x trick, and a companion is_finite. Compare them to the standard ieee_is_nan from the ieee_arithmetic module (Chapter 20). When would you prefer the intrinsic module?
  • Option C. Add the dx > 0 precondition to your solver's validate_config, then write a one-paragraph "incident report" for this bug: root cause, how it reached output, and the two guards that now prevent it. This is the habit that turns a painful afternoon into institutional memory.

Key Takeaways

  • A silent, late failure (NaN in the output, exit code 0) is defeated by a method: reproduce it, make it loud with -ffpe-trap, read the backtrace, inspect with gdb, then fix the cause and install a guard.
  • The tools falsify hypotheses fast. -fcheck's silence ruled out an index bug; -ffpe-trap located the arithmetic; gdb named the zero variable. Each step replaced an opinion with a fact.
  • The root cause was mundane and typical — an unset, unvalidated parameter — which is why validation at the boundary (§13.5, the Project Checkpoint) is worth more than any single debugging session: it stops the bug from ever forming.
  • x /= x is .true. only for NaN; it is the one-line postcondition that keeps a solver from writing poison to disk.