FluidSimCpp/verification/FLUID_PHYSICS_AUDIT.md

# Fluid physics audit — 23 October 2026

This pass repairs conservation and boundary coupling, establishes repeatable
physical tests, and measures kernel and threading performance. It is not a
certification that every engine or acoustic subsystem is physically validated.

## Findings and changes

1. **Pipe–volume mass/energy mismatch.** Previously, the volume used an orifice's
   requested discharge while the pipe used a different Rusanov ghost-cell flux.
   Each boundary now supplies one integrated mass, energy and fresh-air flux;
   the volume receives the exact opposite exchange. Stationary walls transfer
   momentum but no mass or energy. Open-end audio reports the mass flux actually
   applied by the finite-volume update.
2. **Substep integration.** Components previously received only the final
   boundary flow after several pipe substeps. Every component now updates after
   every pipe substep. CFL is recalculated as the solution evolves.
3. **Profiled pipes.** Geometry was stored but absent from the conservation
   update. Shared face fluxes now include face area and divide by cell volume;
   momentum includes the quasi-1D pressure source `p dA/dx`. Constant-pressure
   rest remains stationary while closed profiled pipes conserve mass/energy.
4. **Orifice choking.** The old formula used the upstream sound speed with
   throat Mach number, overestimating discharge. The nozzle calculation now
   uses the isentropic throat temperature and satisfies the sonic-flow test.
5. **Unstable small-signal boundary.** Evaluating `sqrt(dp)` using only the
   adjacent cell-centre pressure destabilized a small Helmholtz excitation.
   The nozzle discharge is now matched to the outgoing duct characteristic
   using a safeguarded analytic Newton solve of a squared-flow residual.
   Bisection remains the fallback. Subsonic, isentropic characteristic and
   quasi-steady nozzle assumptions still limit this model.
6. **Inertance units.** For mass flow, `d(mdot)/dt = A dp/L`, not
   `A dp/(rho L)`. Dimensionless quadratic resistance is integrated with a
   dissipative semi-implicit term. Removed the additional lumped inertance in
   the Helmholtz scenario: its resolved neck already carries inertance.
7. **Moving volumes.** Volume0D previously applied work without changing its
   geometric volume. Closed adiabatic work now follows `E V^(gamma-1)=constant`
   and updates volume. Crankcase geometry work is applied once per geometry
   change rather than again for each solver substep. Removed its five-bar cap,
   which silently deleted gas during high compression.
8. **Invalid configurations.** Reject invalid geometry, timesteps, duplicate or
   missing endpoints and incompatible gas constants. A pressure-only boundary
   cannot specify a supersonic inlet; that case now fails explicitly. Direct
   pipe joints support equal endpoint areas and all orientation combinations;
   unequal-area direct joints are rejected. Use a profiled transition instead.
9. **Performance.** Cache primitive states/wave speeds, evaluate each shared
   interior face once, and hoist uniform relaxation outside the cell loop.
   Compiler SIMD is enabled without fast-math. No per-step allocations were
   introduced into the pipe/boundary kernels.

## Validation

`fluidsim_fluid_validation` contains 29 named cases, in addition to the existing
physics/audio/ECU regression executable. Coverage includes:

- Closed pipes, variable-area equilibrium/conservation, stationary hot contacts.
- Sealed volume–pipe exchange, several volumes and a shared tiny donor.
- Direct joints in four orientations and split/unsplit equivalence.
- Fresh-air species conservation, actual open-end mass accounting on both ends.
- Sonic choking against the analytic ideal-gas formula; inertance units and
  steady quadratic resistance.
- Adiabatic moving volumes, thermal exponential relaxation and sealed crankcase
  compression beyond five bar without mass deletion.
- Closed/open acoustic reflection signs, closed-duct mode frequencies at
  50/100/200 cells, and small-signal cavity resonance.
- Sod shock-tube density error at 100/200/400 cells, strong pressure pulses,
  expanding-flow near-vacuum states, positivity-repair checks.
- Invalid geometry/configuration, incompatible gas coupling, unsupported
  supersonic inlet and missing/duplicate endpoint rejection.

Representative measured results:

| Check | Result |
| --- | --- |
| Sod density L1, 100 / 200 / 400 cells | 0.02332 / 0.01595 / 0.01052 |
| Closed-duct fundamental | approximately 171.9 Hz; analytic 171.91 Hz |
| Cavity fundamental | approximately 136.2 Hz; lumped estimate 138.51 Hz |
| Closed mass/energy relative tolerances | 2e-5 to 5e-5, depending on case |
| Cavity frequency tolerance | 5%, with growth and positivity checks |
| Closed-duct frequency tolerance | 1% across three resolutions |

The cavity estimate assumes a compact, lossless neck without radiation end
correction; the discretized model includes neck compressibility and numerical
losses. These comparisons are verification of stated assumptions, not a
comparison with measured engine hardware.

