# waterfall: Validation Report

**Date:** 2026-07-27
**Hardware:** AMD Ryzen 9 5900X (12C/24T), 64 GB RAM, Radeon RX 6800
**Scope:** Phase 1: establish whether the solver's mathematics is trustworthy
before anything is built on top of it.

---

## 1. Verdict

The solver that existed at the start of this work contained **no fluid
mechanics** (§2). It has been replaced with an incompressible Navier–Stokes
solver.

| | Before | After |
|---|---|---|
| Momentum equation | none | ✅ solved |
| Mass conservation (∇·u = 0) | none | ✅ enforced to 1e-16 |
| Reynolds-number dependence | **none at any Re** | ✅ correct trend |
| Blasius boundary layer | cannot run | ✅ L2 = 0.7–1.3 % |
| Kármán vortex shedding | structurally impossible | ✅ St within 2.6 % |
| Wake structure / separation | no wake region exists | ✅ L_r/D converges into range |
| Bluff-body Cd | geometry heuristic | ⚠️ +3 to +12 %, **does not converge** |
| **Absolute vehicle Cd** | meaningless | ❌ **2.7× high (§5.6)** |

**Verified and usable:** the flow solution. Velocity fields (Blasius L2
0.7–1.3 %), wake geometry (recirculation length converges monotonically into
the reference band), shedding dynamics (Strouhal within 2.6 %), exact mass
conservation (∇·u ~ 1e-16 in every case run).

**Not verified:** force magnitude. The grid-convergence study (§5.8) shows Cd
moving *away* from a limit under refinement (1.610 → 1.622 → 1.779 across
16/24/32 cells/D, observed order −6.3) while the flow field in the same runs
converges properly. Four candidate causes were tested and eliminated by
measurement (§6.1); the evidence points at the force estimator itself, and the
recommended fix is a control-volume momentum balance. **Treat drag numbers as
±10 % and do not compare across resolutions.**

**Do not use for absolute vehicle Cd.** The Model S returns 0.557 against a
published 0.208. That run is stable and divergence-free; it is correctly solving
a problem at an *effective* Reynolds number near 850, because subgrid viscosity
on an affordable uniform grid reaches 1.2 × 10⁴ times molecular (§5.6). This is
a property of the method at achievable resolution, not a bug.

The honest summary: **the physics engine is verified; the force estimator is
not, and the resolution a 30-minute budget permits is not sufficient for
absolute automotive drag.** The rewrite replaced a solver with no fluid
mechanics with one whose remaining errors are measured, bounded, and traced to
a specific component.

---

## 2. What the previous solver actually computed

This matters because old `summary.txt` files exist and their numbers are not
what they appear to be.

