Engine design: config integrity, intent/output split, and chug stability fixes

[swissskimmilk]

Scope of this writeup

This branch carries two independent bodies of work. This description covers the config-integrity
and chug-stability changes only. The thermal / Bartz / combustion-efficiency commits on the same
branch are a separate effort and will be documented in a follow-up section.

---

1. Chug stability — the time-series path never saw a real injector pressure drop

Issue. The array / time-series evaluation path (Layer 2) built its stability-diagnostics dict
with only mdot_O/F and the tank pressures — the injector-ΔP keys were never forwarded. So the chug
model fell through to its or 0.30 * Pc fallback on every timestep of every evaluation, scoring
stability against an assumed 30%-stiff injector that never reflected the actual design. Two related
fail-open masks compounded it:

• _chug_gate_margin / _acoustic_gate_margin returned a neutral-pass 1.10 on non-finite input,
  so a failed root-find silently passed the gate.
• the injector-ΔP fallback substituted a healthy value on missing data — and because it used or,
  a genuine ΔP = 0.0 (a completely unstiffened injector, the most unstable case possible) also
  read as 30%.

Fix.

• Forward delta_p_injector_O/F and delta_p_feed_O/F into the array-path stability call.
• Injector ΔP now fails closed: missing/non-finite → NaN → treated as unsafe, never an
  assumed value. Feed ΔP (a secondary term) routes through the assumptions registry so its fallback
  is declared, not hidden in a literal.
• The gate-margin helpers fail closed on non-finite input.
• Surface injector stiffness η = ΔP_inj / Pc as a first-class per-timestep result.

Evidence. η = ΔP_inj/Pc is the standard, calibration-free chug criterion:

• MIT Rocket Team, Injector Design: "a good rule of thumb … that the drop in pressure across an
  orifice be 20% of the chamber pressure … also referred to as the stiffness of the injector."
• Huzel & Huang / Sutton design range ΔP_inj/Pc ≈ 15–25% of Pc; NASA throttling review (NTRS
  20090037061) gives the same band and the "prevents communication between chamber gases and the
  feed system" rationale.
• Confirmed the underlying physics against this repo's own model: since ΔP ∝ ṁ² and Pc ∝ ṁ, stiffness
  should scale as η ∝ ṁ. Measured η/ṁ is constant to three significant figures across a pressure
  sweep (0.2283 / 0.2291 / 0.2301), and stability_source: native confirms the C kernel is the live
  gate.

This removes the same fail-open class addressed elsewhere on the branch. Layer 1's single-point path
already forwarded the real ΔP, so this change does not alter Layer 1 behaviour.

---

2. Config integrity — staleness was structurally invisible

Issue. Configs mixed hand-authored intent (design_requirements) with optimizer-produced
output (chamber_geometry, injector geometry, tank pressures, pressure curves) in one file, with
nothing marking which half the optimizer owns. Consequences found in the audit:

• No committed config satisfied its own declared gates, and nothing reported it. Canonical's
  geometry was solved for 7000 N / MR 2.55 while its requirement asked for 8000 N / MR 2.8 — the
  geometry block literally records a different spec than the requirements. Root cause: requirements
  edited without re-running Layer 1, plus the physics moving under frozen geometry.
• The backend booted on default.yaml, which sets no propellant preset → fluid identity None
  → non-finite solver residual → negative thrust (F = −334 N, Isp = −7.27 s). A fresh session
  evaluated a physically impossible engine until the user loaded something else.

Fix.

• tests/test_config_design_validity.py: loads every config, evaluates it at its own design point.
  Fails only on genuine breakage (won't load / won't evaluate / non-physical). Off-spec configs
  (miss their own target) are reported, not gated — gating them would block work-in-progress
  designs and create pressure to widen bands, the exact anti-pattern that produced the current
  bands. A KNOWN_BROKEN ratchet that can only shrink.
