Case Study 2: Building the Numerical Foundation

"Choose your precision on purpose, or the hardware will choose it for you — and it will not choose well."

Executive Summary

Where the first case study repaired someone else's arithmetic, this one asks you to design the numerical bedrock of the heat solver from scratch and get it right the first time. Every serious scientific code begins with two decisions that are almost never revisited and almost always regretted if made carelessly: what precision do we compute in? and how do we encode the physical constants so they are correct, portable, and tamper-proof? You will build the kinds module the whole book depends on, declare the plate's physical constants as parameters with their units, derive the grid quantities without falling into the integer-division trap, and write a small validation program that prints the foundation back so you can eyeball it for physical sense. The output is not glamorous — it is nine lines of numbers — but those nine lines are the load-bearing footing that a parallel, validated simulation will eventually stand on.

Skills applied: requesting portable precision with selected_real_kind (§3.2); designing a single-purpose module with public/private; declaring named constants with parameter (§3.6); building a chain of derived parameters through correct mixed-mode arithmetic (§3.3–§3.4); formatted reporting (§3.7). Project advanced: this is the Chapter 3 project increment, built and validated end to end.

Background

Recall the physical problem from Chapter 1: a square metal plate, held at fixed temperatures on its edges, whose interior temperature evolves over time. To simulate it we lay a grid over the plate and track the temperature at each grid point. Before any physics can happen, the code needs to know a few numbers: how quickly the material conducts heat (the thermal diffusivity $\alpha$), how big the plate is, how fine the grid is, and — derived from those — how far apart the grid points sit. Get these constants and their precision right, once, in a place the whole program shares, and every later chapter builds on solid ground. Get them wrong — hardcode a diffusivity in single precision here, recompute a grid spacing with integer division there — and you will chase phantom bugs through code that is arithmetically sound but numerically poisoned at the root.

Phase 1 — Decide the Precision, and Own It

The first decision is precision, and we make it explicitly, in one place, so that no other file ever has to think about it again. Create kinds.f90:

module kinds
  implicit none
  private
  public :: dp
  integer, parameter :: dp = selected_real_kind(15, 307)
end module kinds

Three design choices are packed into these five lines. First, selected_real_kind(15, 307) states the scientific requirement — at least fifteen significant digits and a range to $10^{307}$ — rather than a representation like real(8); if this code is ever compiled on hardware where "kind 8" means something different, our request still resolves to a type that meets the requirement, or fails loudly at compile time. Second, the module is private by default and exports only dp through the explicit public statement, so it has exactly one job and cannot accidentally leak anything else. Third, it is a module, so any file can use kinds, only: dp and share the identical definition — no copy-paste, no drift. (You will study modules properly in Chapter 8; here we borrow just this one pattern.)

The design principle: a foundational decision like precision should be expressed once, as a requirement, in a shared location. The alternative — repeating selected_real_kind(15, 307) at the top of every program, as we did earlier in the chapter — works, but it scatters a decision that should live in a single sentence, and scattered decisions drift apart over a large codebase.

Phase 2 — Declare the Physical Constants as Parameters

With precision settled, encode the physics. These are values fixed for the run, so every one of them is a parameter — a compile-time constant the code cannot accidentally overwrite (§3.6):

real(dp), parameter :: alpha  = 1.0e-4_dp   ! thermal diffusivity, m^2/s
real(dp), parameter :: length = 1.0_dp      ! plate side length, m
integer,  parameter :: nx     = 101         ! grid points along a side

Note what each declaration carries beyond its number. The _dp suffix on 1.0e-4_dp guarantees the literal is formed in double precision, not rounded to single first — the trap of Exercise 3.12. The unit in the comment (m^2/s, m) is not decoration: a diffusivity is meaningless without its units, and half of all scientific-computing disasters trace to a quantity that was correct in the wrong unit. And nx is an integer, because a count of grid points is a whole number — using real for it would be both dishonest and an invitation to the very bugs this chapter warns against.

💡 Intuition: think of the parameter attribute as soldering a value into the circuit. A plain variable is a value on a breadboard — convenient, and one stray hand rearranges it. A parameter is soldered: fixed, cheap (folded into the machine code at compile time, §3.6), and impossible to knock loose by accident fifty lines later.

Phase 3 — Derive the Grid Quantities Without the Trap

Now the derived quantities — and this is exactly where a careless line quietly destroys everything. The grid spacing dx is the plate length divided by the number of intervals, which for nx points is nx - 1:

real(dp), parameter :: dx = length / real(nx - 1, dp)   ! grid spacing, m

Read that denominator carefully. nx - 1 is 100, an integer; if we wrote length / (nx - 1) the compiler would see real / integer, promote the integer, and — this particular case would actually work, because length is already real. But the habit that never fails is to convert the integer to real(dp) explicitly with real(nx - 1, dp), so the reader can see at a glance that the division is real-valued and no future edit (say, changing length to an integer plate-count) can silently turn it into integer division. Defensive arithmetic means making the real-ness visible, not relying on one operand happening to be real today.

From dx we can derive more, each a parameter because each is a compile-time function of constants:

real(dp), parameter :: tau    = length**2 / alpha        ! diffusion timescale, s
real(dp), parameter :: dt_max = 0.25_dp * dx**2 / alpha  ! stability-limited step, s

tau is a rough diffusion timescale — the order-of-magnitude time for heat to cross the plate — and dt_max is the largest timestep an explicit scheme can safely take. Why the safe step scales with $\mathrm{d}x^2 / \alpha$ is a genuine piece of numerical analysis that this chapter does not owe you; it is the stability story owned by Chapter 24. For now we are only encoding the arithmetic, correctly and in double precision, so that the reasoning has solid numbers to stand on when it arrives.

