26 March 2026 18 min read

Measuring a Fluid from Many Truncations at Once

A consistency check across Galerkin truncations, the solver bug it found, the cross-check on the fix, and what the corrected solver shows an inviscid truncation doing.

The Navier–Stokes equations describe viscous fluids. On a computer they are solved by truncation: keep the Fourier modes with |k| ≤ N, drop the rest, and integrate the finite system that remains. Every such run is a different dynamical system, and none of them is the PDE. The question that matters for any computation is how much of what you see is the fluid and how much is the truncation.

The scaffold array is a diagnostic for that question. This post is about what it measured: a defect in my own solver, the correction, a cross-check against an independent solver architecture, and what the corrected solver shows an inviscid truncation doing.

The Method: One Datum, Many Truncations

Take a single initial velocity field u₀. Project it onto the Galerkin truncation at N = 3, at N = 4, and so on up to N = 8 or more, with the new modes at zero. Run each truncation with the same viscosity and the same time step. Record the same observables from each: total energy E(t), enstrophy Ω(t), the per-shell balance between nonlinear transfer into shell k and viscous drain out of it, and the distribution of energy across shells.

Now compare across N. If an observable settles down as N grows, the truncations agree about it. If it does not, they disagree, and the pattern of which observables disagree first is informative. Totals and per-shell balances integrate different things. An error that lives at the truncation boundary shows up in the shell-resolved quantities before it shows up in the totals; an error that lives in the operator itself shows up everywhere at once.

That is all the scaffold array is: the same measurements taken at several truncation levels and read side by side. It does not tell you what the PDE does. It tells you whether your truncations are measuring the same thing.

What the Comparison Showed

Run on my original triad kernel, the array showed a pattern I could not attribute to physics. Enstrophy grew with N without settling. The per-shell balances were well-behaved at every shell individually. The disagreement was in the aggregate, not at any particular scale, and that pointed at the operator.

The decisive experiment was the plain one. Set ν = 0. The nonlinear term of Navier–Stokes conserves energy exactly, so an inviscid truncation must hold E(t) constant up to time-stepping error, and time-stepping error shrinks as the step shrinks. I ran the inviscid truncation at step sizes spanning eight orders of magnitude:

N	Modes	Drift at dt = 10⁻⁴	at dt = 10⁻⁸	at dt = 10⁻¹²
3	122	−0.37%	−0.37%	−0.37%
5	514	+1.09%	+1.09%	+1.09%
7	1,418	+14.87%	+14.87%	+14.87%
8	2,108	+12.25%	+12.25%	+12.25%

The drift did not depend on the time step. Whatever was injecting energy, up to 15% at N = 7, was in the right-hand side, not in the integrator.

The Defect

In Fourier space the Navier–Stokes nonlinearity is

dû_k/dt = −i Σ_p+q=k (û_p · q) P_k û_q (plus the viscous term).

The kernel stored the coefficients as real triples and computed Σ (û_p · q) P_k û_q with no −i. That is a different trilinear form. It has the same triads, the same Leray projection and the same coupling graph, and it does not conserve energy: the antisymmetry that makes Σ Re(conj(û_k) · NL_k) vanish depends on the factor −i acting on complex coefficients. Without it the sum is not zero, and the amount by which it fails grows with the number of interacting triads, which is why the drift grew with N.

The correction, made on 24 March 2026, is to store complex coefficients and apply the −i. With it the kernel conserves energy at ν = 0 at every N, energy decreases at ν > 0 at every N, and the Taylor–Green vortex decays as exp(−6νt). An independent NumPy implementation and a scipy RK45 integration of the same truncation agree with the C kernel on E(t).

The Cross-Architecture Check

A triad kernel sums over the coupling graph explicitly. A pseudo-spectral solver computes the same nonlinear term by transforming to physical space, multiplying, and transforming back with an FFT, with dealiasing. The two share no code and no data structures. Run from the same datum at N = 4 over T = 0.01, the FFT solver and the corrected triad kernel agree to 2.8 × 10⁻¹⁷. That is the check that the corrected kernel computes the Navier–Stokes truncation and not some third system.

What the Corrected Solver Shows

With the corrected kernel, run the inviscid truncation for a long time. Enstrophy grows, and the growth is now a property of the truncated system rather than of a defective kernel. It has a well-known explanation. A finite set of modes that conserves energy and interacts through a cascade goes to equipartition: the Lee–Kraichnan absolute equilibrium, in which every retained mode holds the same energy. In that state Ω/E equals the mean of |k|² over the retained lattice modes, ⟨k²⟩.

The runs reach it. Ω/E lands on ⟨k²⟩ to within 1% at every N from 3 to 8, and the value does not change when the time step is halved. Enstrophy growth in a truncated Euler run is thermalisation, and its end state is a number set by N, not by the flow. That is what the corrected solver measures, and it is a statement about truncations.

A Note on Validation

The original kernel passed the usual tests. It converged under time-step refinement. It preserved incompressibility. It matched a known solution at low N. It produced a plausible-looking cascade. The one test it failed is the one that is not usually run: inviscid energy drift across several time steps. If the drift is independent of dt, the error is in the operator. The test takes minutes on any spectral solver, and I now run it first.

The general principle: a conservation law the equations satisfy exactly is a test the solver must pass exactly, and the way to tell an operator error from an integrator error is to vary the step and see whether the violation moves.

Try It

If you have a C compiler and Python, the checks take about twenty minutes:

Clone the repo: github.com/senuamedia/lab-code
Run the energy audit: build the v3 triad kernel, run at ν = 0. Energy is conserved at every N.
Run the Python check: an independent NumPy implementation. Same result, no shared code.
Reproduce the defect: build the same audit against the v2 kernel, which is kept in the repository, and see the dt-independent drift.

Everything compiles with gcc or clang, and Python needs only numpy and scipy. The README has the commands.

If you work with spectral fluid solvers, for weather, aerodynamics, ocean modelling or research, run the inviscid drift test on yours. If you work with Lean 4, Coq or Isabelle/HOL, the energy identity of the corrected kernel is a clean algebraic property of an 800-line C file, and an independent formalisation would be a useful contribution.