HELIOS

GPU shape optimisation for horn-loaded loudspeakers. Stellarator optimisation architecture, acoustic physics kernel.

v1 scope: one mid-range horn. Throat to mouth, rigid walls, single device. No cabinet, no crossover, no array.

Picking this up? Read HANDOFF.md first. It has the traps, the three open problems with everything measured about them, and the ordered list of what to do next.

Status

Milestone	State
M0 — Webster TMM screen	passing
M1(a) box eigenfrequencies	passing
M1(b) duct cross-mode	passing — 0.09 % vs the 2 % asked
M1(c) CFS-PML	passing — −85 dB normal, −75 dB grazing vs −55 dB asked
M1(d) fp32 energy drift	passing — 2.8×10⁻⁶ % vs the 0.1 % asked
M1 final, FDTD vs TMM	partial — meets 1.5 dB above 1.5 kHz, not below. See below
M2 geometry / voxeliser / CAD	passing, with one shortfall recorded below
M3 objective + report	passing
M4 CMA-ES loop	first real run failed honestly — the optimiser exploited the grid; caught by held-out regressions and grid convergence. See HANDOFF §3.4
M5 physical validation	cannot be done in software — see below
M6 adjoint / topology	not started; explicitly gated on M5

103 tests, all green. pytest -m "not slow" runs the fast 83 in ~45 s.

The three results that are not simply "passing"

M1 final: the Tier 1 mouth model, not the solver

The FDTD and the Webster TMM agree within 1.5 dB above about 1.5 kHz (2.5 dB worst case, 1.2 dB rms) and disagree by up to 10 dB below it.

The gap is Tier 1's, and the evidence is that everything checkable in closed form checks out:

• an anechoic duct returns Z = ρc/S to within 0.02 dB at every grid spacing tried, which exercises the whole measurement chain;
• the horn answer is unchanged to 0.00 dB across a 9× domain-volume range, so the PML is not truncating the mouth near field;
• unchanged to 0.01 dB from a 60 ms record out to 400 ms, so nothing is truncated mid-ring.

Below 1.5 kHz this horn is sharply resonant, and a lumped flat-piston mouth termination cannot place its resonances better than a few percent in frequency — which at a steep |Z| peak is many dB. Consequence: Tier 1 stays sound for its actual job in §5 — rejecting on cutoff and gross misloading — and is not trusted for fine low-frequency ripple. Sharpening the mouth model is the obvious next improvement.

M4: it passed, and that is the problem

A 14-generation, 168-evaluation run (13 min) improved J from 495 to 305 — 38 % over the analytic exponential seed, with no held-out metric degraded. That is M4's acceptance criterion, met literally. The winning shape is nonetheless worse as a horn.

Per-term, the whole gain came from one place:

term	seed	best	Δ
group_delay_std	271.3	66.9	−204.5
throat_resistance_min	36.1	22.6	−13.6
directivity_error	101.1	98.5	−2.6
p0_mean	36.9	53.5	+16.5
throat_resistance_ripple	29.6	35.3	+5.7
mfg_draft_angle	3.8	8.9	+5.1
mfg_mouth_area	0.0	2.1	+2.1

group_delay_std was 55 % of J at the seed, because its σ of 30 µs is several times tighter than any horn achieves. The optimiser did the rational thing: it bought a halving of group delay with plane-wave purity — the dominant physics term, the quasi-symmetry analogue this whole project is built on — plus throat ripple, polar smoothness, mouth area, and a draft angle that collapsed from 1.04° to 0.011°, which is unmouldable.

Two conclusions.

The §4.5 held-out mechanism has a blind spot. Held-out metrics catch the optimiser exploiting something nobody measured. They cannot catch it trading measured terms against each other in a way the weights permit but the designer never intended — every term involved is scored, so the held-out set never moves and reports "no regressions" while the design degrades. scored_regression() and dominant_term() in optimize/cmaes.py close that gap; the second flags any weight set where one term exceeds 40 % of J, which is worth checking before committing eight hours to it.