• engine/pipeline/run_record.py + scripts/record_runs.py: emit one JSONL row per (commit, config)
  with git provenance, a resolved-config sha256, and F/Isp/Pc/MR/η metrics. config_sha
  unchanged + git_sha changed ⇒ the physics moved; the reverse ⇒ someone edited a config. Makes
  regressions attributable to a commit; intended for a CI-written metrics branch (wiring authored,
  not yet switched on).
• Backend boots on configs/canonical/impinging.yaml — already the documented source of truth via
  POST /api/config/switch — instead of the unmaintained default.yaml.

Evidence. The guard caught a mislisted config on first run; the recorder immediately surfaced that
the tanh gate saturates in real data (two engines with GM 3.15 vs 5.75 both reporting a squashed
1.2999…). Negative-thrust default reproduced (F = −334 N) and the canonical boot verified clean.

---

3. Config structural cleanup — separate what a human types from what the code writes

Four behavior-preserving refactors, each verified by a 0-diff on the fully-resolved config vs its
pre-change version (fields silently reconciled at load are compared, not just the raw YAML):

• material_preset (configs/materials/phenolic_graphite.yaml): material constants (heat of
  ablation, densities, oxidation kinetics) are properties of the substance, not the engine, but were
  duplicated verbatim — 61 identical fields across the two canonical configs alone. Now shared
  via the same mechanism as the existing propellant_preset; explicit-wins, every override logged.
• Intent / output split: configs/canonical/<name>.yaml holds intent; <name>.outputs.yaml
  holds optimizer output. Merged at load, re-split on save; hand-setting an output field in the
  intent file is a hard load error. Membership is checkable — io._OUTPUT_FIELDS lists exactly the
  fields something in the codebase assigns.
• Burn time was stored in three slots (thrust.burn_time, pressure_curves.target_burn_time_s,
  design_requirements.target_burn_time) force-synced to the requirement, so two of the three file
  values were dead. Now single-sourced to the intent field; the synced copies live in the outputs.
• Design-point stamp (chamber_geometry.design_thrust/MR/pressure) moved to the outputs sidecar,
  where "solved for 7000 N vs want 8000 N" reads as a legible staleness signal instead of a confusing
  duplicate in the intent block. geometry.py now reads the required ChamberGeometryConfig fields
  directly instead of getattr(..., <magic default>) — the defaults were dead (required Field(gt=0)
  fields) and fail-open.

