pressure-wave — ray-traced blast overpressure on body armor

pressure-wave is a GPU ray tracer that estimates how much blast overpressure reaches a spherical sensor seated inside a soldier’s helmet, under five armor configurations, for a cube of candidate blast points around the body. It is a small research codebase aimed at quantifying the protective contribution of a MICH-style ballistic helmet and SAPI torso plates against the impulse of a free-air blast, and at producing a tidy CSV that downstream R scripts can run an ANOVA and paired comparisons against.

The numerical core runs on NVIDIA Warp, one GPU thread per ray, each thread marching its own bounce path against a hardware BVH. The five configurations are full_armor (helmet plus front and back plates), front_plate_helmet, helmet_only, vest_only (the front plate alone), and no_armor. Every parameter of a run lives in config.yaml, with repeatable dotted --set key.path=value flags for one-off overrides.

Goals

Estimate, for a given peak pressure, the total impulse delivered to the sensor at the origin from blast points scattered across a configurable cube.
Compare all five armor conditions on the same set of blast points so each point acts as its own paired “subject” for the statistical analysis.
Stay simple enough to read end-to-end: a few hundred lines of Python for the physics, a few dozen lines of R for the stats, and two .obj meshes for the geometry.
Be reproducible. A run is described completely by its YAML file plus any --set overrides, and writes a timestamped CSV next to the others in data/.
Be fast enough that ray resolution stops being a budgeting decision.

How it works

Host-side preparation lives in sampling.py, which is plain float64 NumPy. It builds the blast grid, discarding points closer to the sensor than blast.min_dist_m, derives an orthonormal frame per blast point aimed at the sensor, and lays out the padded azimuth/elevation offsets that describe that point’s ray fan. By default the fan spans only the angular extent of the sensor disc (rays.fan_mode: sensor_cone); the geometry mode widens it to cover the armor so that reflected paths are sampled too. The same module precomputes K, the Friedlander shape integral over the positive phase, since the waveform is linear in peak pressure and a ray’s impulse is just its effective pressure times that constant.

The device core (warp_core.py) takes those compact arrays and generates the rays inside the kernel, so the full ray set is never materialized in memory. Each thread owns one (explosion, ray) pair: it queries the armor meshes with Warp’s BVH wp.mesh_query_ray, tests the sensor sphere analytically over the same segment, and takes whichever hit comes first. A geometry hit reflects specularly, scales the carried energy by physics.reflection_loss (0.8 by default), and continues until the ray reaches the sensor, runs past rays.max_length_m, or exhausts rays.max_bounces. A ray that reaches the sensor contributes the initial peak p₀ attenuated by an inverse-square distance law and by the accumulated reflection factor, multiplied by K, weighted by the ray’s differential solid angle and atomically added into the per-explosion accumulator.

simulate_blasts.py drives the run, launching the kernel once per armor configuration. The configurations differ only in which meshes the kernel may query, so no_armor is the same code path with the mesh count set to zero. The front and back plates are merged into one BVH, which keeps full_armor inside the kernel’s two-mesh limit.

One consequence of the specular model is worth stating plainly. The sensor sits inside the helmet shell, so armor can amplify the impulse as well as attenuate it, with the shell reflecting rays back onto the sensor. No ordering between configurations is guaranteed at a given blast point, which is why the test suite asserts invariants that hold either way instead of a ranking.

Edge diffraction

Pure ray tracing casts sharp geometric shadows: a blast point the armor occludes contributes exactly zero. The pulse is acoustically large, though. Its positive phase spans roughly 3.4 m, about ten times the size of the helmet, so real blast energy diffracts around the shell into those shadows. Setting diffraction.enabled: true adds first-order edge diffraction (diffraction.py) so the geometric shadows fill in physically. It exists mainly as a robustness check: does the ranking of armor configurations survive the model’s single largest omission?

The method is Biot–Tolstoy–Medwin (BTM), chosen because impulse is the zero-frequency (DC) component of pressure. The high-frequency GTD/UTD diffraction coefficient diverges there, scaling as 1/√k, while the BTM time-domain response is integrable and yields a finite DC edge gain. The diffracted impulse is the incident impulse times that gain, folded into the same accumulator as the direct and reflected rays, so geometric optics plus diffraction is continuous across the shadow boundary by construction. Because MICH.obj is a smooth watertight shell with no rim edge, the helmet’s diffracting edge is sourced procedurally as the intersection of the shell with a horizontal plane; the torso plates use their perimeter feature edges.