The weights have since been retuned from measurement and probe-validated (see HANDOFF §3.3 for the full story): group_delay_std's σ is now 200 µs — the seed measures 267 µs, so the spec's 30 µs was unreachable — and the draft-angle, mouth-area and stored-energy caps are steep walls. Two 168-evaluation probes under the new weights show the failure mode is gone: the first bought its gains correctly (P₀ 0.903 → 0.936, draft rising), and with every wall up the second could no longer beat the analytic seed at probe scale at all — which is honest. Beating it now takes a real overnight run, and that is the next step for M4.

M2: sub-cell weighting is 8.6 dB, not 10 dB

Sub-cell volume weighting drops the spurious higher-order-mode floor by 8.6 dB against the 10 dB §7 asks for. It works and works substantially — the floor sits near −26 dB, well below the spread between candidate shapes — but it does not clear the bar. A linear smeared Heaviside is a first-order estimate of a cell's fluid fraction and is biased for curved walls; exact polyhedral clipping is the next step, and amounts to the "cut-cell / conformal treatment" the spec prefers anyway.

M5 cannot be done here

M5 needs a printed 1:4 horn, a compression driver, a turntable and a measurement mic. No software can produce that measurement. What exists is the harness on both sides of it — validate/measured.py reads REW/ARTA exports, applies the ka scale mapping, gates the impulse, and applies the 2 dB / 3 dB tolerances — tested against synthesised data. Until a real horn is printed and measured, M5 is unproven, and §11.4 applies: nothing downstream is trustworthy on the strength of the solver alone.

M6 is explicitly gated on M5 in the spec, so it has not been started.

Setup

python -m venv .venv && .venv/Scripts/activate && pip install -e ".[dev,gpu,opt,cad]"

Requires CUDA 12.8+ for Blackwell (sm_120). Verified on Warp 1.16 / CUDA 12.9.

pytest -m "not slow"

python -m helios.physics.webster --flare exponential --cutoff 300 --mouth-radius 0.23

Measured on this workstation (RTX 5060 Ti 8 GB)

Quantity	Spec estimate	Measured
Reference grid	18.0 M cells	16.9 M (336×224×224, dx 3.57 mm)
Raw kernel throughput	450–670 GB/s	312 GB/s (6.24 G cell-updates/s)
Wall time per candidate	10–20 s	~10.5 s (4.35 s simulating, 5.95 s voxelising)
VRAM per candidate	0.6 GB	0.44 GiB → 12 concurrent
fp32 energy drift / 8 000 steps	< 0.1 % required	2.8×10⁻⁶ %

Two deviations from the spec's hardware assumption, neither blocking: the card is the 8 GB 5060 Ti, and throughput is ~70 % of its 448 GB/s peak.

The probe sampling now happens on the device (sampling="device", the default on Warp). The loop used to download the whole pressure field every step to read a few thousand fixed points — on the M4 candidate, 11.6 MB moved to extract 19 KB, which cost 3× the kernel. Gathering those points in a Warp kernel took the step loop from 3.02 to 0.679 ms per step (4.45×), against a kernel-only floor of 0.466 ms. Every recorded array is bit-identical either way; sampling="both" reads both paths from the same field in one run, which is how test_device_and_host_probe_sampling_are_identical asserts exact equality.

The §1 free speedups are otherwise taken: quarter symmetry is on by default, and early termination at −60 dB is implemented — though note the M4 winner never triggers it and runs the full 6 400 steps.

What is now the bottleneck: voxelisation, on the CPU. One candidate splits 5.95 s of setup against 4.35 s of simulation and 0.2 s of post-processing, and 5.85 s of that setup is voxelize_horn — almost all of it bernstein_basis and bernstein_derivative_basis re-evaluated over every voxel coordinate on each of four SDF calls. That, not the GPU, is the next thing worth attacking.