Evidence. Full suite green throughout (360 passed / 75 skipped at the tip); new tests lock each
invariant (preset resolution + a no-reinline guard, split partitioning, collision-is-error, and a
save round-trip that proves the API can't collapse a split config back into one file).

---

Reviewer notes

• Behavior-preserving is verified for every config refactor by comparing fully-resolved configs
  (post preset-overlay and post-sidecar-merge), not raw YAML — the resolved config the solver sees is
  byte-identical before and after.
• No Layer 1 optimization logic is touched. The chug fix affects the Layer 2 array path only; the
  single-point path already forwarded real ΔP.
• Deferred, logged (not in this branch): having Layer 1 write the design stamp on solve + a
  staleness check comparing intent-target vs output-stamp; an independent achieved-vs-target burn
  time for blowdown; the chug-gate re-centring (blocked on hardware measurement).
• Uncommitted native/thermal working-tree changes are not part of this push.

Thermal / gas-side heat transfer (review item 1.4)

Follow-up to the config/chug writeup above. This covers the thermal, Bartz, and
combustion-efficiency commits on the branch. Review item 1.4 was "heat load
~3–5× overestimated, and the checks that matter are dead." The commits below
work through the chain that produces the gas-side heat load, from transport
properties up to the gas-side film coefficient. All Huzel references are to
Huzel & Huang, Design of Liquid Propellant Rocket Engines, NASA SP-125, ch. IV.

Net effect on the canonical impinging chamber (throat/chamber heat flux):
convective flux went from a collapsed 0.009 MW/m² to a physical 4.36
MW/m², radiative flux from an inflated 5.97 to 0.28 MW/m², and the two
happened to sum to a plausible-looking total the whole time — which is why none
of it was caught (see §2).

---

1. Gas-side transport properties — three contradicting sources and a 386× unit error

Issue. The Bartz/Dittus-Boelter group needs four properties of the combustion
products (μ, cp, k, Pr). They were computed in three separate places that had
drifted into disagreeing by ~3× in both k and Pr, in opposite directions, and
one produced Pr = 2.25 — thermodynamically impossible for a gas. All three
derived cp as γR/(γ−1) using CEA's GAMMAs, which is the isentropic exponent
of a dissociating mixture, not cp/cv — overshooting cp by ~40%. Worst of all,
the viscosity conversion multiplied Huzel eq. (4-16)'s output by 6894.76
(a psi→Pa factor) when eq. (4-16) returns viscosity in lb/(in·s), not
lbf·s/in² — collapsing the convective coefficient by exactly g_c = 386.088.

Fix.

• One provider, thermal/gas_transport.hot_gas_transport, replaces the three
  copies. Rule: none of μ/cp/k/Pr is derived from the others — that
  derivation is what let the copies disagree.
• μ/cp/k/Pr taken from CEA's frozen transport block per propellant (LOX/CH₄:
  k=0.394, Pr=0.644, cp=2463), not derived. Frozen, not equilibrium — Bartz's
  eq. (4-16) is a diatomic-air fit with no reaction contribution, so pairing it
  with equilibrium properties double-counts recombination (~2.95× on the group).
• Viscosity constant corrected to LB_PER_IN_S_TO_PA_S = 17.858 and renamed
  (the old name/value was correct for lbf·s/in², merely misapplied). Mirrored in
  the C kernel (ED_LB_PER_IN_S_TO_PA_S).

Evidence. Huzel eq. (4-16), μ = 46.6e-10·M^0.5·T^0.6, with the nomenclature
under eq. (4-12) giving units of lb/(in·s). Sample Calculation 4-3 carries
μ = 4.18e-6 lb/in-sec, which the corrected code reproduces. The 386.088 ratio
is exactly g_c in lbm·in/(lbf·s²), confirming the two unit systems that were
conflated. CEA's own header labels conductivity MILLICALORIES/(CM)(K)(SEC)
(factor 0.41842 → W/m·K), caught via a Pr = μcp/k consistency check.

---

2. The broken values were invisible — silent dict-key typos and unrendered fields

Issue. Two .get() calls in the time-varying solver read keys that never
existed ("h_g" vs the produced "h_hot", "heat_flux_radiative" vs
"heat_flux_rad"), so they silently returned hardcoded defaults and discarded
the real computed values. Combined with the fact that the only heat-flux number
surfaced to the UI is the conv+rad total, the 386× viscosity error was
undetectable: collapsed convection (0.15% of total) was almost exactly offset by
inflated radiation, so the total read a plausible ~6 MW/m². h_g, throat flux,
and surface temperature are rendered nowhere, and zero tests imported any
thermal module.

Fix. Both keys now direct-index (dict["h_hot"]), so a producer/consumer
rename raises instead of silently defaulting. Added tests/test_gas_side_heat_transfer.py
as the first test to exercise the thermal modules directly.

Evidence. With the radiative key fixed (previously always 0), the degenerate
wall solver in §3 produced −16,650 K, which is what forced §3 — a latent bug
the typo had been masking.

---

3. The wall-temperature solver returned its own initial guess

Issue. calculate_steady_state_temperature_profile iterated
Ts_new = Tg − (h(Tg − Ts) + q_rad)·(1/h), which algebraically simplifies to
Ts_new = Ts − q_rad/h — an identity in Ts that never references the wall
conduction stack. With q_rad = 0 (its state until §2) it "converged" instantly
on its initial guess Ts = 0.8·Tg for any h, conductivity, or thickness;
with q_rad > 0 it marched down by q_rad/h for a fixed 10 iterations,
reaching −16,650 K.