The feature is off by default, and its caveats are real: first order only, no double diffraction or reflection×diffraction, edge diffraction rather than smooth-surface creeping waves, impulse-only, linear acoustics. The diffraction term is signed, so a transition cell whose near-zero geometric field the ray fan under-resolves can sum slightly negative. Those cells are clamped at zero and counted in the log, and the effect shrinks with finer rays.resolution_arcmin. Diffraction-on results are best read as a comparative bound, well short of absolute dosimetry.

Validation

uv run pytest runs the suite in tests/ against whatever compute.device resolves to, so it works on the Warp CPU backend when no GPU is present. test_physics.py checks the physical invariants and cross-checks the no_armor column against an independent NumPy computation of the same quantity. test_diffraction.py is a validation ladder for the diffraction coefficient: an exact static rigid-wedge Green’s function (a Legendre-Q modal series) is checked against image sums to 1e-9, the BTM coefficient is checked against that oracle in the deep shadow across wedge angles including the screen case, and the shadow-boundary continuity condition independently pins the coupling constant 1/(16π²). That constant is fixed by those checks rather than fitted to the data.

Analysis

The R side (analysis.R) reads a results CSV, splits blast points into near and intermediate-range bands, pivots to long form with each blast point as its own subject, and runs a repeated-measures ANOVA on impulse vs. armor condition. When the omnibus test is significant, it follows up with paired t-tests across every pair of conditions the file actually carries, then applies a Holm correction: six pairs on the older four-column CSVs, ten once front_plate_helmet is present. It also reports the helmet’s effect size as both an absolute impulse reduction (with 95% CI) and a Cohen’s d on the per-subject percent reduction. visualize.R contains the box-, violin-, and 3D scatter plots used to eyeball the same data, and Exploration.ipynb walks the whole pipeline interactively, from configuration and geometry through the ray fan to a small Warp run with plots and a NumPy cross-check.

Inputs and outputs

Geometry. Two .obj meshes in assets/: a MICH-style helmet and a SAPI plate. geometry.py rotates and scales them into the project frame, which is Z-up, X-left, and Y-rearward with the helmet facing −Y. The SAPI asset is posed twice, as a front plate at [0, -0.12, -0.5] and as a back plate mirrored about Z at [0, 0.18, -0.5], 0.3 m behind the front plate with its concave face toward the body.
Parameters. All of them in config.yaml, grouped and commented: the blast cube’s center, side length, samples per axis, and minimum stand-off distance; ray spacing in arc-minutes, max bounces, max trace length, and fan mode; the Friedlander constants, including the initial peak pressure (default 50,000 kPa) and the per-bounce reflection loss; the sensor diameter; the mesh paths; the compute device; and the output path. config.py validates the whole tree on load, with a specific message for each constraint it enforces.
Output. A timestamped blast-results-YYYY-MM-DD-HHMM.csv with one row per blast point and columns bp_x, bp_y, bp_z, full_armor, front_plate_helmet, helmet_only, vest_only, no_armor. The files committed in data/ predate the back plate, so they carry four armor columns and their full_armor corresponds to today’s front_plate_helmet.

Performance notes

A 15×15×15 cube, 3,352 points after the min-distance filter, takes about 1.2 seconds end to end on a DGX Spark (GB10). That figure holds from 8 arc-minutes down to 1 arc-minute, where each blast point casts a 150×150 fan (warp-times.log).
Rays are generated in the kernel, so memory stays trivial: about 0.5 GB even at 9,000 blast points and 2 arc-minute resolution.
Warp’s mesh queries run in float32 while the per-explosion impulse accumulates in float64, so expect five or six significant digits against a float64 reference. The last couple of CSV decimals, and cone-edge rays, can differ.
Warp uses a software BVH on the CUDA cores, not the GB10’s RT cores. The RT-core path (Mitsuba 3 / Dr.Jit) is worth another 3 to 9x if peak throughput at very fine resolution ever matters.
compute.device: cpu selects the Warp CPU backend, which is also where a run lands automatically when no CUDA device is visible.