## Benchmarks and SIMD

Machine: Intel Core i5-12600K, GCC 16.2.1, Linux x86-64, Release. These are local
wall-clock microbenchmarks; CPU frequency and background load affect timings.
The pipe benchmark advances an equilibrium lossless duct with reflected ends;
the twin benchmark runs one simulated second including its nonlinear boundaries.

Before this audit, isolated kernels measured about 49–52 million cell updates/s
and the twin about 3.84 times real time. Cached shared-face evaluation raised
pipe throughput to roughly 110–140 million updates/s. The final twin is around
1.7 times real time: it does substantially more boundary work than the original,
physically inconsistent solver. Early faster twin numbers from before the
characteristic/nozzle correction must not be treated as final results.

A SIMD comparison measured roughly 111–139 million updates/s with compiler
vectorization versus 84–102 million with vectorization disabled. GCC's
vectorization report confirms 16-byte SIMD in the primitive-state and face-flux
loops. `FLUIDSIM_ENABLE_SIMD=ON` permits vectorization and disables math errno;
it does **not** enable fast-math, approximate reciprocal/square-root intrinsics,
or machine-specific instruction sets. OFF disables compiler vectorization for
repeatable comparison. Scalar and SIMD builds run the same physics tests.
Windows/MSVC settings are provided but were not executed in this Linux audit.

## Multithreading investigation

`fluidsim_thread_benchmark` measures persistent workers with two barriers per
iteration. Boundaries/components must remain in ordered phases; pipes can be
advanced independently only after all boundary fluxes are resolved.

| Independent pipe workload, 1000 steps | Serial | 2 workers | 4 workers |
| --- | ---: | ---: | ---: |
| 6 pipes × 10 cells | 0.00093 s | 0.00272 s | 0.00625 s |
| 64 pipes × 256 cells | 0.1275 s | 0.0700 s | 0.0410 s |

This is an optimistic stage-only benchmark with equilibrium data, excluding
worker startup and a real graph's serial boundary work. It supports further
work on persistent pools for large networks; it does not justify threading the
current small twin. No fluid worker pool is enabled in production. The existing
simulation/audio/UI separation remains.

## Reproduce

```sh
cmake -S . -B build-audit -DFLUIDSIM_BUILD_DESKTOP=OFF -DCMAKE_BUILD_TYPE=Release
cmake --build build-audit -j 4
ctest --test-dir build-audit --output-on-failure
./build-audit/fluidsim_physics_benchmark 44100
./build-audit/fluidsim_thread_benchmark
```

For a scalar comparison configure a separate directory with
`-DFLUIDSIM_ENABLE_SIMD=OFF`. On GCC/Clang, use another directory with
`-DCMAKE_CXX_FLAGS="-fsanitize=address,undefined -fno-omit-frame-pointer"` to
reproduce the sanitizer run. Final Linux desktop regressions, all 29 scalar fluid cases, and both test
executables under AddressSanitizer/UndefinedBehaviorSanitizer passed. The final
sanitizer run completed in about 13 seconds.

Compiler vectorization diagnostics can be produced
with GCC `-fopt-info-vec-optimized`.

## Remaining limits and next work

- First-order Rusanov flux is dissipative. Spatial/time resolution convergence
  and attenuation across the full audible band are not yet established. A
  higher-order limited scheme is a separate change, not silently enabled here.
- The default four-substep floor was retained: removing it failed a tiny-volume
  stress case. CFL alone is not a positivity guarantee for arbitrary stiff
  pipe/volume networks. Conservative donor budgets can limit transfer; inspect
  `BoundarySystem::limited_fluxes` and refine time resolution when they activate.
  Budgets for shared donors are allocated in registration order. Pipe positivity
  repairs are counted by `Pipe1D::positivity_repairs`; a nonzero count invalidates
  a claim of exact conservation for that run.
- Pipe geometry/cell storage remains public for initialization. Changing its
  topology or resizing cells after connections have been built is unsupported.
- The pressure-only open boundary extrapolates interior entropy. It is an ideal
  pressure-release model, not a complete external atmosphere, anechoic boundary,
  or frequency-dependent radiation impedance. Strong inlet flows need richer
  thermodynamic boundary data.
- Junctions and orifices are quasi-1D approximations. Jet separation, vena
  contracta, discharge coefficients and strong multidimensional area changes
  need validation; equal-area joint tests do not validate arbitrary fittings.
- The optional lumped inertance law is dimensionally checked, but its stored
  neck kinetic energy is not a separately tracked thermodynamic state. Use
  resolved pipe cells when evaluating physical neck acoustics and energy.
- Core pipe damping and volume thermal relaxation default to zero. Explicit
  scenario damping/relaxation still represents empirical source terms and must
  be calibrated before claiming real dissipation or heat transfer.
