Case Study 1: First Day on a Hundred-Thousand-Line Code

"I don't need to read the whole code. I need to find the one routine that matters, and I need to find it before lunch."

Executive Summary

You have joined a research group and been handed the URL of thermex, a mature open-source thermal-transport code — a hundred thousand lines of Fortran you have never seen, written by twenty people over fifteen years. Your first task is small and specific: a colleague reports that the diffusion coefficient looks wrong in one material region, and asks you to find where the diffusivity is set and how it reaches the solver. You will not — cannot — read the code to answer this. Instead you will use the architecture of §36.1–§36.5 as a map, and grep and tags as your feet, to walk straight to the answer having read a fraction of a percent of the code. This case study is that walk, step by step, and its lesson is the one every maintainer of a large scientific code has internalized: you navigate by structure, not by reading.

Skills applied: reading a source tree and module hierarchy (§36.1); classifying modules by role (§36.2); reading the build file (§36.3); finding the entry point and building a module map (§36.4); locating where state lives and following the data flow (§36.5). It exercises the reading discipline that Chapter 17 applies to old code, on a large code.

Background

thermex is a composite stand-in for the real thing, and everything you do here transfers directly to an actual open-source Fortran code — WRF, CESM, or Quantum ESPRESSO, all large, modular, open, and laid out in recognizable layers (Chapter 1). Clone any one of them alongside this case study and run the same commands; the specifics of the tree will differ, but the moves will not. We keep thermex's internals generic on purpose — the transferable skill is the method, not any one code's module names.

Your task, precisely: find where the diffusivity alpha is set, and trace how it flows into the solver.

Phase 1 — Orient: README and Build File (Five Minutes)

You do not open src/ yet. You read README.md and the build file first, because together they tell you what the code is and what it is made of. The README says thermex solves a heat-transport equation on a structured grid and is built with fpm. The fpm.toml confirms the shape:

name = "thermex"
version = "3.4.1"
[build]
auto-executables = true

