README.md

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
# Pressure Wave

Simulates blast overpressure and impulse on a spherical sensor (a human brain)
placed inside MICH-helmet and SAPI-plate geometry (front and back torso plates),
using Friedlander waveforms and ray tracing with specular reflection. It is a
simplified model: it does **not** approximate human body composition around the
sensor.

The numerical core runs on **[NVIDIA Warp](https://nvidia.github.io/warp/)** — one
GPU thread per ray, each marching its own bounce path with real early-exit over a
hardware BVH (`wp.mesh_query_ray`). Rays are cast toward the sensor for every
point in a configurable blast cube, through five armor configurations
(`full_armor` = helmet + both plates, `front_plate_helmet` = helmet + front
plate, `helmet_only`, `vest_only` = front plate, `no_armor`); the
solid-angle-weighted total impulse is written to CSV.

# Getting Started

Environment management is via [`uv`](https://docs.astral.sh/uv/); all simulation
parameters live in a YAML file (`config.yaml`).

## Install

Host dependencies (NumPy, trimesh, PyYAML, SciPy) are the same everywhere:

    $ uv sync

The GPU backend is NVIDIA Warp, installed via the `cuda` extra:

    $ uv sync --extra cuda

**On the DGX Spark / GB10** (aarch64 + CUDA 13), if the PyPI wheel (CUDA 12.9)
doesn't cover your stack, install the official aarch64 CUDA-13 wheel directly:

    $ uv pip install \
        https://github.com/NVIDIA/warp/releases/download/v1.15.0/warp_lang-1.15.0+cu13-py3-none-manylinux_2_34_aarch64.whl

**Smoke-test Warp before anything else** — the GB10 is new hardware:

    $ uv run python -c "import warp as wp; wp.init(); print(wp.get_cuda_devices())"
    [Device('cuda:0')]

An empty list means Warp fell back to CPU (no CUDA device seen); set
`compute.device: cpu` to run on CPU deliberately.

## Run

    $ uv run simulate_blasts.py --config config.yaml

Every parameter lives in the YAML file. Override individual keys on the command
line with repeatable dotted `--set` flags:

    $ uv run simulate_blasts.py --config config.yaml \
          --set blast.cube_segments=15 \
          --set blast.cube_extents_m=20 \
          --set rays.resolution_arcmin=2

## Testing

    $ uv run pytest

The suite checks the physical invariants and cross-checks the `no_armor` column
against an independent NumPy computation. It runs on whatever `compute.device`
resolves to — the GPU, or the Warp CPU backend as a fallback.

## Configuration

See `config.yaml` for the full, commented schema. The knobs, by group:

| Group | Key | Meaning |
|-------|-----|---------|
| `blast` | `cube_center_m`, `cube_extents_m`, `cube_segments`, `min_dist_m` | blast-cube center, **side length**, **samples per axis (sampling density)**, and the minimum blast-to-sensor distance |
| `rays` | `resolution_arcmin`, `max_bounces`, `max_length_m`, `fan_mode` | **angular ray spacing (ray density)**, reflection depth, trace length, and whether the ray fan spans the sensor disc (`sensor_cone`) or widens to the armor (`geometry`) |
| `physics` | `initial_peak_pressure_kpa`, `positive_phase_duration_s`, `decay_alpha`, `reflection_loss` | Friedlander/blast constants |
| `sensor` | `diameter_m` | sensor sphere diameter (radius = diameter / 2) |
| `geometry` | `helmet`, `vest` | mesh paths |
| `compute` | `device` | `cuda` (GPU) or `cpu` (Warp CPU backend, for testing) |
| `output` | `file` | output CSV path (`null` → timestamped) |

Number of explosions = `cube_segments³` (minus points closer than `min_dist_m`).
Output CSV columns: `bp_x, bp_y, bp_z, full_armor, front_plate_helmet,
helmet_only, vest_only, no_armor`.

> **Data-compatibility note:** CSVs generated before the back plate was added
> (including those committed in `data/`) have four armor columns, and their
> `full_armor` corresponds to today's `front_plate_helmet` (helmet + front
> plate only).

## Edge diffraction (optional, off by default)

Pure ray tracing casts sharp geometric shadows: a blast point the armor occludes
contributes exactly zero. But the pulse is acoustically *large* — the positive
phase spans `c·t_d ≈ 3.4 m`, roughly ten times the helmet — so real blast energy
diffracts around the shell into those shadows. Setting `diffraction.enabled: true`
(or `--set diffraction.enabled=true`) adds **first-order edge diffraction** so
the geometric shadows fill in physically, which is primarily a robustness check:
does the armor-configuration ranking survive the model's single largest
omission?

The method is **Biot–Tolstoy–Medwin (BTM)**, not GTD/UTD. Impulse is the
zero-frequency (DC) component of pressure, and the high-frequency GTD/UTD
diffraction coefficient *diverges* at DC (it scales `1/√k`); BTM's time-domain
response is integrable, giving a finite DC edge gain. The diffracted impulse is
the incident impulse times that gain, folded into the same accumulator as the
direct/reflected rays so **geometric-optics + diffraction is continuous across
the shadow boundary** by construction.

The DC coefficient is validated in `tests/test_diffraction.py` against the exact
static rigid-wedge Green's function (a Legendre-Q modal series, itself checked
against image sums to 1e-9), across wedge angles including the screen (the rim
case), and independently by the shadow-boundary continuity condition. The
coupling constant `1/(16π²)` is fixed by those checks, not fitted.

