The simulator

Point-scatterer and Fourier-wavefield backends, scenes, heightfields, and objects.

A GPU-accelerated synthetic aperture sonar data simulator. Produces POSSM-compatible HDF5 files that the beamform_tdbp pipeline consumes directly.

What You Get Out

run_simulation(...) writes a single HDF5 at output_path with the layout expected by beamformer/sas/data_loader.py:

/ReceiverInformation (attrs: SampleRateHZ, CenterFrequencyHZ) and ElementPositionsM/{x,y,z}
/TransmitterInformation/Transmitter1/ElementPositionsM/{x,y,z}
/SimEnvironment attr waterSpeedMetersPerSecond
/Pings/<N> (1-based) with Data (complex baseband, shape (n_channels, n_samples)), sensorPosition, timeSeconds, and {roll,pitch,heading}Degrees

The baseband data is directly usable by the beamformer; no post-hoc conversion required.

Two Backends

run_simulation(scene, motion, cfg, out, method='ps' | 'fw', ...) dispatches to one of:

Point-Scatter (`method='ps'`, default)

Per-ping CUDA pipeline:

Upload scatterer positions + reflectivities + heightfield + normals to GPU once (invariant across pings).
For each ping, compute TX/RX positions in NED, optionally include intra-ping platform velocity and body-frame angular velocity for per-channel motion during TOF.
gpu_raytrace produces delays, amplitudes, visibility (ray-heightfield LOS shadows), and grazing angle per scatterer per channel.
gpu_echo accumulates complex echoes into a (n_channels, n_samples) time series via an atomicAdd CUDA kernel convolved with the LFM chirp.

Accurate but slow: roughly a minute per ping at 1000 scatterers/m^2. Requires CuPy. Supports shadows. shadow_enabled=False disables the LOS check.

Fourier Wavefield (`method='fw'`)

Implements Sanford et al. 2024 (references/sanford2024.pdf) in the paper's four steps:

A. Geometric rendering: the seafloor heightfield is rendered once from nadir (world grid, Sanford Eq. 11); each aspect (spacing fw_aspect_spacing_deg, default 1 deg since 2026-09-10; the regression ladder pins 0.5) rotates the normals (Eq. 4) and runs a per-aspect line-of-sight occlusion pass sheared along y tan(phi) with the true per-pixel grazing angle. Placed CAD meshes are rendered per aspect with a tilted orthographic camera along the sonar ray, first hit per pixel with the facet normal (simulator/fw_mesh.py, FW_MESH_FACETS), so facing walls, self-occlusion and layover come from the geometry; the mesh's surface samples are then not injected as points. Procedural boxes, cylinders, wedges and truncated cones register their own closed mesh and take the same path. Whatever is still injected (cables, manual points, porous meshes) is occluded per aspect against the render's horizon (FW_INJECT_OCCLUDE).
B. Scattering intensity: Lambert from the normal maps with the same BSS calibration as PS, projected onto the slant-range plane; mesh facets add (TS / A_mesh) cos(angle) per camera pixel.
C. Wavefield: random phase screen, 2-D FFT, Stolt mapping to (k_q, ω), system window, inverse FFT. The result is ONE monostatic, stop-and-hop, straight-track wavefield e(q, t). The Stolt gather interpolates along k_r (Lanczos-5), so the range FFT is padded to keep the far edge of the swath inside the kernel's passband (_compute_fw_pad_dims, docs/fw_calibration_tests.md fix 10); a 100 m swath at 18.75 mm needs 16384 range bins.
D. Data generation: per ping and per receiver: gather e at the phase-centre along-track coordinate (FFT upsample fw_upsample_factor, default 8), then apply one exact delay map dτ(ping, ch, t) holding the bistatic excess, sway/heave with range-dependent grazing, attitude and (unless simulation.stop_and_hop) intra-pulse drift (fw_engine.compute_step_d_delays). There is no post-processing after the HDF5 is written; the stop_and_hop attribute is the config value, as for PS.

Not modelled: the angle dependence of the bistatic excess and the off-broadside (Doppler) part of drift / wide-beam sway; both are applied at their broadside value. Validation ladder and tolerances: docs/fw_calibration_tests.md. Design: docs/superpowers/specs/2026-09-09-fw-step-d-exact-design.md.