Layout

Follows spec §10, under src/ so tests import the installed package.

src/helios/
  physics/    webster.py  fdtd_kernel.py  boundaries.py  farfield.py  simulate.py
  geometry/   param.py  voxelize.py  export.py
  objective/  terms.py  assemble.py  presets/{balanced,f1_character,constant_directivity}.yaml
  optimize/   cmaes.py  checkpoint.py
  validate/   measured.py
  report/

simulate.py is the one addition to §10's layout: M1's cross-check, M3's objective and M4's loop all need the same "put a horn in a box and drive the throat" setup, and a leak in it is far easier to find in one place than in three copies — which is exactly what happened.

Bugs worth knowing about

All four presented as something other than what they were, which is the reason each is recorded rather than silently fixed.

1. Sub-cell weighting was a porous membrane, not a cut cell. Face fractions alone throttle flux but leave the cell filling its whole volume, so boundary cells respond too slowly and modes land flat — the duct cross-mode came out at 534 Hz against 1005. The cell volume fraction has to divide the pressure update too.

2. Faces could open into zero-volume cells, draining energy into frozen pressure. Now an invariant inside set_geometry.

3. The horn body's axial truncation smeared into the throat plane, leaving the cell layer carrying the piston a third open across the whole domain. Whether it happened depended on where the cell centre fell relative to z=0, so it came and went with grid spacing and read for a long time as an erratic resolution requirement.

4. plane_wave_purity was computing the wrong integral. P₀ needs ∫|P(x,f)|² dA — the area integral of the squared spectrum at each point — and it was transforming the time-domain ∫p² dA, which is a convolution and unrelated. The symptom was P₀ = 1.000 for every shape at every frequency: the dominant term of the objective measuring nothing, and looking perfect while doing it.

Deviations from the spec, and why

• The reference file the-bench-soundsystem-physics.html was not in the workspace. webster.py is implemented from the §5 equations rather than ported from it. It validates against the exact analytic infinite-exponential horn, but the mouth-termination and wall-loss coefficients are engineering choices rather than derivable — worth diffing if that file turns up.

• SIMSOPT is not imported. What is actually reused is the objective algebra, which is a page of code; the rest of SIMSOPT is an MPI-oriented parallel layer and a fusion dependency tree. The assembly shape is kept the same so a later switch is mechanical.

• The §4.4 time-domain terms use a 700 Hz–8 kHz sub-band. Group delay through a horn diverges at its cutoff — that is intrinsic to the flare, not a bad shape. Measured across the full band a perfectly good 339 Hz horn shows 1.25 ms of spread against the 30 µs asked for, and the term then contributes 98 % of J and drowns out everything else. The sub-band is the range §4.3 already uses for the directivity target.

• CFS-PML defaults to alpha_max = 0. The CFS pole is implemented and parameterised, but made no measurable difference in either the normal or the grazing configuration, so the default reflects the measurement rather than the intent.

• The wall-loss coefficient in §5 is ~18× the physical thermoviscous value. Implemented as specified and exposed as alpha_coeff; the FDTD comparison sets it to zero, since the FDTD carries no wall loss in v1.

• 48 TMM segments is not converged at 8 kHz. Tier 1 is second-order; max error in normalised throat resistance is 0.076 at 50 segments and 0.018 at 200. Fine for screening, which is all it does; constructors default to 200.

Next

1. Move probe reductions into Warp kernels — the candidate loop is host-bound.
2. Exact polyhedral cell clipping, to clear M2's 10 dB bar and sharpen P₀.
3. A better Tier 1 mouth model, to close the sub-1.5 kHz M1 gap.
4. Recalibrate the objective weights before any long run. The demo showed one term at 55 % of J; a run that unbalanced spends the night on it.
5. Print a 1:4 horn and measure it. Everything above is unvalidated until then.