Why this codebase
The Lomax–Pulliam–Zingg two-book CFD curriculum sets the theory and the algorithms in one consistent voice: Book 1 is the analytical spine — operator analysis, modified wavenumbers, stability regions, modified-PDE character — and Book 2 is the algorithmic spine, three structurally different solver families on the same governing equations. Reading the books is one thing; building the codebase that runs every result they describe is the thing that makes the theory operational.
CFDLab is that codebase — mine, written from scratch in Python, still developing.
What it is
An analytical-probe layer, three structurally different quasi-1-D Euler solvers, two 2-D lanes in verification-gated development, and a shared library of the numerical fundamentals everything else stands on — mesh generation and quality, spatial discretization, implicit solves, viscous operators, Riemann and dissipation methods, exact solutions, multigrid transfer — each guarded by tolerance-gated checks. Over a hundred numerical baselines are pinned by per-project MD5 hashes and enforced by the repository’s checks on every push (100+ as of July 2026; the repo is private, so the counts live here rather than behind a link). A 74-entry worked Exercise Guide spans both books.
The reading order on this page mirrors the codebase. Analytical foundations first — the modified wavenumber, the TVD region, modal decay rates that everything else has to obey. Then the three 1-D solver projects, each on the canonical problem it was built to handle, and the cross-project JST-vs-Roe shock-tube comparison that closes the 1-D arc. The two-dimensional lanes come last, because that is where the work now lives.
Analytical foundations
Three figures, all from the analytical-probe layer. The modified wavenumber sets the reading frame: every scheme below the line is a choice in phase and amplitude error per resolved wave. The Sweby diagram is the bridge from pure FD analysis to shock-capturing — limiters live inside that region or they are not TVD. The diffusion modal decay closes the lead trio with the time-domain side of the same machinery: discrete decay rates against the analytical reference, mode by mode.
These are not anchored to any one solver. They are the language the solvers below speak.
Three solver projects
The three Book 2 solvers are structurally different. ARC1D is the Ch.4 lane — implicit finite difference with JST scalar artificial dissipation in delta-form, solved by block-Thomas. FLOMG is Ch.5 — explicit finite volume with five-stage Van der Houwen Runge–Kutta time-marching wrapped in a FAS multigrid driver. upwind1d is Ch.6 — Roe approximate Riemann + MUSCL piecewise-linear reconstruction + TVD slope limiting + SSP-RK2.
Each runs on the canonical problem the project was built to handle. ARC1D does Sod and the transonic nozzle; FLOMG does steady transonic with multigrid; upwind1d does Sod with shock-capturing. The figures here are per-solver — each solver, on its own, doing its own work.
A note on FLOMG: the explicit-FV multigrid lane is purpose-built for steady flows with a single shock, not for the Sod shock-tube transient. Showing FLOMG on transonic rather than on Sod is the honest figure for that lane.
Cross-project shock-tube comparison
With the analytical frame and the three solver lanes on the screen, the JST-vs-Roe Sod overlay reads correctly: the rarefaction → contact → shock structure across the whole domain (full overlay), then the contact zoom that makes Book 2’s central algorithmic argument visible in one image. Measured from each solver’s own output at the sample time, on grids of equal 0.025-wide cells, the scalar-pressure-sensor JST artificial dissipation spreads the contact discontinuity over roughly 1.5 to 1.8 times as many cells as characteristic-decomposition Roe + MUSCL upwinding does — 16 cells against 9 when the contact is counted as the cells lying more than 10 % of the density jump from either plateau, and 20 against 12 at a 5 % tolerance. Neither scheme holds the contact in a single cell, and the ratio moves with the tolerance chosen. That is one 1-D model problem at one grid spacing under one width metric rather than a general ranking of the two schemes, and it is the result the 1-D arc was built to reach.
Into two dimensions
Moving past quasi-1-D changes the problem more than it multiplies it. Geometry arrives: the grid becomes curvilinear, and every derivative picks up metric terms that are themselves computed quantities with their own error. The implicit operator that was a single block-tridiagonal solve in 1-D is too large to invert directly, so it is approximately factored into direction-wise sweeps — a decision with stability consequences of its own. Boundaries become curves with their own discrete treatment, viscous terms bring their own operators and Jacobians, and the mesh itself becomes a numerical artefact that has to be generated, smoothed, and graded for quality before the flow solver ever sees it.
CFDLab’s two 2-D lanes take that step from opposite directions. ARC2D reconstructs the implicit curvilinear method of the NASA Ames lineage — Beam–Warming approximate factorization with the Pulliam–Chaussee diagonal form on a C-grid. These are the general-purpose implicit codes Ames introduced in 1977–78, and the diagonal form earned its keep by cutting the work and storage per step on the vector supercomputers of that era. CFDLab does not reproduce the historical code; it rebuilds and studies the numerical method itself, in a modern Python setting where every milestone is checked before the next is attempted. The lane is built and milestone-gated: the transonic airfoil case below captures the supersonic pocket and shock with the residual driven to a few parts in a thousand of its starting value. It is a qualitative reconstruction — the figure says so on its face — and it is not presented as validated CFD. The reference-data comparison that would earn that word is still ahead.
upwind2d extends the Ch.6 upwind family instead: HLLE and Roe fluxes through cell faces, direction-by-direction MUSCL reconstruction, SSP-RK2 in time — a finite-volume method that is genuinely two-dimensional rather than a rotated 1-D scheme. Its discipline is verification against cases with known answers: reduced to one dimension it reproduces upwind1d’s results, and on the regular oblique-shock reflection below the computed region states land on the derived analytic solution — the dashed lines over the field are the analytic shocks. Where no exact answer exists, results are labelled as comparisons against published figures, never called validated.
What it earns
Building three structurally different solvers against the same governing equations, with the analytical layer right there to consult, is what lets me read a method-property probe before running a case and know what the case will say. It is the same discipline that makes a workflow elsewhere — an OpenFOAM run on the cluster, a high-order SBP scheme on a model PDE — give up its pathologies before a dashboard hides them. The 2-D lanes are that discipline meeting more machinery: curvilinear metrics, factored implicit operators, real meshes — and the same rule holding throughout, that nothing gets called more than what its checks have earned.
The tagged plateau is v0.6-two-book-exercise-guide (April 2026). The work is now well past that tag — the 2-D lanes and the shared fundamentals library are where it is moving — with no newer formal release yet.
Analytical foundations
Three solver projects
Into two dimensions