- Engine cylinders retain pressure/temperature safety clamps, approximate
  combustion/species bookkeeping and unvalidated mechanical work integration.
  ECU behavior and the broader engine model are not certified by this audit.
- Sound radiation remains the documented compact-source approximation. These
  fluid tests do not establish microphone SPL, high-frequency radiation,
  thermoviscous losses or full mechanical/chemical energy conservation.

## Physical references

- [Clawpack: approximate Euler Riemann solvers](https://www.clawpack.org/riemann_book/html/Euler_approximate.html)
  — conservative Euler fluxes and total-energy equation of state.
- [Clawpack: Euler equations and exact Riemann solutions](https://www.clawpack.org/riemann_book/html/Euler.html)
  — shock/contact/rarefaction verification and the Sod problem.
- [NASA Glenn: mass-flow choking](https://www.grc.nasa.gov/www/k-12/BGP/mflchk.html)
  — isentropic throat relations and maximum mass flow at Mach one.
- [EPFL: quasi-1D Laval nozzle model](https://www.epfl.ch/labs/mathicse/wp-content/uploads/2018/10/25.2016_MP-FN-PL.pdf)
  — variable-area Euler fluxes and geometric pressure source.


## Sustained-load follow-up

A one-second run missed a runtime deadline problem: at full throttle the heated
gas raised average CFL work from 4 to approximately 7 substeps per audio frame.
By simulated second 5, each second needed about 1.04 seconds of CPU time, so a
finite startup buffer could not prevent eventual underruns. Idle stayed near
0.6 seconds per simulated second.

The solver now reuses its CFL primitive-state/wave-speed preparation in the
finite-volume update. Boundary Newton residual tolerance is matched to the
single-precision conserved state (`4e-7` in the normalized squared-flow
residual), rather than solving to `1e-12` relative residual each time. The same
CFL limits and four-substep floor remain. On the same machine, a nine-second
full-throttle run measured approximately 0.74–0.81 seconds per simulated second
in the hot portion; idle measured approximately 0.49–0.52 seconds. All 29 fluid
cases and existing regressions still pass. This is measured headroom, not a
hard real-time guarantee on other workloads or machines.

Reproduce sustained-load timings with:

```sh
./build-audit/fluidsim_physics_benchmark --soak 12
```

The benchmark reports each simulated second, average substeps/frame and
boundary iterations/frame. Desktop audio also counts underrun blocks; the
counter is shown with the audio-buffer status when a UI font is available.

## Intake reservoir and repeated-blip regression

The twin's upstream pressure-only boundary inherited pipe entropy during
inflow. Repeated blips could therefore admit progressively hotter gas rather
than fresh ambient air: a captured sequence reached approximately 1140 K in
the plenum at wide-open throttle and lost fresh charge density. The intake now
connects to an explicit infinite `Atmosphere` reservoir (101325 Pa, 300 K,
fresh-air fraction one), updated on every fluid substep. Intake radiation uses
that reservoir's actual signed mass exchange. The remaining pressure-release
open-end index now refers only to the exhaust.

The twin also disables pipe relaxation toward ambient **pressure** energy.
That source can heat low-density evacuated cells without a physical wall heat
transfer law. This does not impose any pressure floor or artificial manifold
refill. Hot residual-gas reversion, closed-throttle pumping and pressure waves
remain possible; their engine-model accuracy still needs further validation.

`fluidsim_intake_regression` runs 16 simulated seconds: idle, an 80% blip,
recovery, a full-throttle opening, and return to idle. It checks pressure
recovery, fresh-air temperature/composition, finite audio, zero pipe positivity
repairs and final idle recovery. Representative results after the fix:

- Before the second opening: 35.76 kPa manifold absolute pressure.
- After opening: 85.45 kPa near 9950 RPM; manifold 298.96 K, inlet 296.92 K.
- After return to idle: 34.51 kPa at 1820 RPM.

Manifold vacuum at small/closed throttle is expected and is not removed by the
fix. An actual obstruction/regression is distinguished by failure to recover
pressure and fresh-charge density when reopening the throttle. The inlet
reservoir controls the upstream gas state; it does not thermostat the manifold.
See [AEM's engine-management manual](https://documents.aemelectronics.com/techlibrary_30-6905_universal_programmable_engine_management_system.pdf)
for the qualitative throttle-angle/manifold-pressure relationship. These are
scenario regression bounds, not universal pressure limits for real engines.

All three CTest executables pass after this change (existing regressions,
29-case fluid suite, and the new intake transient). A six-second sustained
full-throttle follow-up stayed at four substeps/frame and approximately
0.51–0.53 seconds of wall time per simulated second, with RPM around the ECU
rev limit. This also identifies the earlier hot-gas CFL growth as a consequence
of the incorrect inlet thermodynamic state, not solely an unavoidable engine
workload increase. The earlier timing sections are historical measurements of
those intermediate implementations.