Fix. Replaced with a hot-surface energy balance solved by bisection:

h(Tg − Ts) + q_rad,in  =  εσ(Ts⁴ − T_amb⁴) + (Ts − T_amb)/R_wall

The residual is strictly decreasing in Ts, so bisection is unconditionally safe —
no initial guess to get wrong, no iteration-count artefact.

Evidence. Validated against Huzel Sample Calculation 4-7 (radiation-cooled
wall, eq. 4-38): ours 2666 °R vs Huzel's printed 2660 °R (0.24%), heat flux 259.4
vs 260.0 kW/m². The surface re-radiation term flagged as "a modelling choice" is
in fact Huzel's own eq. (4-38).

---

4. Gas radiation used the charred wall's emissivity, not the gas's

Issue. The gas→wall radiation term q_rad = ε·σ·(T_gas⁴ − T_wall⁴) used
ε = ablative_config.surface_emissivity — the charred wall's 0.85 — where
the physically correct value is the gas emissivity. With ε=0.85 the model put
radiation at ~39% of total wall heat flux, roughly 5× too high.

Fix.

• q_rad now uses a gas emissivity eps_gas, carried on the hot_gas_transport
  dict so the lumped and axial paths cannot diverge and cannot reach for a wall
  surface property by mistake.
• DEFAULT_GAS_EMISSIVITY_ND = 0.10, applied across the config schema and 8 YAMLs
  (was 0.8).

Evidence. Göbel, Kniesner, Frey, Knab & Mundt (EUCASS 2013 / CEAS Space
Journal 6:79, 2014): for CH₄/O₂, the radiative share of total wall heat flux is
~8% at its local peak and ~2.5% integrated; H₂/O₂ 9–10%. Our ε=0.10 maps to a
7.0% radiation share at the canonical operating point (rule of thumb: share ≈
0.7·ε here), i.e. near the local peak — defensible for a lumped chamber number
and within the cited band. Stated plainly in the code: 0.10 is an engineering
estimate, not a derived value. Huzel carries no gas-emissivity term at all;
he instead notes (p. 101) that eq. (4-13) under-predicts "if a substantial
fraction of the combustion gases are strong radiators" — i.e. the term we add is
correcting a deficiency he flags, not contradicting him. Net effect: canonical
total flux 10.21 → 6.76 MW/m², radiation share 38% → 6.8%.

Related — kerosene carbon-deposit warning. For LO₂/RP-1 Huzel eqs.
(4-17)/(4-18) model a carbon deposit as a thermal resistance R_d in series with
the gas film (fig. 4-25: R_d dominates, cutting gas-side conductance ~5×). We do
not model it — the published data is a single curve at one propellant/Pc/MR, and
extrapolating would invent precision. Instead warn_if_carbon_depositing fires
once per kerosene-class fuel when a wall model is enabled, stating that heat load
is over-predicted and liner sizing is unreliable for that propellant.

---

5. Combustion inefficiency never reached the chamber temperature

Issue. CEA reports the chamber temperature of a perfect chamber. A real
chamber releases less energy (incomplete mixing/kinetics/finite L*) and runs
cooler. The code applied the combustion efficiency η to c* (cstar_actual = η· cstar_ideal) but never to the temperature — the thermal chain consumed CEA's
raw Tc_ideal. Heat flux is steeply temperature-dependent, so this over-predicted
the heat load. An orphaned calculate_actual_chamber_temp existed but (a) was
called only by a plotting script, and (b) used an uncited, underivable formula
that disagreed with Huzel by ~9% on the multiplier.

Fix.