The entire physics kernel was ([legacy/ray_solver.py:401](legacy/ray_solver.py#L401)):

```python
angle = angle_between(direction, nearest_face_normal)
loss  = angle**2
pressure[i,j] = loss              # dimensionless
energy       *= exp(-loss)
velocity[i,j] = energy            # a monotonically decaying scalar
```

Findings from the audit:

1. **No governing equations.** No advection scheme, no pressure solver, no
   continuity, no momentum balance, no turbulence model. Every grid cell was
   independent: there was no neighbour stencil anywhere, hence no pressure
   gradient, no viscous coupling, and no mechanism by which a wake could form.

2. **No Reynolds-number dependence.** `viscosity` appeared in `config.json`
   and was read by nothing. Output was identical at Re = 10 and Re = 10⁷.

3. **The "rays" never moved.** Direction was updated each slice but position
   was never integrated, so each sample point stayed at a fixed (x, y) for the
   whole traverse. The "4D ray carry-over" was a decaying scalar per fixed
   column, not a streamline.

4. **Sampled the bounding box, not the surface.** `sample_surface_normal` ran a
   nearest-point query for every polar grid point, including points deep inside
   the body and far out in free stream, and discarded the returned distance.
   Normals were assigned to empty air.

5. **Drag scaled with slice count.** The force integral summed
   `pressure × q × cell_area` over the full footprint *for every slice*, so
   total integrated area was footprint × n_slices. Halving `slice_depth`
   doubled the drag.

6. **Negative drag: root cause identified.** The per-cell drag map was clamped
   ≥ 0, but the reported force used `−dot(Σ pressure·n̂, flow)` over unclamped
   normals. Two different formulas; the sum could face upstream. This was a
   known issue in the legacy code.

7. **`C_d ≠ F_d /(q·A)` in its own output.** Two dead assignments
   (`legacy/ray_solver.py:548` and `:551`) meant the reported force was scaled
   by one factor and the reported coefficient by a different one.

8. **Reference area was the bounding box**, overstating a car's frontal area by
   15–25 %, applied directly to the denominator of C_d.

**Conclusion:** not "inaccurate physics" but *no* physics. The benchmarks
requested in Phase 1 could not be run against it even in principle: Blasius
requires viscosity and a wall-normal velocity profile (neither existed), and
Kármán shedding requires time integration with neighbour coupling (the scheme
was single-pass with decoupled cells).

---

## 3. What replaced it

```
∂u/∂t + ∇·(u u) = −∇p/ρ + ν∇²u + f_ibm
∇·u = 0
```
Incompressible, isothermal, Newtonian, constant density.

| Component | Method | Rationale |
|---|---|---|
| Grid | Uniform staggered MAC | Enables a direct pressure solve |
| Pressure | **DCT-II Poisson, all-Neumann** | Direct, non-iterative; 57 ms @ 2.4M cells |
| Time | Adams–Bashforth 2, one projection/step | 2nd order at one Poisson solve, not three |
| Convection | 2nd-order central, divergence form | Non-dissipative (see below) |
| Diffusion | Explicit 7-point Laplacian | Convective CFL binds first, so it is free |
| Geometry | Immersed boundary, fractional η ∈ [0,1] | No meshing; STL swap costs seconds |
| Forces | IBM momentum integral (Newton's 3rd law) | Avoids integrating pressure over a staircase surface |

**Why central differencing, not upwind.** Upwinding introduces numerical
dissipation that is mathematically indistinguishable from physical viscosity.
It would corrupt exactly the quantity the Blasius case measures, and it damps
the vortex shedding the cylinder case depends on, producing a plausible,
steady, *wrong* wake. Central differencing is non-dissipative; when
under-resolved it oscillates visibly rather than failing quietly.

**Why the FFT pressure solve is the load-bearing decision.** It is what makes
a 3D vehicle case fit in 30 minutes. Its cost is that the grid must be uniform,
so there is no local refinement near the body. That single constraint accounts
for most of the residual error in §5.

---

## 4. Exact verification (`benchmarks/smoke.py`)

These isolate individual operators. All must pass before any physics result is
meaningful.

| Test | Result | Status |
|---|---|---|
| Poisson vs manufactured solution | max err 8.9e-16 | ✅ machine precision |
| Freestream preservation (50 steps) | max\|u−1\| = **0.0** | ✅ exact |
| ∇·u after projection | 2.2e-16 | ✅ round-off |
| 2D collapse (nz=1 ⇒ w ≡ 0) | max\|w\| = 0.0 | ✅ exact |
| IBM drag sign | F_x = +0.62 N | ✅ positive |

Freestream preservation being *exactly* zero is the strongest single check
here: it means the convective operator, the ghost cells and the projection are
mutually consistent. A solver that fails it can still produce plausible drag
numbers, which is why it is tested first.

---

## 5. Benchmark results

### 5.1 Blasius flat-plate boundary layer

`benchmarks/blasius.py`: 600×192 cells, ν = 2e-4, U = 1, leading edge inside
the domain with a slip section upstream. Reference: similarity solution of
f''' + ½ f f'' = 0, solved by shooting (recovered f''(0) = 0.332057 against the
exact 0.332057).

| x (from LE) | Re_x | pts in BL | L2(u/U_e) | max err | θ err | δ* err | H | cf err |
|---|---|---|---|---|---|---|---|---|
| 0.199 | 1036 | 20 | **0.0073** | 0.0122 | −2.13 % | −9.18 % | 2.41 | +3.75 % |
| 0.401 | 2101 | 28 | **0.0096** | 0.0159 | −2.24 % | −7.84 % | 2.44 | +5.42 % |
| 0.699 | 3687 | 37 | **0.0133** | 0.0218 | −2.86 % | −7.71 % | 2.46 | +7.73 % |

Shape factor H = δ*/θ against the exact 2.59. Momentum thickness (which *is*
the drag integral) is within 3 % at every station. ∇·u held at 2.3e-16
throughout 12,367 steps.

Quantities are normalised by the **local edge velocity**, not U∞. This is
standard boundary-layer practice and it is necessary here: the finite domain
height accelerates the freestream by 4.1–5.5 % (reported as `blockage` in the
JSON), and normalising by U∞ instead makes θ integrate a spurious negative
contribution large enough to flip its sign.

### 5.2 Circular cylinder: Kármán vortex street

`benchmarks/cylinder.py`: 720×480, 24 cells/D, 5 % blockage, far-field sides,
convective outflow.

| Re | Cd | ref Cd | err | St | ref St | err | Cl_rms | ref | L_r/D | ref |
|---|---|---|---|---|---|---|---|---|---|---|
| 40 | 1.622 | 1.48–1.56 | +6.7 % | steady ✓ | n/a | n/a | 0.000 ✓ | ~0 | **2.25** | **2.13–2.35 ✓** |
| 100 | 1.396 | 1.32–1.38 | +3.4 % | **0.1698** | 0.164–0.167 | **+2.6 %** | 0.231 | 0.30–0.34 | 1.16 | n/a |
| 200 | 1.376 | 1.30–1.34 | +4.2 % | **0.2061** | 0.192–0.197 | +5.9 % | 0.458 | 0.60–0.70 | n/a | n/a |

All trends are correct: Cd falls with Re, St rises with Re, Cl_rms rises with
Re, and the Re=40 wake is steady while Re=100 and 200 shed periodically.
∇·u ≤ 5.9e-16 in every case; IBM slip velocity 2.1–2.3 % of freestream.

**Reading the errors.** Cd and St are biased high *together*, which is the
signature of blockage: confinement raises the effective velocity past the
body, lifting both. Recirculation length at Re=40 lands inside the reference
band, and it is a purely geometric measure of the wake that does not depend on
the force calculation, so the near-wake structure is right even where the
force magnitude is a few percent high.

Cl_rms is low by ~28 % at both shedding Reynolds numbers. The consistency of
that bias across Re points at the immersed boundary smearing the surface
pressure distribution that drives the lift fluctuation, rather than at the
shedding dynamics, which the Strouhal number confirms are correct.

Re=200 in 2D is expected to overpredict Cd and Cl_rms, because the physical
wake becomes three-dimensional above Re ≈ 190.

### 5.3 Grid convergence (Re=40)

| cells/D | Cd | L_r/D |
|---|---|---|
| 16 | 1.610 | 2.23 |
| 24 | 1.622 | 2.25 |

Cd changes by 0.7 % across a 50 % refinement while remaining ~6 % above the
reference band. A resolution-driven error would shrink with refinement; this
one does not. That isolates the residual to the domain, not the discretisation,
as tested directly in §5.5.

### 5.4 Sphere, Re=100: 3D path

| Quantity | waterfall | Reference |
|---|---|---|
| Cd | 1.172 | 1.08–1.10 (+8 %) |
| Cl | −0.0000 | 0 (axisymmetric) ✓ |
| ∇·u | 2.5e-16 | n/a |

At **6 cells across the sphere**. Cl vanishing to 5 decimal places confirms the
3D kernels preserve axisymmetry, a check that catches asymmetric indexing
errors in the v/w momentum routines that a drag number alone would hide.

### 5.5 Blockage isolation (Re=40 cylinder)

| Domain height | Blockage | Cd | L_r/D |
|---|---|---|---|
| 20 D | 5.0 % | 1.6221 | 2.250 |
| 40 D | 2.5 % | 1.6165 | 2.259 |

Halving blockage changes Cd by 0.3 %. See §6.1: this falsifies blockage as
the cause of the Cd bias.

### 5.6 Test geometries

#### 2023 Tesla Model S: the important negative result

| Quantity | waterfall | Published |
|---|---|---|
| **Cd** | **0.557** ± 0.011 | **0.208** |
| Cl | +0.090 ± 0.023 | n/a |
| Frontal area (rasterised) | 2.724 m² | ~2.4 m² |
| Re_L | 9.9 × 10⁶ | n/a |
| Drag force | 836 N @ 30 m/s | n/a |
| Grid / runtime | 776 k cells, 36 cells/car | **20.3 min** |
| ∇·u | 3.1e-16 | n/a |
| **max ν_t/ν** | **11,716** | n/a |

**Cd is 2.7× the published value. Do not use this number.**

The diagnostic that explains it is the last row. The Vreman model is producing
an eddy viscosity nearly twelve thousand times molecular; `sgs.py` documents
that above ~100 the model is doing more work than the resolved scales. At
ν_t/ν ≈ 1.2e4 the effective Reynolds number is
30 × 4.97 / (1.17e4 × 1.5e-5) ≈ **850**, not 10⁷.

So the run is not a Re=10⁷ simulation with some modelling error. It is
converged, divergence-free, stable, and solving a Reynolds number four orders
of magnitude below the one requested. A bluff body at Re ≈ 10³ genuinely does
have Cd near 0.5–0.6, which is why the answer looks plausible and is not.

This is not a bug; it is the arithmetic of LES on a coarse grid. Subgrid
viscosity scales as (CΔ)²|S|, and at Δ = 138 mm that is simply a large number.
Two things follow, and they matter more than the drag value:

1. **A streamlined body is the worst case for this solver.** The Model S gets
   Cd 0.208 from smooth pressure recovery over a long rear taper. Represent
   that taper as a 138 mm staircase and drown it in eddy viscosity, and the
   recovery is destroyed: the flow separates early and the body behaves like
   a bluff one. The benchmark results in §5.2–5.4 are all *bluff* bodies, where
   separation is pinned by geometry; they do not transfer to this case.
2. **Absolute vehicle Cd is outside this solver's validated envelope**, at any
   resolution reachable in 30 minutes. What remains valid is comparison: two
   variants at identical settings share the bias.

Honest accounting of the runtime target: this case does meet it (20.3 min
against an estimate of 27.3), but meeting a runtime budget on a case whose
answer is not trustworthy is not a success. The budget and the accuracy are in
direct conflict here, and the budget won.

The first attempt at this case, with **no** subgrid model, went unstable at
t = 2.6 s (see §6.7). The pressure/wake *fields* from the successful run are
physically sensible (stagnation Cp → +1 at the nose, suction over hood and
roof, coherent wake structure in the Q-criterion iso-surface), which is why the
visualisation is still useful even where the coefficient is not.

#### Front bumper example

Pre-flight completed. The case file and geometry for this part are not shipped;
its run summary and stills are in `results/front_bumper_example/`. A deliberate
caveat applies: **a bumper in isolation has
no meaningful absolute Cd.** There is no reference-area convention for a
component, and most of a real bumper's aerodynamic effect is how it feeds the
underbody and wheel wakes, none of which exists in a domain containing only
the bumper. Use it to compare variants of the part under identical conditions.

One pre-flight catch worth recording, because it would have been invisible in
the result: the silhouette method reported a frontal area of 0.637 m² against
a true rasterised 0.206 m², because it fills through-gaps. Left uncorrected
that would have divided Cd by 3.

### 5.7 Square cylinder: sharp-edged, grid-aligned

`benchmarks/bluff_body.py`: 720×480, 24 cells/D.

| Re | Cd | ref Cd | err | St | ref St | err | Cl_rms |
|---|---|---|---|---|---|---|---|
| 100 | 1.521 | 1.44–1.51 | +3.1 % | 0.1547 | 0.144–0.150 | +5.2 % | 0.555 |
| 200 | 1.598 | 1.37–1.48 | +12.1 % | **0.1535** | 0.148–0.156 | **in range ✓** | 0.555 |

This case exists to separate two things the circular cylinder confounds, and
it does so decisively.

A square's faces align exactly with the Cartesian grid, so it is the most
favourable geometry a Cartesian immersed boundary will ever see: there is no
staircase approximation of curvature at all. It nonetheless shows **+3.1 % at
Re=100, statistically the same as the circle's +3.4 %**.

That is strong evidence against staircasing being the cause of the systematic
bias of §6.1. If representing a curved surface as a stair-step were the
mechanism, the square would be markedly better and it is not. What both bodies
share is the diffuse η layer (roughly one cell of smeared surface), which
points to the effective-body-size error as the mechanism, and predicts it
should shrink as h¹. That is what the convergence study in §5.8 measures.

Re=200 at +12.1 % is the expected 2D penalty: the physical wake behind a square
cylinder is three-dimensional by that Reynolds number, so a 2D computation
overpredicts. Note that St stays inside the reference band even there: the
shedding *dynamics* remain right while the force magnitude degrades, which is
the same pattern as the circular cylinder.

### 5.8 Grid convergence, Re=40 cylinder: the force does not converge

Fixed domain (30D × 20D), fixed Re, current force treatment, three resolutions.

| cells/D | Cd | L_r/D | IBM slip | Cd std | ∇·u | steps |
|---|---|---|---|---|---|---|
| 16 | 1.6100 | 2.224 | 0.031 | 0.0001 | 2.5e-16 | 7,813 |
| 24 | 1.6221 | 2.253 | 0.021 | 0.0001 | 5.9e-16 | 12,606 |
| 32 | **1.7785** | 2.260 | 0.013 | 0.0001 | 2.5e-16 | 20,562 |
| reference | 1.48–1.56 | 2.13–2.35 | n/a | n/a | n/a | n/a |

Observed order of accuracy: **−6.3**. A negative order means the quantity is
moving *away* from a limit under refinement. The Richardson extrapolation the
script prints (1.5918) is therefore meaningless and must not be quoted.

Read the columns against each other, because they disagree in an informative
way:

- **The flow field converges.** Recirculation length rises monotonically
  2.224 → 2.253 → 2.260, settling inside the reference band. Immersed-boundary
  slip falls 0.031 → 0.021 → 0.013, roughly as h¹, exactly as it should.
  Divergence is at round-off throughout.
- **The force does not.** Cd moves +0.7 % from 16→24 and then +9.6 % from
  24→32. `Cd std` is 1e-4 at every resolution, so all three are genuinely
  steady converged states, not unsettled runs or noise. The 32-cell value
  reproduces an anomaly observed in an independent earlier run (1.778), so it
  is not a one-off.

This falsifies the last surviving hypothesis of §6.1. A first-order
effective-body-size error would make Cd *decrease* toward the reference under
refinement. It increases, and accelerates.

**What this means for the deliverable.** The claim "validated drag prediction"
is not supported and is not made anywhere in this report. What *is* supported:

| Claim | Status |
|---|---|
| Momentum and continuity solved correctly | ✅ exact checks, §4 |
| Velocity fields accurate | ✅ Blasius L2 0.7–1.3 % |
| Wake structure and separation accurate | ✅ L_r/D converges into range |
| Shedding dynamics accurate | ✅ St within 2.6 % |
| **Force magnitude** | ❌ **does not converge; treat as ±10 %** |

**Recommended fix, highest priority.** Replace the direct-forcing momentum
integral with a control-volume momentum balance: integrate momentum flux and
pressure over a box enclosing the body, well away from the immersed surface.
That formulation never touches the forced cells, so it is insensitive to the
η layer and to the 1/dt amplification in the present estimator. The force is
currently a small residual of a near-cancellation divided by a timestep that
shrinks under refinement, which is a plausible mechanism for error growth and
is worth testing directly.

Until that is done, use waterfall for flow structure and for relative
comparison at **fixed resolution**. Changing resolution between two cases
being compared invalidates the comparison, since the bias is resolution
dependent.

---

## 6. Where the solver diverges from truth, and why

Ranked by how much they will affect a vehicle drag number.

### 6.1 The systematic +3 to +8 % Cd bias: cause NOT yet established

Every bluff-body case overpredicts Cd by a similar margin: cylinder +6.7 %
(Re=40), +3.4 % (Re=100), +4.2 % (Re=200), sphere +8 % (Re=100). The
consistency says it is one mechanism, not four coincidences.

Two candidate causes were tested and **both are ruled out by measurement**:

**Blockage: ruled out.** Halving blockage from 5 % to 2.5 % (domain height
20D → 40D, resolution held fixed) moved Cd from 1.6221 to 1.6165, a change of
0.3 %, and left recirculation length unchanged at 2.26. Blockage does
demonstrably accelerate the freestream (it is measured at 4.1–5.5 % in the
Blasius case), but it is not what is inflating cylinder drag. An earlier draft
of this report asserted it was the dominant term; that was wrong.

**Immersed-boundary inner-mass term: ruled out.** The direct-forcing force is
F = −ρ∫f dV + ρ d/dt∫_solid u dV, and the second term is only negligible at
statistical steady state. It was measured: momentum inside the solid drifts at
5.6e-10 per step, i.e. zero. The first term alone is the correct force.

Worth recording because it looks alarming and is not: **67.5 % of the reported
force comes from cells deep inside the body**, not from the surface band.
That is correct behaviour, not a bug. Inside the solid region the projection
establishes the front-to-back pressure difference across the body, and the
forcing that opposes it *is* pressure drag. Restricting the force integral to
the surface band gives Cd = 0.539 against a reference of ~1.5; the interior
contribution is the physics, not an artefact.

**Staircasing of curved surfaces: ruled out.** The square cylinder (§5.7) has
faces exactly aligned with the grid, so it suffers no staircase approximation
whatsoever. It shows +3.1 % at Re=100 against the circle's +3.4 %, the same
bias. If stair-stepping a curved surface were the mechanism, the square would
be markedly better. It is not.

**Effective body size / first-order η layer: ruled out.** This was the last
surviving hypothesis and the convergence study of §5.8 kills it. A half-cell
oversized body produces an error that *shrinks* under refinement; the measured
Cd instead rises 1.610 → 1.622 → 1.7785 across 16/24/32 cells/D, an observed
order of −6.3. Meanwhile recirculation length converges monotonically into the
reference band and IBM slip falls as h¹, so the flow field is converging while
the force is not.

**Conclusion: the cause is in the force estimator, not the flow solution.**
All four geometric and boundary explanations are eliminated by measurement, and
the one quantity that fails to converge is the one computed by a different
route from everything else. The present estimator sums the immersed-boundary
forcing and divides by dt: a small residual of a near-cancellation, divided by
a number that shrinks under refinement. That is a credible mechanism for error
that grows with resolution, and it is consistent with every observation here,
but it has not yet been confirmed directly.

**Recommended fix:** control-volume momentum balance over a box enclosing the
body, away from the forced cells. See §5.8. This is the highest-value next
piece of work on the solver.

The earlier two-point resolution data (16 → 1.610, 24 → 1.622) cannot settle
this: those runs predate the force filter of §6.5, so their means are
contaminated by the Nyquist mode and are not comparable to each other.

**Practical consequence, independent of cause:** absolute Cd carries a
systematic few-percent-high bias. Differences between two geometries run at
identical settings do not; the bias is common-mode and cancels. Use it for
ranking, not for quoting.

### 6.2 Immersed boundary is first-order

The solid fraction η smears the surface over roughly one cell, so the
effective body is slightly larger than the CAD body and the error decays only
as h¹.

This is why the reference area is taken from the **rasterised** mask rather
than the CAD silhouette: numerator and denominator then describe the same
object, so refinement converges instead of chasing two different bodies. The
effect is large enough to see: the Model S rasterises to 2.72 m² at 36
cells/length against 2.43 m² at 40, a 12 % change in the denominator of Cd
from a 10 % change in resolution.

**Mitigation:** ≥ 20 cells across the body minimum, ≥ 40 preferred. Below 20
the pre-flight warns.

### 6.3 No wall model: friction drag is not physical

The first cell off the body sits far outside the viscous sublayer at any
vehicle Reynolds number. Wall shear is therefore under-resolved.

This is tolerable for bluff bodies, where drag is pressure-dominated
(separation is fixed by body-line breaks, not by a delicately balanced
boundary layer). It is **not** tolerable for a streamlined body, a wing at
low angle of attack, or any case where friction is a large share of total
drag. The Blasius case is the honest measure of this limitation: with the
boundary layer properly resolved (20–37 points), cf is still +4 to +8 %.

### 6.4 Uniform grid: cannot refine near the body

A direct consequence of the FFT pressure solve. Resolution near the body
cannot be raised without raising it everywhere, at O(N³) cost. This is the
structural limit on accuracy per unit runtime, and it is the first thing to
change if much higher fidelity is ever needed (at the price of the 30-minute
budget).

### 6.5 Direct-forcing produces a 2Δt force oscillation

The instantaneous IBM force alternates between two values every step, at a
few percent amplitude. Diagnosed explicitly: 398 of 399 consecutive
differences flipped sign.

It does **not** bias the mean (averaging cancels it exactly), but before it
was handled it corrupted two derived quantities: the reported force standard
deviation, and the Strouhal number, where the FFT locked onto the Nyquist
spike and returned St ≈ 57 instead of ≈ 0.16. Now removed by a two-point
boxcar, which has an exact zero at Nyquist and leaves the shedding band
untouched.

### 6.6 Boundary conditions that look equivalent and are not

Two failures found during validation, both worth stating because both
produced plausible-looking wrong answers rather than obvious errors:

- **`slip` as an open boundary.** A slip plane is impermeable. A growing
  boundary layer displaces flow that then has nowhere to go, driving u above
  U∞. On the flat plate this flipped the *sign* of the computed momentum
  thickness (−221 %), because the near-wall deficit was cancelled by the
  excess above it. Fixed by adding a `farfield` type that lets the normal
  velocity pass.
- **Re-applying boundary values after the projection.** The projection takes
  boundary-normal velocity as given. Overwriting it afterwards silently
  reintroduces divergence in the outermost cell layer: the residual jumped
  from 2e-16 to 9e-5 while the interior still looked clean. Fixed by splitting
  the BC layer into `set_boundary_values()` (pre-projection only) and
  `fill_ghosts()` (safe any time).

### 6.7 Reynolds-number ceiling

The verification cases are laminar and fully resolved. A vehicle at Re ≈ 10⁷
on a few million cells is not: the resolved scales stop well short of the
dissipative range. Without a subgrid model, dissipation comes from truncation
error (implicit LES). It works, but the amount is an accident of the scheme
rather than a stated assumption. A Vreman model is available
(`waterfall_core/sgs.py`) and is the right choice above Re ≈ 10⁵.

Vreman rather than Smagorinsky because constant-coefficient Smagorinsky does
not vanish at a wall or in laminar shear, so it damps the attached boundary
layer it should leave alone. Known approximation: the subgrid stress is
applied as div(ν_t ∇u), dropping the transpose term.

### 6.8 What is *not* a solver defect

- 2D results at Re ≳ 190 overpredict Cd and Cl_rms because the real wake
  becomes three-dimensional there. That is a property of running a 2D case,
  not an error in the discretisation.
- A stationary floor under a moving car thickens the floor boundary layer and
  inflates drag. Use `rolling_road: true` (the default); a fixed floor is a
  wind-tunnel artefact, not reality.

---

## 7. Runtime

Measured on the Ryzen 9 5900X (12C/24T). No GPU path: the RX 6800 is AMD, so
CUDA/CuPy is unavailable; everything runs on 24 CPU threads via Numba.

| Case | Cells | Steps | Wall time | ns/cell/step |
|---|---|---|---|---|
| Blasius (2D) | 115 k | 12,367 | 3.6 min | ~150 |
| Cylinder Re=40 (2D) | 346 k | 12,700 | 8.0 min | ~109 |
| Cylinder Re=100 (2D) | 346 k | 19,900 | 10.9 min | ~95 |
| Sphere Re=100 (3D) | 432 k | 1,234 | 1.1 min | 127 |
| Tesla Model S (3D) | 1.05 M | ~21,500 | ~38 min | ~105 |

**Where the time goes** (2D, 230 k cells, per step):

| Stage | Before tuning | After |
|---|---|---|
| Pressure projection | 15.0 ms | 8.6 ms |
| Ghost-cell fill | 6.3 ms | 0.13 ms |
| Convection + diffusion | 0.6 ms | 0.6 ms |
| **Total step** | **40.8 ms** | **19.1 ms** |

Two changes did most of that. The DCT now transforms only non-degenerate axes,
since a transform along a length-1 axis is the identity but scipy still walks the
array for it. And in 2D the z-derivatives are switched off in the kernels
(`zfac = 0`), which makes the z ghost cells unreadable and lets the boundary
layer skip them: with `nz = 1` a z-face plane spans the *entire* domain and is
a stride-3 write, so filling it cost more than convection and diffusion
combined.

### Against the 30-minute target

The Model S case at 40 cells/length overshoots at ~38 min. The runtime
estimator was the thing at fault, and it has been corrected: it assumed peak
local velocity of 1.4×U∞, where the measured timestep implies ~2.5×U∞. Sharp
immersed features accelerate flow far more than smooth-body intuition
suggests, and CFL is set by the single worst cell, not a representative one.

Runtime scales roughly as (cells per length)⁴: finer cells cost more cells
*and* a smaller timestep. Dropping 40 → 36 therefore buys about a third of the
runtime back and lands near 27 min. `cases/tesla.yaml` ships at the setting
that meets the budget, with the higher-resolution result documented alongside.