That one word — fpm — already tells you the layout (§36.3's convention): the library is in src/, the program in app/, tests in test/. You now know where to look without having looked.

Phase 2 — Find the Entry Point

Every program begins in exactly one place. Find it:

$ grep -rin "^\s*program " app/ src/
app/main.f90:6:program thermex

You open app/main.f90 and — because a healthy driver is thin (§36.2) — the whole shape of a run fits on one screen:

program thermex
  use config_io,   only: read_config
  use material,    only: set_material
  use solver,      only: run_solution
  use output_io,   only: write_results
  ! ... read the config, set up materials, run, write ...
end program thermex

In sixty seconds you have the skeleton of a hundred-thousand-line code: it reads a config, sets up materials, runs a solution, writes results. The word material leaps out — that is where diffusivity almost certainly lives.

Phase 3 — Build the Module Map, Follow the Name

You do not read material. You query for your name, alpha. First, where is it declared?

$ grep -rin "alpha" src/material.f90 | head
src/material.f90:8:  real(dp) :: alpha = 0.0_dp     ! diffusivity (module variable)
src/material.f90:22:    alpha = table(region)%diffusivity

There it is, and the comment tells you something crucial about where the state lives (§36.5): alpha is a module variable, not an argument. It is global state. That single fact reshapes your search — because a module variable is a hidden input to every routine that reads it, you must now find everyone who assigns it. Grep for the assignment across the whole tree, case-insensitively:

$ grep -rin "alpha *=" src/
src/material.f90:22:    alpha = table(region)%diffusivity   ! set here, from the region's material
src/solver.f90:41:    a_local = alpha                        ! read here, copied into the step

Two hits, and they are the whole story: alpha is written in material.f90 (from a per-region table) and read in solver.f90 (copied into the stepping routine). You have found the data flow in two greps.

Phase 4 — Reproduce the Mechanism in Miniature

To be sure you understand the mechanism — a module variable set by one routine and read by another — you reproduce it as a small, compilable model. This is a habit worth keeping: when an unfamiliar mechanism matters, rebuild the shape of it in ten lines you fully control.

module material                       ! the "where state lives" mechanism, distilled
  implicit none
  private
  integer, parameter :: dp = selected_real_kind(15, 307)
  real(dp) :: alpha = 0.25_dp         ! GLOBAL diffusivity; a default until set_material runs
  public :: set_material, diffusivity
contains
  subroutine set_material(a)          ! the ONE writer (thermex reads it from a table)
    real(dp), intent(in) :: a
    alpha = a
  end subroutine set_material

  real(dp) function diffusivity()     ! a reader: its result depends on hidden state
    diffusivity = alpha
  end function diffusivity
end module material

program trace_alpha
  use material, only: set_material, diffusivity
  implicit none
  print '(a, f6.2)', 'alpha (default)  : ', diffusivity()   ! before any set_material
  call set_material(0.50_dp)                                 ! a distant line changes it
  print '(a, f6.2)', 'alpha (after set): ', diffusivity()   ! same call, new answer
end program trace_alpha

$ gfortran -std=f2018 -Wall -O2 trace_alpha.f90 -o trace && ./trace
alpha (default)  :   0.25
alpha (after set):   0.50

The output is exact and it is the bug's mechanism: diffusivity() returns $0.25$ then $0.50$ from identical calls, because set_material reached into the shared state between them. In the real thermex, that is precisely how a wrong region-to-material mapping in the table silently poisons the solver: the solver reads whatever the last set_material left in alpha, and nothing at the read site reveals it. You have localized the fault to the assignment in material.f90:22 — the table lookup — without ever reading the solver's numerics.

Phase 5 — Report, and the General Lesson

Your report to the colleague is now precise: alpha is global module state, written once per region from a material table (material.f90:22) and read by the solver (solver.f90:41); if a region maps to the wrong table row, every cell in that region gets the wrong diffusivity, invisibly. The fix is a table correction (or, better long-term, threading alpha explicitly per region so the dependency stops being hidden — §36.5). You found it in five commands.

Tally what you actually read: the README, a paragraph of fpm.toml, the ~30-line driver, and about ten lines across two modules. Call it 50 lines of a 100,000-line code — 0.05%. That is the thesis of §36.4 made concrete: you were productive in a code you have overwhelmingly not read, because you navigated its architecture and followed one name instead of reading front to back.

Discussion Questions

1. In Phase 3, the comment ! diffusivity (module variable) was decisive — it told you alpha was global state. If the comment had been absent, what single grep would have told you the same thing (that alpha is a module variable, not a local or an argument)?
2. You found alpha was assigned in exactly two places. Suppose the second grep had returned twelve assignments across eight files. How does that change your confidence, and what does a large number of writers to one global variable tell you about the code's design (§36.5)?
3. The driver in Phase 2 was thin and readable in sixty seconds. What would you conclude about the code's maintainability if app/main.f90 had been 3,000 lines that did the setup, the stepping, and the output inline, with no solver module at all?

Your Turn: Extensions

• Option A. Clone a real open-source Fortran code (WRF, CESM, or Quantum ESPRESSO). Without reading it, run the Phase 2–3 commands: find the program, list the modules with grep -rin "^\s*module ", and pick one physical quantity (a coefficient, a boundary value) and grep for where it is assigned. Time yourself. Report how much you learned and how little you read.
• Option B. Extend the trace_alpha model: add a second reader function in a different module that also depends on alpha, and demonstrate the "hidden output" hazard — call set_material between the two readers and show that changing state for one silently changed the other. This is "spooky action at a distance" in fifteen lines.
• Option C. Take the material model and remove the global state: rewrite diffusivity to take alpha as an argument, so its result depends only on its inputs. Write the one-line compile-time consequence a caller now faces (it must supply alpha), and argue why that inconvenience is exactly the point.

Key Takeaways

• The first moves in any large code are always the same: README and build file, then the program entry point, then the module map — orient before you read.
• Follow a name, not a file. grep for where a value is declared, then for where it is assigned and read; two or three searches usually trace the whole data flow.
• A comment or a grep revealing that a value is a module variable is a turning point: global state has hidden inputs, so you must find every writer to know what a reader sees.
• You can be productive in a code you have read almost none of. Navigating 0.05% of thermex answered the question — because architecture, not reading, is how large scientific software is understood.