• calculate_actual_chamber_temp rewritten as Tc · η_c*², and the thermal
  chain now runs on the corrected temperature. Requires a solver reorder:
  combustion efficiency is computed before the cooling models (it depends
  only on geometry/Pc/mixing/kinetics — verified non-circular by sweeping the
  one downstream input, Tc_kinetics, over 1000–4500 K for zero change in η).
• eta_cstar split into its combustion and cooling components so a temperature
  correction consumes combustion-only η, never the cooling-blended value (which
  would deduct cooling from temperature twice — it already reaches temperature
  via effective_Tc). Mirrored in the C kernel.

Evidence. Huzel Sample Calculation 4-3: design (Tc)ns = theoretical (Tc)ns × (c* correction factor)² = 6460 × 0.975² = 6140 °R. The exponent is not fitted:
c* = √(γRTc)/Γ, so at fixed γ, R an efficiency measured on c* lands on
temperature squared. Measured on canonical: Tc_ideal 3013 K → Tc_combustion
2683 K (−336 K, −11.2%), worth ≈ −14% convective and −38% radiative flux once
wired through (radiation goes as T⁴). Layer 1 is bit-identical (it disables
cooling, so the reorder changes nothing there — measured, not asserted).

---

6. Cooling's energy deficit was applied to c* linearly, not via √

Issue. The cooling efficiency is an energy fraction,
1 − Q_removed/(ṁ·cp·Tc). It was multiplied straight into the c* efficiency. But
c* ∝ √Tc, so an energy deficit of f costs c* only √(1−f), not (1−f) —
roughly doubling cooling's penalty on c*.

Fix. Cooling's contribution to the c* efficiency is now √(cooling_eff);
the energy fraction is still applied directly to temperature (via effective_Tc),
which is its correct currency. Both currencies are exposed on the components dict.
Mirrored in the C kernel.

Evidence. Verified against the CEA cache: reconstructing cstar_ideal from
CEA's own Tc/γ/R reproduces it to 0.07% mean over MR 2.4–4.0, and scaling Tc by
f scales c* by √f to five decimals. At the canonical cooling_eff = 0.96,
the c* penalty drops from 3.9% to 2.0% → Pc +1.15%, c* +2.1%, Isp +2.3%.

---

7. Gas-side film coefficient: Dittus-Boelter → Bartz

Issue. The gas-side coefficient was plain Dittus-Boelter, 0.023·Re^0.8·Pr^0.4
on chamber diameter — a smooth-pipe, moderate-ΔT, constant-section correlation.
A thrust chamber violates every one of those: Twg/(Tc)ns ≈ 0.35, accelerating
flow, varying cross-section. It has no throat concept, so it cannot produce the
peak flux at the throat where liners actually fail.

Fix. New thermal/bartz.py implementing Huzel eq. (4-13) with the eq. (4-14)
σ property-variation correction, exposed as one provider bartz_chamber_h_g that
all cooling models share (mirroring the single-provider pattern from §1). Wired
into the recession-driving path and mirrored expression-for-expression in the C
kernel. The throat radius of curvature R derives from the nozzle's own Rao arc
coefficients (0.941·R_t, the mean of the 1.5 and 0.382 fillets), coupled to the
drawn contour by a circle-fit test so the two can't drift.

Evidence.

• Term-by-term against Huzel Sample Calculation 4-3. Every printed
  intermediate reproduces: 0.026/D_t^0.2 = 0.01366, μ^0.2·cp/Pr^0.6 = 0.046,
  (Pc·g/c*)^0.8 = 4.02, (D_t/R)^0.1 = 1.078, product 0.0027, and all three
  station h_g values (0.00185 chamber / 0.0027 throat / 0.000507 exit).
• σ is eq. (4-14) verbatim (both ½ factors, 0.68/0.12 exponents), pinned by
  an independent re-transcription and cross-checked against fig. 4-24 at the three
  points Sample Calc 4-3 reads off (1.6/3.1/2.5% — the chart being the imprecise
  side, read to two significant figures).