fw_pixel_size <= lambda/4 (~1.25 mm at 300 kHz).

Both backends write identical HDF5 layout.

Coordinates and the Altitude Rule

NED (x=North / along-track, y=East / cross-track, z=Down). The heightfield stores positive-up elevations; when converted to NED scatterer positions, z = -elevation. Platform altitude above the seafloor is -pos_z; P_0 = [0, 0, -altitude[0]].

run_simulation(..., max_range=...) is a slant range in meters. By convention, SAS platform altitude is roughly 10% of max ground range. E.g. ground range 60 m → altitude 6–8 m → slant range ~60 m.

Design-Rule Warnings and Preflights (2026-09-05)

simulator/scenes/run.py::instantiate runs three checks before any heavy work, on every path (CLI, cloud worker, studio, GUI), not only in the GUI validation layer:

Design rules (simulator/design_rules.py, the formulas of docs/reference/sas_sonar.md): R2 PRF vs range ambiguity, R6 altitude band, R7 along-track advance per ping ≤ L/2, R8 channel spacing vs λ/2, R17 TX azimuth width ≥ channel spacing (ghost rule), R18 baseband Nyquist, R19 elevation null. Printed as [sim-prep] design rule ... lines; never blocks a run. APERTURE_DESIGN_RULE_WARNINGS=0 silences them.
Host-RAM preflight (mem_preflight.preflight_host_ram): raises MemoryError with the density and a fix when the predicted build peak exceeds MemAvailable. APERTURE_SKIP_RAM_PREFLIGHT=1 bypasses.
GPU preflight (mem_preflight.preflight_gpu_ram, PS backend, just before the first device upload): the resident elevation + normals grids (32 B/cell) plus the echo accumulators are checked against the device total (raise) and what is free right now (warn). Same bypass flag.
Object placement contract (run.validate_objects, also run by the studio's validate_scene): every key on an object dict must be a parameter of that kind's placer or one of OBJECT_METADATA_KEYS (kind, scour, annotation, name, id, _scatterer_range); PLACERS in run.py is the single kind -> placer -> renames table. A bad key is a one-line objects[i] (kind 'name'): unexpected key ... error at prep, never a TypeError inside placement (2026-09-05: the studio's name identity key killed a run that way). Guard tests: simulator/tests/test_object_placement_contract.py tie the studio's documented object keys and the GUI's dropped keys to the signatures.

Related loader rule: a bottom.zone_grid_path that does not exist beside the YAML is a FileNotFoundError; the inline zone_grid next to it is only a 1x1 stub, and the old silent fallback to it built zone-less scenes. APERTURE_ALLOW_MISSING_ZONE_GRID=1 restores the fallback (with a warning) for deliberate use.

Scatterer Density

Scene-level parameter (Scene(scatterer_density=...), units scatterers/m^2). Controls speckle fidelity and GPU memory:

Density	Use
`200`	Debug. Debug at 200 first.
`2000`	Standard production runs.
`5000`	High-fidelity speckle / publication figures.
`33000`	Physics-accurate but memory-intensive.

Memory scales linearly: a 3000 m^2 flat_seafloor at 2000/m^2 yields ~6 M scatterers (~3 GB at 32 channels). Iterate at 200/m^2; only bump density when you know the setup is right.

Workflow

from simulator import Scene, PlatformMotion, SonarConfig, run_simulation

cfg = SonarConfig()                               # 300 kHz / 75 kS/s / 60 kHz BW defaults
motion = PlatformMotion.straight_line(
    speed=1.5, altitude=8.0, n_pings=128, ping_rate=3.3,
    sway_amplitude=0.15, sway_period=10.0,        # optional 6-DOF perturbations
    roll_amplitude=0.01,
)

scene = Scene(
    x_min=-5, x_max=65, y_min=0, y_max=70,
    cell_size=0.05, bottom_type='sand',
    scatterer_density=200.0,                       # debug first!
)
scene.set_bottom_zone_map(                         # optional zoned seafloor
    [['sand'], ['rock']],
    rms_heights={'sand': 0.001, 'rock': 0.1},
    seed=42,
)
scene.add_procedural_rocks(                        # optional rock population
    zone=(1, 0), density=2.0, d_min=0.2, d_max=1.0, seed=543,
)
scene.add_box_object(                              # objects raise the heightfield
    x_center=20, y_center=20,
    x_extent=2.0, y_extent=1.0, height=0.5, ts_db=-15.0, seed=43,
)
scene.build(seed=42)                               # must be called before run_simulation

run_simulation(scene, motion, cfg,
               output_path='sim_output/scene.h5',
               method='ps', max_range=60.0,
               shadow_enabled=True, verbose=True)

The Heightfield

2D elevation grid (n_x, n_y) with cell_size. Owns:

elevation[n_x, n_y]: positive-up heights (m). Objects raise local cells to max(existing, object_top).
bottom_type (uniform) or bottom_type_grid (per-cell index into _zone_bottom_types) when SeafloorBuilder runs a zone map.
object_scatterers: explicit list of (x, y, z_ned, complex_refl) tuples added by place_*.
Surface normals via compute_normals() with a 60 deg max-slope clamp (prevents single-cell object edges from producing near-horizontal normals that silence seafloor scatterers).

Primitives on the heightfield (all exposed via Scene.add_*):

place_box_object(x_center, y_center, x_extent, y_extent, height, ts_db)
place_cylinder_object(x_center, y_center, length, diameter, orientation='along_track'|'cross_track', ts_db): horizontal cylinder
place_rockan_mine(x_center, y_center, heading, ts_db): tapered wedge (~1.0 x 0.5–0.8 x 0.4 m)
place_manta_mine(x_center, y_center, ts_db): truncated cone (~0.98 m base x 0.44 m tall)
place_mesh_object(x_center, y_center, mesh_path, scale, rotation_z, ts_db) : loads STL/OBJ/PLY via trimesh, rasterizes into heightfield via downward ray shooting, then samples surface scatterers uniformly
place_superellipsoid_rock(a, b, c, n_shape, ts_db, ...) ; |x/a|^n + |y/b|^n + |z/c|^n = 1. n_shape<2 pinched, >2 boxy
place_procedural_rocks(x_min, x_max, y_min, y_max, density, d_min, d_max, alpha) : Poisson-placed superellipsoids, log-normal axis ratios, size-dependent n_shape (large rocks angular, small rocks rounded)
add_point_scatterer(x, y, z, ts_db, n_sub=10, jitter_m=0.01) ; Rayleigh sub-scatterers inside a small jitter disc (avoids the single-pixel coherent halo)

The SeafloorBuilder

Driven by Scene.set_bottom_zone_map(grid, rms_heights, seed) and Scene.add_procedural_rocks(...). During scene.build():

Divide the heightfield XY extent into equal-area rectangular zones, one per cell of grid (rows index along-track).
For each zone, synthesize FFT power-law roughness S(k) ∝ k^(-β) using the bottom type's spectral_exponent, normalized to rms_height_m (or the rms_heights override). A cosine taper at zone edges avoids boundary discontinuities.
Cells already occupied by placed objects are protected via protect_mask so object tops stay smooth.
Rocks auto-populate by bottom_type.default_rock_density unless disabled; explicit add_procedural_rocks overrides apply last.

Optional Tang ripple overlays (scene.set_tang_ripples(ripple_wavelength, rms_height, direction_deg, skewness_b)) layer on top: skewness_b>0 uses the Tang eq. (5) exponential transform for sharp peaks / wide troughs (b=18 m^-1 matches Traykovski data).

Bottom Types

Defined in simulator/bottom.py. BSS at 20 deg grazing calibrates the amplitude; spectral params drive FFT macro-roughness.

name	mean_reflectivity	roughness_std (m)	bss_db_at_20deg	β (spectral_exponent)	spectral_strength	rms_height_m	default_rock_density
`sand`	0.30	0.002	-28.0	3.25	1.41e-4	0.02	0.0
`mud`	0.15	0.001	-32.0	3.0	5.0e-5	0.005	0.0
`rock`	0.60	0.010	-15.0	2.0	6.0e-4	0.15	15.0
`gravel`	0.45	0.005	-22.0	2.5	3.0e-4	0.08	3.0
`silt`	0.20	0.001	-38.0	3.0	3.0e-5	0.003	0.0

roughness_std is sub-wavelength per-scatterer z-jitter; rms_height_m drives FFT macro-roughness from SeafloorBuilder. spectral_strength is currently reserved: the FFT path normalizes output to rms_height_m regardless. bss_db_at_20deg=None means uncalibrated; Scene will refuse to build unless you pass allow_uncalibrated_bottom=True.

Platform Motion

PlatformMotion.straight_line(speed, altitude, n_pings, ping_rate, heading=0, ...) generates a straight track. Sinusoidal perturbation pairs (sway_amplitude+sway_period, heave_*, roll_*, pitch_*, yaw_*) inject 6-DOF wiggle, all zero by default. Angles are radians. motion.altitude returns -pos_z; motion.time is arange(n_pings)/ping_rate.

motion.kind: spline_track builds a curved track instead, through simulator/track.py: a centripetal Catmull-Rom curve through track.control_points, a per-point speed law (piecewise-linear in arc length), and a kinematic attitude derivation with a first-order vehicle response per axis (track.dynamics). The FW backend requires a straight aperture, so its control points must be colinear at one altitude (rule R12d, appcore/validation/rules.py); the PS backend has no such restriction.

SonarConfig Defaults

Match POSSM scene066: fc=300 kHz, fs=75 kHz, BW=60 kHz, pulse_length=4 ms, sos=1500, 36 channels at 33.33 mm spacing, 3 cm square TX/RX elements. create_baseband_chirp() returns the LFM waveform used for matched filtering; beamwidth is computed on demand via tx_beamwidth_rad / rx_beamwidth_rad (never hardcode).

Point-Scatter Performance (GPU path)

Beam-wedge tiling of a streamed cloud (2026-09-13)

When the cloud is too large to stay GPU-resident, the prefilter used to stream the whole cloud through PCIe on every ping (at 80k/m^2 on a 71 x 100 m scene: 580 M scatterers, 11.6 GB per ping, about 3 s of each 4.3 s ping, GPU at 30 to 40 %). simulator/cloud_window.py now sorts the cloud once into 2 m (x, y) tiles (stable radix sort on int16 keys, RAM-gated) and, per ping, streams and scans only the tiles the TX beam wedge can reach: a conservative per-tile test derived from the prefilter's own keep rule, with the tile's circumscribed radius and the cloud's z extent folded in, so the survivor set is unchanged (property tests over random attitudes, both sonar sides and elevated scatterers; byte-identical HDF5 and DRC output on a straight scene and on the arcing testing_fw_linear_track scene at 5000/m^2). Because the wedge is a thin fan around broadside, the selected share is similar at every heading: 19 % of the cloud per ping on the straight scene (6.3 to 15.4 pings/s, streaming forced), 33 % on the 90-degree arc (3.75 to 6.13 pings/s). The log reports the tiling time, ping 1's selection and the run's average streamed share. PS_X_WINDOW=0 disables it; the sort is skipped with a log line when free RAM cannot cover its transient (22 bytes per scatterer).

Resident clouds tile on the GPU (2026-09-13, later that day). A cloud that will stay device-resident (<= 25 % of free VRAM) is uploaded as is and sorted there by cloud_window.tile_cloud_gpu: same tile grid, same stable permutation as the host path (pinned by test_gpu_tiling_reproduces_the_host_permutation_exactly), so the HDF5 is unchanged, and there is no second host copy of the cloud. The host copy had pushed the build-to-sim peak RSS 0.9 GB over the preflight budget at 96 M scatterers (test_build_and_sim_handoff_peak_rss_budget). Streamed clouds still tile on the host, page-locked, because the sorted copy is what the per-ping copies read from.

Parallel tile sort (2026-09-13, evening). The host tiling is a counting sort run block-wise on the shared thread pool (simulator/cpu_pool.py, PS_BUILD_THREADS, 1 = serial): each block of 4 M rows computes its tile keys and a per-tile count, an exclusive prefix sum over blocks gives every block the destination of its rows for every tile (tile start + that tile's rows in earlier blocks), and each block stable-sorts its own keys and scatters its rows straight into the output. Blocks are contiguous index ranges, so the result is the same permutation as the global np.argsort(keys, kind="stable") it replaces (test_parallel_tile_sort_matches_the_global_stable_sort, serial and threaded, odd block sizes). The page-locked output is allocated on the pool first so the locking overlaps the key passes. 100 M rows with page-locked output (WSL): 21.6 s serial -> 5.0 s threaded, i.e. about 104 s -> 30 s for the 580 M cloud.

Single-pass, double-buffered streaming (2026-09-13, same day)

The streamed prefilter used to count survivors in one pass and fill them in a second, so every selected chunk crossed PCIe twice per ping, and each copy was synchronous. It is now one pass: two staging slots on the GPU, host-to-device copies issued on a dedicated non-blocking stream so chunk k+1 lands while chunk k's keep mask and compaction run, and survivors appended to persistent output buffers that grow on demand (no per-ping allocation, no concatenate). The tiled cloud is allocated page-locked when it is at most 40 % of free RAM (PS_PIN_CLOUD=0 keeps it pageable), so .set() streams straight from it with no host memcpy; a pageable cloud still bounces through two chunk-sized pinned buffers. Output order is chunk order, identical to before: the straight scene's HDF5 and DRC are bit-identical to the original two-pass unwindowed run. Streaming forced at 5000/m^2: straight 15.4 to 25.6 pings/s, the 90-degree arc 6.1 to 13.0 pings/s. PS_TIMING=1 prints the prefilter and whole-ping milliseconds with every progress line; on the arc scene the prefilter is now about 40 ms of an 85 ms ping.

Reference: the File -> New Scene default (60 m x 80 m image, HISAS-class 36-channel array, 150 pings) at 10,000 scatterers/m^2 = 57.6M scatterers, RTX 4500 Ada. Sim stage 226 s -> 19 s (per ping 1.56 s -> 0.11 s) on 2026-09-02; the beamformed image correlates 0.9998 with the old path (differences are fp32 summation-order speckle, ~53 dB down).

Per-ping stages, in the order they run, with the switch that restores the previous implementation. All are ON by default.

Stage	What it does now	Off switch
Cloud residency	The scatterer cloud is uploaded once and stays on the GPU when it needs <= 25 % of free VRAM; larger clouds stream from host chunk-wise as before (`simulate._resident_or_pinned_cloud`).	`PS_RESIDENT_CLOUD=0`
TX-beam prefilter	One mask kernel + one compaction over the resident cloud (`gpu_engine._prefilter_resident`) instead of a two-pass CuPy chain.	`PS_PREFILTER_KERNEL=0`
Line of sight	`check_los_hard_mip_kernel`: the ray march skips any segment (<= 8 cells long) whose highest point clears the 3x3-dilated max elevation of its coarse cell; the fine samples that do get tested are exactly the legacy ones, so visibility is bit-identical.	`PS_LOS_MIP=0`
Raytrace	Fused CUDA kernel (`fused_raytrace_kernel`, one thread per channel x scatterer) writes delays + complex64 amplitudes directly; replaces the CuPy elementwise chain (~10 (n_ch, n_sc, 3) temporaries per chunk).	`PS_FUSED_RAYTRACE=0`
Cull + sort	`raytrace_ping_fused_culled`: a per-scatterer score kernel yields max-channel amplitude and mean delay; the cull mask and delay argsort run on those two vectors; survivors' inputs are compacted and only then does the fused kernel write the (n_ch, n_keep) arrays, already sorted. Same keep rule as `_compact_active_scatterers`.	`PS_FUSED_CULL=0`
Echo (FFT path)	`accumulate_deltas_c64_kernel` reads the complex64 amplitudes directly and uses `sincospi` on the fractional carrier cycles (186 dB vs the split real/imag + `cos`/`sin` kernel).	`PS_DELTA_C64=0`
Raytrace + echo in one pass	`raytrace_echo_kernel`: one thread per delay-sorted survivor (indexed through the cull's index list, no compaction), looping over channels, scattering each pair's delta straight from registers into the FFT delta planes. Runs of equal bin address within a warp are summed by a segmented shuffle scan before ONE fp64 atomicAdd pair. The carrier phase is float32: fc·τ split with an FMA two-product so the fractional cycle is exact to ~6e-8 (fp64 `sincospi` on the raw ~3e4-cycle argument ran at 1/64 rate and was 3/4 of the kernel). Chunks are sized for this path's ~64 B/scatterer working set, so the default scene is one chunk per ping and the cull runs against the true per-ping peak. 146 dB vs the two-kernel path at equal chunking.	`PS_FUSED_ECHO=0`

None of this depends on the amplitude cull being enabled: with echo_amplitude_cull = 0 (Preferences -> Run) the fused pass keeps every visible scatterer and skips the score kernel. Only the TX-beam prefilter stays tied to the cull, since it is itself a cull. On the 4125i scene at 20k/m^2 (145M scatterers, 103 pings) the sim stage is 90 s at cull 0, 44 s at 1e-4, 36 s at 1e-3; before the decoupling cull 0 took 5 min. The default is 1e-4 since 2026-09-03 (simulate.run_simulation, the GUI preference run/echo_amplitude_cull, and therefore browser/AWS runs): 1e-3 culled the dim far-range ripple-trough and shadow-edge scatterers, so ripple fields beyond ~75 m broke into along-track dashes and point targets grew sidelobe crosses; 1e-4 is visually identical to culling off.

The ping loop also trims CuPy's memory pool whenever its cached total passes 60 % of the free VRAM measured at sim start (PS_POOL_TRIM_FRACTION). The fused path's per-ping transients change size with the survivor count, so cached blocks stop matching and the pool only grows; on the 145M-scatterer 4125i scene it reached 46 GB "total" on a 24 GB card. WSL's driver spills that to host memory rather than failing, and mid-track pings ran 2-10x slower (worse from the GUI, which holds VRAM of its own). With the trim: 67 s -> 36 s for that scene.

Second round, same scene at 20,000/m^2 (115M scatterers, 2.3 GB resident): per ping 215 ms -> 98 ms, sim 34.6 s -> 18.7 s, wall 71 s -> 56 s. What remains per ping (~98 ms): echo/raytrace kernel 30 ms, LOS 21 ms, score 14 ms, prefilter 12 ms, argsort 6 ms, the rest in compaction/finalize. LOS coarse factor is 16 (measured fastest, bit-identical at every factor). Outside the sim, instantiate() is ~10 s (the single-zone scatterer draw now writes through slices instead of a 57.6M-element fancy index; the zone-feather gaussian is 2.3 s of it) and the two PNG products ~7 s (zlib level 1; speckle does not compress and the row filter dominates).

Profiling recipe: simulator/perf/bench_echo.py (CUDA-event timers on a synthetic scene) or wrap gpu_engine / gpu_raytrace entry points with event timers around generate_ping_data_gpu on a real YAML.

Object placement hands arrays to the store (2026-09-13)

Every place_* method used to feed ObjectScattererStore.append one (x, y, z_ned, refl) tuple per scatterer. On the 9obj scene at 80k/m^2 that loop ran 39 M times and was 77 s of a 641 s run (the third-largest stage after echo simulation and the tiling sort). The placers now call ObjectScattererStore.extend_arrays(xs, ys, z_ned, refl) with the arrays they already had (z_ned may be a scalar); the store keeps the same float64 / complex128 chunks and the same order, so the object cloud is bit-identical to the tuple path (checked for every placer, including the Manta top disc whose per-point cos/sin became one vectorised call). The legacy append face stays for the texgraph placement drain and tests. Parsed mesh geometry is also cached per (path, mtime, size) as read-only arrays, so a scene that places the same OBJ many times parses it once. Measured on the 9obj scene (WSL): placement 77 s -> 20 s; what remains is trimesh's surface sampler and OBJ parsing of distinct files.

Mesh objects prepared in parallel (2026-09-13)

place_mesh_object is now two steps. prepare_mesh_object is pure: it loads the mesh, transforms it, rasterises its top surface over the footprint window and draws the surface scatterers from the object's own seed, reading only grid geometry and density. apply_mesh_prep writes the heightfield stamp (or feather), the protect and footprint masks, the facet record and the scatterers. instantiate() prepares every mesh object at once with prepare_mesh_objects on the shared thread pool (PS_BUILD_THREADS, 1 = inline), then applies them in YAML order between the other kinds, so overlapping stamps (max against the current bed) and per-object scatterer ranges are exactly what the one-by-one placer produced (test_prepared_mesh_placement_matches_ sequential_placement; the full 9obj placement pass was also compared against the committed code: scatterers, elevation, footprint, facet records and ranges identical). Distinct mesh files are parsed once, serially, before the threads start: trimesh's OBJ parse holds the GIL, and letting 16 threads contend for it made a low-density scene slower than serial. 9obj scene (WSL): placement 18.7 s serial -> 11.1 s at 80k/m^2, 4.3 -> 3.7 s at 500/m^2. What remains is mostly that serial parse (about 9 s for the nine distinct OBJ files); a process pool for the parse alone is the obvious next step if it matters.

Texture-graph nodes evaluated in parallel (2026-09-13)

evaluate_graph used to walk the topological order one node at a time. Nodes are pure functions of the context, their connected inputs and a seed derived from (base_seed, node_id), so independent nodes now run concurrently: a node is submitted to the shared thread pool (PS_BUILD_THREADS, 1 = the serial loop) as soon as every node it links from has finished, results are stored on the calling thread, and a producer is released only once every consumer has finished. The material palette keeps its last-writer-in-topological-order rule. Same node code on the same inputs, so the height, reflectivity, gains, material ids, palette and masks are byte-identical to the serial order (simulator/tests/texgraph/test_parallel_eval.py, plus the real 9obj graph applied to its heightfield compared against the committed evaluator). The 9obj graph is three generators (tang_ripple 4.6 s, musgrave 2.7 s, voronoi 2.6 s in WSL) and cheap operators, so the parallel time is the slowest generator: 10.5 s -> 7.1 s (WSL). The ripple node is a serial 10-iteration Hilbert/FFT loop; multi-threaded FFTs inside it are the next step if that stage matters.

Threaded scatterer sampling (2026-09-13)

Heightfield.generate_scatterers draws one PCG64 stream in a fixed order (uniform x, uniform y, then per zone normal / rayleigh / uniform phase), and cached references pin the cloud bit for bit to that stream, so the draws cannot be parallelised. They are only about a third of the sampler though: at 580 M scatterers roughly 30 s of draws against 55 s of cos/sin, cell indexing, grid gathers and float32 casts. Those now run on a bounded thread pool (_BlockRunner, PS_BUILD_THREADS, default min(16, cores), 1 = the serial reference path) while the calling thread keeps drawing the next block; numpy releases the GIL for all of it. Every draw still happens on the calling thread in the same order with the same block size, and each block gets the same operations it always did, so the output is byte-identical to the serial and to the historical unblocked path (test_threaded_build_is_bit_identical compares all three on a multi-zone scene with both texture gains and on a single-zone scene). The x and y draws are now blocked too, which removes the last full-length float64 transient (4.6 GB at 580 M). Measured on the 9obj scene at 40k/m^2 (290 M, WSL): scene.build 41.0 s serial -> 21.3 s threaded.

Calibration Harness

calibrate_simulator.py validates simulator fidelity against measurement targets. It builds three scenarios: level (flat seafloor, checks speckle CV), PSF (single calibration target, checks point-spread function), and ghosts (targets at multiple ranges, checks ghost suppression): runs the simulator + beamformer, and writes sim_output/calibration/calibration_report.txt with CV, dynamic range, contrast, and sidelobe metrics. Design spec: docs/superpowers/specs/2026-04-15-simulator-calibration-design.md.

Quick Presets

simulator/presets.py ships ready-to-run (scene, motion, cfg) tuples: flat_seafloor(), point_targets(), calibration_target(range_m=50), rippled_seafloor(), complex_scene(). All call scene.build() before returning.

← Signal processing Scene framework →