Because MICH.obj is a smooth watertight shell with no rim edge, the helmet's
diffracting edge is sourced procedurally — the intersection of the shell with a
horizontal (constant-Z) plane (`diffraction.edge_mode: rim_plane`); the torso
plates use their perimeter feature edges. See the `diffraction:` block in
`config.yaml`.

**Scope and caveats.** First-order only (no double diffraction, no
reflection×diffraction); *edge* diffraction, not smooth-surface creeping waves;
impulse-only (DC), linear acoustics (no Mach stem). The diffraction term is
signed, so where the discrete ray-traced geometric field under-resolves a
near-zero-field transition cell the naive sum can dip slightly negative (worst
case a few percent of the data scale, at cells whose geometric field is ~0); the
physical impulse is clamped at zero and the count logged, and the effect shrinks
with finer `rays.resolution_arcmin`. Treat diffraction-on results as a
comparative robustness bound, not absolute dosimetry.

## Precision and performance notes

- Warp's mesh queries run in **float32**; the per-explosion impulse is accumulated
  in **float64**. Expect ~5–6 significant digits vs a float64 reference — the last
  couple of CSV decimals, and cone-edge rays, can differ.
- Memory is trivial (~0.5 GB even at 9k points / 2 arcmin) — rays are generated
  in-kernel, never materialized.
- Warp uses a **software BVH on CUDA cores** (not the GB10's RT cores). If peak
  throughput at very fine resolution is needed, the RT-core path is
  [Mitsuba 3 / Dr.Jit](https://www.mitsuba-renderer.org/) — a further ~3–9×.

# Analysis

`analysis.R` performs an ANOVA followed by paired t-tests with a Holm adjustment:

    $ R -q -s --no-save < analysis.R

`visualize.R` contains ggplot/rgl snippets for looking at the data (source
`analysis.R` first — it prepares the data frames those snippets use).

`Exploration.ipynb` is an interactive walkthrough of the full pipeline —
configuration → geometry → ray-fan sampling → Friedlander physics → a small
Warp run with plots and an independent NumPy cross-check. It runs in a couple
of minutes on the Warp CPU fallback (no GPU required):

    $ uv sync --group notebook --extra cuda
    $ uv run jupyter lab Exploration.ipynb

# Project Notes

The geometry is oriented so the helmet faces the negative Y direction. The
environment is Z-up, X-left, Y-rearward. The front SAPI plate sits at
`[0, -0.12, -0.5]` (just in front of and below the origin sensor, torso
height); the back plate is the same asset mirrored about Z at `[0, 0.18, -0.5]`
— 0.3 m behind the front plate, concave face toward the body.

# License

Copyright (C) 2025 SWGY, Inc

This project — including source code, the geometry assets in `assets/`, the
simulation outputs in `data/`, and the `Exploration.ipynb` notebook — is
licensed under the GNU Affero General Public License, version 3 or (at your
option) any later version. See the [`LICENSE`](./LICENSE) file for the full
text, or <https://www.gnu.org/licenses/agpl-3.0.html>.

`SPDX-License-Identifier: AGPL-3.0-or-later`