Phase 4 — Validate by Reporting the Foundation

A foundation you have not inspected is a foundation you are trusting on faith. Write a small driver that uses the module, assembles the constants, and prints them — a "does this look physical?" report you run once and read with your eyes:

program foundation
  use kinds, only: dp
  implicit none

  real(dp), parameter :: alpha  = 1.0e-4_dp
  real(dp), parameter :: length = 1.0_dp
  integer,  parameter :: nx     = 101
  real(dp), parameter :: dx     = length / real(nx - 1, dp)
  real(dp), parameter :: tau    = length**2 / alpha
  real(dp), parameter :: dt_max = 0.25_dp * dx**2 / alpha

  print '(a, i0)',     'precision (digits)  = ', precision(alpha)
  print '(a, i0)',     'range (exponent)    = ', range(alpha)
  print '(a, es10.3)', 'alpha  (m^2/s)      = ', alpha
  print '(a, f6.3)',   'length (m)          = ', length
  print '(a, i0)',     'nx (points)         = ', nx
  print '(a, f8.5)',   'dx     (m)          = ', dx
  print '(a, es10.3)', 'dx^2   (m^2)        = ', dx**2
  print '(a, es10.3)', 'tau = L^2/alpha (s) = ', tau
  print '(a, f6.3)',   'dt_max (s)          = ', dt_max
end program foundation

$ gfortran -std=f2018 -Wall foundation.f90 -o foundation && ./foundation
precision (digits)  = 15
range (exponent)    = 307
alpha  (m^2/s)      =  1.000E-04
length (m)          =  1.000
nx (points)         = 101
dx     (m)          =  0.01000
dx^2   (m^2)        =  1.000E-04
tau = L^2/alpha (s) =  1.000E+04
dt_max (s)          =  0.250

Read the report as a physicist, not a programmer. The precision line confirms we are computing in fifteen digits, as requested — the foundation is double precision, not some default we forgot to check. With nx = 101 points there are 100 intervals across a one-metre plate, so dx = 0.01 m; the program agrees. Squaring, dx^2 = 1.0 \times 10^{-4} m². The diffusion timescale tau = length^2 / alpha = 1 / 10^{-4} = 10^4 seconds — about three hours for heat to diffuse across the plate, which is physically sensible for a poor conductor. And the stability-limited step dt_max = 0.25 \times dx^2 / alpha = 0.25 \times 10^{-4} / 10^{-4} = 0.25 s. Every number is one you can reproduce with a pencil, which is the whole point of a validation report: it turns "the code ran" into "the numbers are right."

Phase 5 — Read the Foundation for What It Foretells

Two of those numbers, side by side, already tell you something sobering about the simulation to come. The plate needs on the order of tau = 10^4 seconds to reach steady state, and each explicit step may advance time by at most dt_max = 0.25 s. The ratio is the number of steps the simulation must take: $\tau / \mathrm{d}t_{\max} \approx 10^4 / 0.25 \approx 4 \times 10^4$ — roughly forty thousand timesteps, each sweeping every one of the $101 \times 101$ grid points. That is a few hundred million grid-point updates for a single, modest run, and it is exactly why the performance (Part VII) and parallel (Part VIII) chapters of this book exist. You have not written the time loop yet — that is Chapter 4 — but your foundation has already sized the computation, and sizing the computation before writing it is what separates a numerical scientist from a hopeful one.

Notice, finally, that refining the grid is not free. Double nx to make the picture twice as sharp and dx halves, so dx^2 quarters, so dt_max quarters — you need four times as many steps, each over four times as many points: sixteen times the work for two times the resolution. That $\mathrm{d}x^2$ scaling is a fact about explicit schemes you will meet rigorously in Chapter 24, but you can read it off your own foundation today, which is the reward for building the numbers carefully and looking at them.

Discussion Questions

1. We defined dp in a module and exported only it. What concrete problem would arise if, instead, each of the book's forty-odd programs declared its own integer, parameter :: dp = selected_real_kind(15, 307) locally? Give a scenario where the copies could silently disagree.
2. Every derived quantity (dx, tau, dt_max) is a parameter, computed at compile time. What are the advantages of that over computing them into ordinary variables at run time? Is there ever a reason you couldn't make such a value a parameter?
3. The report prints alpha and dx^2 in es format but length and dt_max in f format. Justify each choice. When is scientific notation the responsible descriptor, and when does it just make numbers harder to read?

Your Turn: Extensions

• Option A. Add a second material to the report: copper, with alpha ≈ 1.1e-4_dp m²/s. Recompute tau and dt_max for it (you may add a second set of parameters). Which material forces the smaller timestep, and does that match your physical intuition about how fast copper conducts?
• Option B. Make the grid resolution a single knob: change only nx and predict, before running, what dx, dx^2, and dt_max become for nx = 201. Then confirm. State the general rule relating a change in nx to the change in dt_max.
• Option C. Add a compile-time sanity check: a parameter that computes the estimated number of steps as ceiling(tau / dt_max) and prints it with i0. Predict the value, then verify. (Watch the type of ceiling — what does it return, and why is that the right type for a step count?)

Key Takeaways

• The two foundational decisions of a scientific code — precision and physical constants — should be made once, explicitly, in a shared place: dp in a kinds module, constants as parameters with their units.
• State precision as a requirement (selected_real_kind(15, 307)), not a representation (real(8)), so it stays correct and portable across compilers and machines.
• Build derived quantities through visibly-real arithmetic (real(nx - 1, dp)), so no future edit can turn a real division into an integer one behind your back.
• A validation report — print the foundation and check every number by hand — converts "the code compiled" into "the numbers are physically right," and can reveal the shape and cost of the whole computation before a single line of physics is written.