• Pr exponent kept at 0.4, deliberately. Huzel offers eq. (4-12) (Colburn,
  Pr^0.34) and eq. (4-13) (Bartz) as alternatives ("or as Bartz has shown").
  Writing (4-13)'s Nusselt group with k = μcp/Pr gives h ∝ μ^0.2·cp·Pr^(n−1),
  so Bartz's Pr^−0.6 is n = 0.4. A test locks this against the tempting but
  wrong "fix" to 0.34.
• Effect vs Dittus-Boelter at the canonical chamber: −15.7% on h_g. The sign
  flips at the throat, where Bartz peaks and D-B has no throat term — the
  correct axial variation is the actual point of the change.
• Wiring exposed the latent cp divergence flagged in §1: the C kernel still
  derived cp = γR/(γ−1) (dead under Dittus-Boelter, which used cp only via the
  overridden Pr). Bartz uses cp directly, so hot_gas_cp was added to the
  EdCooling struct and the C kernel now reads the config value.

---

8. Test infrastructure — two silent-green failure modes closed

Issue. (a) The A/B parity suite's "Python reference" was runner.evaluate()
with the native kernel on by default — so for Pc, F, and Isp (solved in
C and copied into diagnostics) it was comparing C against itself, bit-
identically, incapable of failing. (b) ctest silently passed residual_golden
on a stale test binary: rebuilding only the shared library (what A/B parity
loads) left the ctest .exe files untouched, so old-binary-vs-old-golden agreed.

Fix.

• Parity reference is now computed with the native path explicitly disabled, so
  every field is a genuine Python-vs-C check. The suite also now runs by default
  whenever the kernel is built (was opt-in via an env var).
• A CTest setup fixture (build_all) rebuilds all targets before any test
  runs, as a FIXTURES_REQUIRED prerequisite of every test — ctest can no
  longer run a test binary that predates its sources, and a broken build marks
  dependents "not run" instead of green.

Evidence. Both proven by negative control: injecting a Python-only divergence
made the fixed parity suite fail on exactly the previously-blind fields (F, Pc
diverging 1.3e-2); breaking the C test and running bare ctest (no manual
rebuild) now fails where it previously passed on the stale binary.

---

Reviewer notes / what is NOT done

• The wall-temperature loop is the largest remaining error and is untouched.
  The recession path still uses surface_temperature_limit (a failure
  criterion, hardcoded 1200 K) as the operating wall temperature. The trustworthy
  §3 solver returns ~2282 K for the canonical chamber but (a) isn't wired into the
  recession path, (b) models no ablation energy sink, and (c) its re-radiation-to-
  cold-ambient term is questionable for an inner surface that radiates to hot
  gas. q_conv swings 5.5× between the two wall temperatures — almost
  certainly the dominant contributor to the "3–5× overestimate." Closing it needs
  a coupled surface energy balance (convection + gas-rad + re-radiation +
  conduction + ablation), and two physics decisions (the radiation sink; char vs
  virgin / two-zone). Deferred deliberately, not overlooked.
• Two Dittus-Boelter sites remain (compute_regen_heat_transfer, chamber-
  solver h_hot_base). Neither drives a shipped scalar today (regen disabled in
  all configs; chamber h feeds film cooling + display only). Should move to
  bartz_chamber_h_g for consistency.
• The gas emissivity (0.10) is an engineering estimate, labelled as such in
  code — within Göbel et al.'s band but not derived for our specific
  propellant/Pc/MR.
• Deferred with decisions recorded: bulk cp at chamber_solver.py:972/1140
  (still γ-derived, but this is a bulk energy balance wanting equilibrium cp —
  opposite direction to the transport-chain fix); Alberti/HITEMP emissivity
  correlation + CEA mole-fraction caching; Huzel's R_d carbon-deposit model.
• Verification at tip: 365 passed / 75 skipped, ruff clean, A/B parity 8/8
  against a freshly-built kernel, ctest 8/8.