{
"cells": [
{
"id": "title-md",
"cell_type": "markdown",
"metadata": {},
"source": [
"# Exploring the blast-pressure model\n",
"\n",
"This notebook is an interactive walkthrough of the simulation in this repository:\n",
"a point blast produces a **Friedlander overpressure wave**, which we propagate by\n",
"**ray tracing with specular reflections** against MICH-helmet and SAPI-vest\n",
"meshes (front and back torso plates), accumulating the **total impulse**\n",
"delivered to a spherical sensor (a brain proxy) at the origin. Every blast\n",
"point is evaluated under five armor configurations — `full_armor` (helmet +\n",
"both plates), `front_plate_helmet` (helmet + front plate), `helmet_only`,\n",
"`vest_only` (front plate), `no_armor` — and the numerical core runs on\n",
"[NVIDIA Warp](https://nvidia.github.io/warp/), one thread per ray.\n",
"\n",
"The arc:\n",
"\n",
"1. **Configuration** — the YAML schema and the override mechanism\n",
"2. **Geometry** — the armor meshes and the world frame, in 3D\n",
"3. **Ray-fan sampling** — how rays are aimed (tangent cones, solid angles)\n",
"4. **Friedlander physics** — the waveform and the precomputed impulse constant\n",
"5. **Running on Warp** — what the kernel does, then a real (small) run\n",
"6. **Results** — heatmaps and comparisons across the four armor configurations\n",
"7. **Validation** — an independent NumPy cross-check of the engine\n",
"8. **Edge diffraction** — filling the geometric shadow (optional)\n",
"9. **Where to go from here**\n",
"\n",
"It runs end-to-end without a GPU: Warp falls back to its CPU backend, and the\n",
"demo grid below completes in a couple of minutes (most of it one-time kernel\n",
"compilation). With a CUDA device the same run takes seconds."
]
},
{
"id": "license-md",
"cell_type": "markdown",
"metadata": {},
"source": [
"<!--\n",
"SPDX-License-Identifier: AGPL-3.0-or-later\n",
"Copyright (C) 2025 SWGY, Inc\n",
"-->\n",
"\n",
"# License\n",
"\n",
"Copyright (C) 2025 SWGY, Inc\n",
"\n",
"This notebook is part of the pressure-wave project and is licensed under\n",
"the GNU Affero General Public License v3.0 or later. See the `LICENSE`\n",
"file at the repository root for the full license text, or\n",
"<https://www.gnu.org/licenses/agpl-3.0.html>."
]
},
{
"id": "setup-md",
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup\n",
"\n",
"From the repository root:\n",
"\n",
"```sh\n",
"uv sync --group notebook --extra cuda\n",
"uv run jupyter lab Exploration.ipynb\n",
"```\n",
"\n",
"`--extra cuda` installs NVIDIA Warp (`warp-lang`). It is needed **even on\n",
"machines without a GPU** — the same wheel contains both the CUDA and CPU\n",
"backends, and `warp_core.configure()` falls back to the CPU automatically when\n",
"no CUDA device is present. On the DGX Spark / GB10, see the README for the\n",
"CUDA-13 aarch64 wheel."
]
},
{
"id": "imports",
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"import logging\n",
"import time\n",
"\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import trimesh\n",
"\n",
"import config as config_mod\n",
"import geometry\n",
"import sampling\n",
"import warp_core\n",
"\n",
"logging.basicConfig(level=logging.INFO,\n",
" format=\"%(asctime)s - %(levelname)s - %(message)s\")\n",
"\n",
"# One fixed color per armor configuration, used consistently in every plot.\n",
"# full_armor and front_plate_helmet share the blue \"helmet + plate(s)\" family\n",
"# hue; plots tell them apart by marker fill and size.\n",
"CONFIG_ORDER = [\"full_armor\", \"front_plate_helmet\", \"helmet_only\",\n",
" \"vest_only\", \"no_armor\"]\n",
"CONFIG_COLORS = {\n",
" \"full_armor\": \"#2a78d6\", # blue (hollow marks)\n",
" \"front_plate_helmet\": \"#2a78d6\", # blue (filled marks)\n",
" \"helmet_only\": \"#eb6834\", # orange\n",
" \"vest_only\": \"#1baf7a\", # aqua\n",
" \"no_armor\": \"#4a3aa7\", # violet\n",
"}\n",
"\n",
"plt.rcParams.update({\n",
" \"axes.spines.top\": False,\n",
" \"axes.spines.right\": False,\n",
" \"axes.grid\": True,\n",
" \"grid.color\": \"#e1e0d9\",\n",
" \"grid.linewidth\": 0.8,\n",
" \"axes.edgecolor\": \"#c3c2b7\",\n",
" \"axes.labelcolor\": \"#52514e\",\n",
" \"axes.titlecolor\": \"#0b0b0b\",\n",
" \"xtick.color\": \"#898781\",\n",
" \"ytick.color\": \"#898781\",\n",
" \"figure.autolayout\": True,\n",
"})"
]
},
{
"id": "config-md",
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1. Configuration\n",
"\n",
"Every simulation parameter lives in `config.yaml` — blast-cube placement, ray\n",
"density, Friedlander constants, mesh paths, compute device. `load_config`\n",
"applies *dotted overrides* on top of the file, and they are exactly the strings\n",
"the CLI accepts via repeatable `--set` flags, so anything you try here\n",
"translates directly to a batch run.\n",
"\n",
"This demo shrinks the problem so the whole notebook re-runs quickly on a CPU:\n",
"\n",
"| Override | Why |\n",
"|---|---|\n",
"| `blast.cube_segments=5` | 5 samples per axis → 125 blast points; an **odd** count means the $z=0$ plane is sampled exactly, which gives us a clean heatmap slice later |\n",
"| `blast.cube_extents_m=6.0` | a 6 m cube around the default center `[0, -4, 0]` |\n",
"| `blast.min_dist_m=0.9` | keeps the nearest grid point `(0, -1, 0)`, exactly 1 m from the sensor |\n",
"| `rays.resolution_arcmin=15` | coarse enough to be fast, fine enough that every blast point still resolves the sensor cone with several rays per axis |\n",
"\n",
"Everything not listed keeps its `config.yaml` default — notably\n",
"$p_0 = 50\\,000$ kPa, `max_bounces=4`, and the 0.111 m sensor diameter."
]
},
{
"id": "config-code",
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"DEMO_OVERRIDES = [\n",
" \"blast.cube_segments=5\",\n",
" \"blast.cube_extents_m=6.0\",\n",
" \"blast.min_dist_m=0.9\",\n",
" \"rays.resolution_arcmin=15\",\n",
"]\n",
"cfg = config_mod.load_config(\"config.yaml\", DEMO_OVERRIDES)\n",
"\n",
"print(f\"blast cube: center={cfg.blast.cube_center_m} \"\n",
" f\"extents={cfg.blast.cube_extents_m} m \"\n",
" f\"segments={cfg.blast.cube_segments} -> {cfg.blast.cube_segments ** 3} points\")\n",
"print(f\"rays: {cfg.rays.resolution_arcmin:.0f} arcmin spacing, \"\n",
" f\"max {cfg.rays.max_bounces} bounces, max length {cfg.rays.max_length_m} m, \"\n",
" f\"fan_mode={cfg.rays.fan_mode}\")\n",
"print(f\"physics: p0={cfg.physics.initial_peak_pressure_kpa:.0f} kPa, \"\n",
" f\"t_d={cfg.physics.positive_phase_duration_s} s, \"\n",
" f\"alpha={cfg.physics.decay_alpha}, \"\n",
" f\"reflection_loss={cfg.physics.reflection_loss}\")\n",
"print(f\"sensor: radius {cfg.sensor.radius_m:.4f} m, at the origin\")"
]
},
{
"id": "geometry-md",
"cell_type": "markdown",
"metadata": {},
"source": [
"## 2. Geometry\n",
"\n",
"The world frame puts the **sensor at the origin** with **Z up, X left,\n",
"Y rearward** — the helmet faces $-Y$. The loaders in `geometry.py` pose the raw\n",
"OBJ assets into that frame once, at load time:\n",
"\n",
"- `load_helmet` rotates the MICH shell $+90°$ about X (asset frame → Z-up) and\n",
" scales it by 0.5, seating it around the origin sensor;\n",
"- `load_vest` flips the SAPI plate $180°$ about X and hangs it at\n",
" `[0, -0.12, -0.5]` — just in front of and below the sensor, at the torso;\n",
"- `load_back_plate` poses the *same* SAPI asset as the back plate: the same X\n",
" flip, then $180°$ about Z so the concave face looks forward (toward the\n",
" body), at `[0, +0.18, -0.5]` — 0.3 m behind the front plate, same height.\n",
" The facing matters even though ray hits are double-sided: curvature decides\n",
" whether reflections focus toward the wearer or scatter away.\n",
"\n",
"Both loaders call `fix_normals()` for consistent face orientation. (The\n",
"reflection formula in the kernel is sign-invariant in the normal, so\n",
"consistency is all that matters.)"
]
},
{
"id": "geometry-code",
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"helmet = geometry.load_helmet(cfg.geometry.helmet)\n",
"vest = geometry.load_vest(cfg.geometry.vest)\n",
"back = geometry.load_back_plate(cfg.geometry.vest)\n",
"sensor = sampling.Sensor(origin=np.zeros(3), radius=cfg.sensor.radius_m)\n",
"\n",
"for name, mesh in [(\"helmet\", helmet), (\"vest\", vest), (\"back\", back)]:\n",
" lo, hi = mesh.bounds\n",
" print(f\"{name:6s} {len(mesh.faces):5d} faces \"\n",
" f\"x [{lo[0]:+.2f}, {hi[0]:+.2f}] \"\n",
" f\"y [{lo[1]:+.2f}, {hi[1]:+.2f}] \"\n",
" f\"z [{lo[2]:+.2f}, {hi[2]:+.2f}] (m)\")"
]
},
{
"id": "scene1-md",
"cell_type": "markdown",
"metadata": {},
"source": [
"### The stage\n",
"\n",
"Everything a ray can interact with: the helmet shell around the red sensor\n",
"sphere, and the front and back plates bracketing the torso. The yellow marker\n",
"is the center of the blast cube, 4 m in front of the face. Drag to orbit.\n",
"\n",
"*The viewer needs internet access to render (three.js from a CDN); the notebook\n",
"executes fine without it.*"
]
},
{
"id": "scene1-code",
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"def tinted(mesh, rgba):\n",
" m = mesh.copy()\n",
" m.visual.face_colors = rgba\n",
" return m\n",
"\n",
"\n",
"sensor_ball = trimesh.creation.icosphere(subdivisions=3, radius=sensor.radius)\n",
"sensor_ball.visual.face_colors = [227, 73, 72, 255] # the brain-proxy sensor\n",
"\n",
"blast_marker = trimesh.creation.icosphere(subdivisions=2, radius=0.08)\n",
"blast_marker.apply_translation(cfg.blast.cube_center_m)\n",
"blast_marker.visual.face_colors = [237, 161, 0, 255] # blast-cube center\n",
"\n",
"trimesh.Scene([\n",
" tinted(helmet, [42, 120, 214, 160]), # translucent helmet\n",
" tinted(vest, [27, 175, 122, 200]), # front plate\n",
" tinted(back, [27, 175, 122, 200]), # back plate\n",
" sensor_ball,\n",
" blast_marker,\n",
"]).show(viewer=\"jupyter\")"
]
},
{
"id": "sampling-md",
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3. Aiming the rays: cone sampling\n",
"\n",
"A blast radiates in every direction, but almost none of those directions can\n",
"ever reach a 5.5 cm sensor several meters away. So instead of sampling the full\n",
"sphere, each blast point only traces the **cone of directions that subtends the\n",
"sensor**. From distance $d$, a sphere of radius $r$ is seen under the half-angle\n",
"\n",
"$$\\theta = \\arcsin(r / d)$$\n",
"\n",
"(the tangent ray touches the sphere where it is perpendicular to the radius, so\n",
"$\\sin\\theta = r/d$).\n",
"\n",
"`sampling.build_grid` prepares, for every blast point, a compact description of\n",
"that cone — the device builds the actual rays later, in-kernel:\n",
"\n",
"1. an orthonormal **frame** `(forward, right, up)` with *forward* aimed at the\n",
" sensor;\n",
"2. an **azimuth × elevation lattice** spanning the cone at\n",
" `rays.resolution_arcmin` spacing (1 arcmin = 1/60°);\n",
"3. all fans **padded** to the largest lattice across blast points, with\n",
" `valid_az` / `valid_el` masks marking the real rays — the Warp launch is one\n",
" rectangular grid of threads, one per *(explosion, ray)*, and a padded thread\n",
" simply exits on its first instruction.\n",
"\n",
"Each ray stands for a small patch of solid angle,\n",
"$d\\Omega = \\cos(el)\\cdot\\mathrm{res}^2$ steradians — the weight used when its\n",
"impulse is accumulated.\n",
"\n",
"Why sample only the sensor cone when rays can *bounce*? In the default\n",
"`fan_mode: sensor_cone`, reflected paths are still explored — the cone rays\n",
"themselves reflect off whatever armor they graze. The alternative\n",
"`fan_mode: geometry` widens each fan to the armor's full angular extent, tracing\n",
"many more rays; useful when reflection-rich paths need denser sampling."
]
},
{
"id": "grid-code",
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"grid = sampling.build_grid(cfg, sensor) # sensor_cone mode needs no mesh vertices\n",
"\n",
"E = grid.blast_os.shape[0]\n",
"gaz, gel = grid.az_off.shape[1], grid.el_off.shape[1]\n",
"print(f\"{E} blast points survive the min-distance filter\")\n",
"print(f\"padded fan: {gaz} az x {gel} el = {gaz * gel} threads per (explosion, config)\")\n",
"print(f\"angular spacing: {np.degrees(grid.res_rad) * 60:.0f} arcmin\")\n",
"\n",
"# The cube-center blast point, used as the running example below.\n",
"center = np.asarray(cfg.blast.cube_center_m, dtype=float)\n",
"e0 = int(np.flatnonzero(np.all(np.isclose(grid.blast_os, center), axis=1))[0])\n",
"d0 = float(np.linalg.norm(center - sensor.origin))\n",
"theta0 = np.arcsin(sensor.radius / d0)\n",
"n_az0 = int(grid.valid_az[e0].sum())\n",
"n_el0 = int(grid.valid_el[e0].sum())\n",
"print(f\"cube center: d = {d0:.0f} m -> theta = {np.degrees(theta0):.2f} deg \"\n",
" f\"= {np.degrees(theta0) * 60:.0f} arcmin -> {n_az0} x {n_el0} rays\")"
]
},
{
"id": "cone-diagram",
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"from matplotlib.patches import Arc, Circle\n",
"\n",
"fig, ax = plt.subplots(figsize=(9, 3.4))\n",
"d_s, r_s = 4.0, 0.7 # schematic: sensor radius exaggerated\n",
"theta_s = np.arcsin(r_s / d_s)\n",
"tangent_len = np.sqrt(d_s ** 2 - r_s ** 2)\n",
"\n",
"ax.add_patch(Circle((d_s, 0.0), r_s, facecolor=\"#cde2fb\",\n",
" edgecolor=\"#2a78d6\", lw=2))\n",
"ax.plot([0.0], [0.0], marker=\"*\", ms=16, color=\"#eda100\", clip_on=False)\n",
"for s in (1.0, -1.0):\n",
" ax.plot([0.0, tangent_len * np.cos(s * theta_s)],\n",
" [0.0, tangent_len * np.sin(s * theta_s)], color=\"#52514e\", lw=1.5)\n",
"ax.plot([0.0, d_s], [0.0, 0.0], color=\"#898781\", lw=1.0, ls=\"--\")\n",
"ax.add_patch(Arc((0, 0), 2.6, 2.6, theta1=0.0, theta2=np.degrees(theta_s),\n",
" color=\"#0b0b0b\", lw=1.2))\n",
"ax.annotate(r\"$\\theta = \\arcsin(r/d)$\", xy=(1.42, 0.2),\n",
" color=\"#0b0b0b\", fontsize=12)\n",
"ax.annotate(\"blast point\", xy=(-0.15, -0.42), color=\"#52514e\")\n",
"ax.annotate(\"sensor (radius $r$)\", xy=(d_s - 0.55, -r_s - 0.45), color=\"#52514e\")\n",
"ax.annotate(\"$d$\", xy=(d_s / 2, -0.32), color=\"#898781\", fontsize=12)\n",
"ax.set_xlim(-0.4, d_s + 1.1)\n",
"ax.set_ylim(-1.5, 1.5)\n",
"ax.set_aspect(\"equal\")\n",
"ax.axis(\"off\")\n",
"ax.set_title(\"The tangent cone from a blast point to the sensor \"\n",
" \"(radius exaggerated)\")\n",
"plt.show()"
]
},
{
"id": "fan-md",
"cell_type": "markdown",
"metadata": {},
"source": [
"### The actual fan\n",
"\n",
"The cube-center blast point sees the sensor under $\\approx 0.8°$, an\n",
"$8\\times8$ lattice at 15 arcmin — small enough to draw every ray, rebuilt here\n",
"with the same frame + offsets arithmetic the Warp kernel uses. The gray points\n",
"are all 125 blast points of the demo grid; the orange bundle converges on the\n",
"sensor and continues 1 m past it."
]
},
{
"id": "scene2-code",
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"bp = grid.blast_os[e0]\n",
"fwd, right, up = grid.frames[e0]\n",
"el = grid.el_off[e0][grid.valid_el[e0]][:, None] # (Nel, 1)\n",
"az = grid.az_off[e0][grid.valid_az[e0]][None, :] # (1, Naz)\n",
"cos_el = np.cos(el)\n",
"dirs = (fwd * (cos_el * np.cos(az))[..., None] # same formula as the kernel\n",
" + right * (cos_el * np.sin(az))[..., None]\n",
" + up * np.sin(el)[..., None]).reshape(-1, 3)\n",
"\n",
"segments = np.stack([np.broadcast_to(bp, dirs.shape),\n",
" bp + dirs * (d0 + 1.0)], axis=1)\n",
"fan = trimesh.load_path(segments)\n",
"fan.colors = np.tile([235, 104, 52, 140], (len(fan.entities), 1))\n",
"cloud = trimesh.points.PointCloud(\n",
" grid.blast_os, colors=np.tile([137, 135, 129, 255], (E, 1)))\n",
"\n",
"print(f\"{dirs.shape[0]} rays in the fan from the cube-center blast point\")\n",
"trimesh.Scene([\n",
" tinted(helmet, [42, 120, 214, 160]),\n",
" tinted(vest, [27, 175, 122, 200]),\n",
" tinted(back, [27, 175, 122, 200]),\n",
" sensor_ball,\n",
" cloud,\n",
" fan,\n",
"]).show(viewer=\"jupyter\")"
]
},
{
"id": "friedlander-md",
"cell_type": "markdown",
"metadata": {},
"source": [
"## 4. The Friedlander waveform\n",
"\n",
"The classic idealized blast overpressure at a fixed point is the **Friedlander\n",
"waveform**: an instantaneous jump to the peak $p_0$, then a decaying positive\n",
"phase of duration $t_d$,\n",
"\n",
"$$p(t) = p_0\\,\\Bigl(1 - \\frac{t}{t_d}\\Bigr)\\,e^{-\\alpha t / t_d},$$\n",
"\n",
"which this model clips at zero (the small negative phase is ignored).\n",
"\n",
"The kernel never integrates this curve. **Impulse is linear in peak pressure**,\n",
"so the shape integral\n",
"\n",
"$$K = \\int_0^{t_d} \\Bigl(1 - \\tfrac{t}{t_d}\\Bigr)\\,e^{-\\alpha t/t_d}\\,dt$$\n",
"\n",
"is precomputed once (`sampling.precompute_waveform_constant`) — for\n",
"$\\alpha = 1$ it even has the closed form $K = t_d / e$. Each ray that reaches\n",
"the sensor then contributes\n",
"\n",
"$$\\text{impulse} = \\underbrace{\\text{refl}}_{0.8^{\\,b}\\ \\text{after}\\ b\\\n",
"\\text{bounces}} \\cdot \\underbrace{\\frac{p_0}{d_{\\text{eff}}^2}}_{\\text{inverse\n",
"square}} \\cdot\\; K,$$\n",
"\n",
"where $d_\\text{eff}$ is the ray's **total path length** including every bounce\n",
"leg, and $p_0$ is referenced at 1 m. Units: $p$ in kPa and $K$ in s give a\n",
"per-ray impulse in kPa·s; after solid-angle weighting the per-explosion totals\n",
"are in kPa·s·sr."
]
},
{
"id": "waveform-code",
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"p0 = cfg.physics.initial_peak_pressure_kpa\n",
"t_d = cfg.physics.positive_phase_duration_s\n",
"alpha = cfg.physics.decay_alpha\n",
"\n",
"t = np.linspace(0.0, 3.0 * t_d, 1000)\n",
"p = np.clip(p0 * (1.0 - t / t_d) * np.exp(-alpha * t / t_d), 0.0, None)\n",
"\n",
"k_const = sampling.precompute_waveform_constant(t_d, alpha)\n",
"\n",
"fig, ax = plt.subplots(figsize=(8, 4))\n",
"ax.plot(t * 1e3, p, color=\"#2a78d6\", lw=2)\n",
"ax.fill_between(t * 1e3, p, 0.0, color=\"#2a78d6\", alpha=0.15)\n",
"ax.axvline(t_d * 1e3, color=\"#c3c2b7\", lw=1, ls=\"--\")\n",
"ax.annotate(\"$t_d$\", xy=(t_d * 1e3, p0 * 0.04), xytext=(t_d * 1e3 + 0.6, p0 * 0.04),\n",
" color=\"#898781\", fontsize=12)\n",
"ax.annotate(f\"shaded area = $p_0 \\\\cdot K$\\n$K$ = {k_const:.5f} s\",\n",
" xy=(t_d * 0.32 * 1e3, p0 * 0.28), color=\"#52514e\")\n",
"ax.set_xlabel(\"time (ms)\")\n",
"ax.set_ylabel(\"overpressure (kPa)\")\n",
"ax.set_title(f\"Friedlander waveform at 1 m ($p_0$ = {p0:,.0f} kPa)\")\n",
"plt.show()\n",
"\n",
"print(f\"K = {k_const:.6f} s closed form t_d/e = {t_d / np.e:.6f} s\")"
]
},
{
"id": "warp-md",
"cell_type": "markdown",
"metadata": {},
"source": [
"## 5. Running it on Warp\n",
"\n",
"`warp_core.simulate` launches one kernel per armor configuration, with **one\n",
"thread per (explosion, ray)** — for this demo, $125 \\times 27 \\times 27 \\approx\n",
"91\\,000$ threads per launch. Each thread:\n",
"\n",
"1. builds its ray from the per-explosion frame and its az/el lattice offsets\n",
" (rays are never materialized on the host — padded lattice slots exit\n",
" immediately);\n",
"2. walks a bounce loop with a real early-exit `break`:\n",
" - queries this configuration's mesh(es) with `wp.mesh_query_ray` — a\n",
" **BVH-accelerated nearest-hit** lookup, double-sided;\n",
" - solves the quadratic for the **sensor sphere** along the same segment;\n",
" - if the sensor comes first: records\n",
" $\\text{refl} \\cdot p_0/d_\\text{eff}^2 \\cdot K$ and stops;\n",
" - if armor comes first: advances to the surface (offset by $10^{-6}$ m to\n",
" avoid re-hitting it), reflects specularly\n",
" ($\\mathbf{d}' = \\mathbf{d} - 2(\\mathbf{d}\\cdot\\mathbf{n})\\mathbf{n}$), and\n",
" multiplies `refl` by `reflection_loss`;\n",
"3. atomically adds $\\text{impulse}\\cdot d\\Omega$ into the per-explosion\n",
" accumulator.\n",
"\n",
"Mesh queries run in **float32**; the accumulator is **float64** — expect ~5–6\n",
"significant digits of agreement with a float64 reference (we quantify this in\n",
"§7).\n",
"\n",
"The five configurations differ only in which meshes are queried — `no_armor`\n",
"queries none, so only direct line-of-sight rays count. The kernel takes at most\n",
"two meshes per launch, so `full_armor` merges the front and back plates into a\n",
"single BVH (nearest-hit over the union is the same either way). Note there is\n",
"**no guaranteed ordering** between configurations: a helmet shadows the sensor\n",
"from some directions but *focuses reflections onto it* from others, so armor\n",
"can amplify as well as attenuate."
]
},
{
"id": "simulate-code",
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"device = warp_core.configure(cfg.compute.device) # \"cuda\" if available, else CPU\n",
"print(f\"Warp device: {device}\")\n",
"\n",
"t0 = time.perf_counter()\n",
"results = warp_core.simulate(cfg, helmet, vest, back, sensor, grid, k_const,\n",
" device)\n",
"elapsed = time.perf_counter() - t0\n",
"print(f\"5 configs x {E} blast points in {elapsed:.1f} s on '{device}' \"\n",
" f\"(first call includes one-time kernel compilation)\\n\")\n",
"\n",
"print(f\"{'config':>20s} {'min':>10s} {'median':>10s} {'max':>10s} (kPa*s*sr)\")\n",
"for name in CONFIG_ORDER:\n",
" v = results[name]\n",
" print(f\"{name:>20s} {v.min():10.4f} {np.median(v):10.4f} {v.max():10.4f}\")"
]
},
{
"id": "results-md",
"cell_type": "markdown",
"metadata": {},
"source": [
"## 6. What the armor does\n",
"\n",
"Below, the $z = 0$ slice of the blast cube — 25 points in the sensor's\n",
"horizontal plane — for each configuration. The sensor sits at the origin\n",
"(marked `+`); $-y$ is in front of the face.\n",
"\n",
"The color scale is **logarithmic** and shared across panels, because total\n",
"impulse falls off as $\\sim 1/d^4$: the wave amplitude decays as $1/d^2$\n",
"*and* the captured solid angle shrinks as $(r/d)^2$. Across 1–8 m that spans\n",
"more than three decades. (A cell where every ray was blocked before reaching\n",
"the sensor would drop off the log scale and show gray.)\n",
"\n",
"Two things to look for:\n",
"\n",
"- **the helmet's shadow**: cells in front of the face are visibly lighter\n",
" (roughly 2–3× less impulse) with a helmet than without;\n",
"- **the panels cluster** — the torso plates hang half a meter *below* the\n",
" sensor and intercept almost nothing arriving from this frontal cube (the\n",
" back plate, behind the sensor, is outside every sensor cone entirely), so\n",
" the helmet does all the work: `full_armor` ≈ `front_plate_helmet` ≈\n",
" `helmet_only`, and `vest_only` ≈ `no_armor`. The scatter below nests the\n",
" markers so the coincidences are visible instead of overplotted. Where the\n",
" plates *do* matter is reflection — sample all around the wearer with\n",
" `fan_mode: geometry` (see §9) and they reappear as small amplifiers."
]
},
{
"id": "heatmap-code",
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"from matplotlib.colors import LinearSegmentedColormap, LogNorm\n",
"\n",
"# Single-hue sequential ramp (light -> dark = small -> large impulse)\n",
"BLUES = LinearSegmentedColormap.from_list(\"pw_blues\", [\n",
" \"#cde2fb\", \"#9ec5f4\", \"#6da7ec\", \"#3987e5\", \"#256abf\", \"#184f95\", \"#0d366b\"])\n",
"BLUES.set_bad(\"#f0efec\") # zero-impulse cells drop off the log scale\n",
"\n",
"slice_mask = np.isclose(grid.blast_os[:, 2], 0.0)\n",
"xs = np.unique(grid.blast_os[slice_mask, 0])\n",
"ys = np.unique(grid.blast_os[slice_mask, 1])\n",
"\n",
"\n",
"def pivot(vals):\n",
" g = np.full((ys.size, xs.size), np.nan)\n",
" for (x, y, _), v in zip(grid.blast_os[slice_mask], vals[slice_mask]):\n",
" g[np.searchsorted(ys, y), np.searchsorted(xs, x)] = v\n",
" return np.ma.masked_where(~(g > 0), g)\n",
"\n",
"\n",
"positive = np.concatenate([results[n][slice_mask] for n in CONFIG_ORDER])\n",
"positive = positive[positive > 0]\n",
"norm = LogNorm(vmin=positive.min(), vmax=positive.max())\n",
"\n",
"fig, axes = plt.subplots(1, 5, figsize=(15.5, 4.2), sharey=True,\n",
" layout=\"constrained\")\n",
"for ax, name in zip(axes, CONFIG_ORDER):\n",
" mesh = ax.pcolormesh(xs, ys, pivot(results[name]), norm=norm, cmap=BLUES,\n",
" shading=\"nearest\")\n",
" ax.plot(0.0, 0.0, marker=\"+\", color=\"#0b0b0b\", ms=10, mew=2)\n",
" ax.set_title(name)\n",
" ax.set_xlabel(\"blast x (m)\")\n",
" ax.set_xlim(-3.9, 3.9)\n",
" ax.set_ylim(-7.9, 0.5)\n",
" ax.set_aspect(\"equal\")\n",
" ax.grid(False)\n",
"axes[0].set_ylabel(\"blast y (m) ($-y$ = in front of the face)\")\n",
"axes[0].annotate(\"sensor\", xy=(0.25, -0.15), color=\"#52514e\")\n",
"fig.colorbar(mesh, ax=axes, label=\"total impulse (kPa*s*sr), log scale\",\n",
" shrink=0.85)\n",
"plt.show()"
]
},
{
"id": "scatter-code",
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"dist = np.linalg.norm(grid.blast_os - sensor.origin, axis=1)\n",
"\n",
"# marker, size, filled? -- hollow marks nest over their filled partners, so\n",
"# coinciding configurations stay visible instead of overplotting each other:\n",
"# full_armor (ring) ~ front_plate_helmet (dot) ~ helmet_only (square), and\n",
"# no_armor (diamond) ~ vest_only (triangle).\n",
"STYLE = {\n",
" \"full_armor\": (\"o\", 11, False),\n",
" \"front_plate_helmet\": (\"o\", 5, True),\n",
" \"helmet_only\": (\"s\", 8, False),\n",
" \"vest_only\": (\"^\", 5, True),\n",
" \"no_armor\": (\"D\", 8, False),\n",
"}\n",
"\n",
"fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4.6))\n",
"\n",
"# Impulse vs distance, all 125 points, all five configurations\n",
"for name in CONFIG_ORDER:\n",
" v = results[name]\n",
" m = v > 0\n",
" marker, ms, filled = STYLE[name]\n",
" ax1.plot(dist[m], v[m], marker, ls=\"none\", ms=ms,\n",
" mfc=CONFIG_COLORS[name] if filled else \"none\",\n",
" mec=\"white\" if filled else CONFIG_COLORS[name],\n",
" mew=0.4 if filled else 1.2,\n",
" alpha=0.85, label=name)\n",
"anchor = int(np.argmax(results[\"no_armor\"]))\n",
"dd = np.geomspace(dist.min(), dist.max(), 50)\n",
"ax1.plot(dd, results[\"no_armor\"][anchor] * (dist[anchor] / dd) ** 4,\n",
" color=\"#898781\", lw=1, ls=\"--\")\n",
"ax1.annotate(r\"$\\sim 1/d^4$\", xy=(dd[28], results[\"no_armor\"][anchor]\n",
" * (dist[anchor] / dd[28]) ** 4), xytext=(4, 4),\n",
" textcoords=\"offset points\", color=\"#898781\", fontsize=12)\n",
"ax1.set_xscale(\"log\")\n",
"ax1.set_yscale(\"log\")\n",
"ax1.set_xlabel(\"blast distance (m)\")\n",
"ax1.set_ylabel(\"total impulse (kPa*s*sr)\")\n",
"ax1.set_title(\"Impulse vs distance (points with a nonzero path)\")\n",
"ax1.legend(framealpha=0.9)\n",
"\n",
"# Armor effect straight ahead: the x=0, z=0 column\n",
"col = np.isclose(grid.blast_os[:, 0], 0.0) & np.isclose(grid.blast_os[:, 2], 0.0)\n",
"col &= results[\"no_armor\"] > 0\n",
"order = np.argsort(-grid.blast_os[col, 1]) # near (y=-1) to far (y=-7)\n",
"front_d = -grid.blast_os[col, 1][order]\n",
"ratio = (results[\"full_armor\"][col] / results[\"no_armor\"][col])[order]\n",
"\n",
"ax2.plot(front_d, ratio, \"-o\", color=CONFIG_COLORS[\"full_armor\"], lw=2, ms=6,\n",
" mec=\"white\", mew=0.4)\n",
"ax2.axhline(1.0, color=\"#c3c2b7\", lw=1, ls=\"--\")\n",
"ax2.text(front_d[-1], 1.02, \"no effect\", color=\"#898781\", ha=\"right\")\n",
"ax2.set_xlabel(\"blast distance straight ahead (m)\")\n",
"ax2.set_ylabel(\"full_armor / no_armor impulse\")\n",
"ax2.set_title(\"Armor effect head-on: < 1 attenuates, > 1 amplifies\")\n",
"plt.show()"
]
},
{
"id": "crosscheck-md",
"cell_type": "markdown",
"metadata": {},
"source": [
"## 7. Validating against plain NumPy\n",
"\n",
"Good practice for any fast numerical engine: check it against an independent,\n",
"slow, *obvious* implementation. The `no_armor` column is perfect for this — no\n",
"meshes are involved, so the whole computation is \"cast the cone rays, intersect\n",
"each with a sphere, sum the weighted impulses\". Below it is rebuilt in ~20 lines\n",
"of NumPy that share **no code** with the Warp kernel (the test suite performs\n",
"the same cross-check in `tests/test_physics.py`).\n",
"\n",
"The tolerance is deliberate: Warp traces in float32 while the reference runs in\n",
"float64, and a ray *exactly* on the cone edge can flip between hit and miss —\n",
"so agreement is asserted against the overall data scale, not point by point."
]
},
{
"id": "crosscheck-code",
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"def no_armor_reference(cfg, sensor, grid, k_const):\n",
" \"Direct-path impulse at each blast point -- independent of warp_core.\"\n",
" p0 = cfg.physics.initial_peak_pressure_kpa\n",
" max_len = cfg.rays.max_length_m\n",
" res = grid.res_rad\n",
" r = sensor.radius\n",
" out = np.zeros(grid.blast_os.shape[0])\n",
" for e in range(out.size):\n",
" bp = grid.blast_os[e] - sensor.origin\n",
" fwd, right, up = grid.frames[e]\n",
" el = grid.el_off[e][:, None]\n",
" az = grid.az_off[e][None, :]\n",
" valid = grid.valid_el[e][:, None] & grid.valid_az[e][None, :]\n",
" cos_el = np.cos(el)\n",
" dirs = (fwd * (cos_el * np.cos(az))[..., None]\n",
" + right * (cos_el * np.sin(az))[..., None]\n",
" + up * np.sin(el)[..., None])\n",
" # ray-sphere quadratic: |bp + s*dir|^2 = r^2 -> s^2 + 2bs + c = 0\n",
" b = dirs @ bp\n",
" c = bp @ bp - r * r\n",
" disc = b * b - c\n",
" sq = np.sqrt(np.maximum(disc, 0.0))\n",
" s_near, s_far = -b - sq, -b + sq\n",
" use_near = (s_near >= 0) & (s_near <= max_len)\n",
" use_far = (s_far >= 0) & (s_far <= max_len)\n",
" s = np.where(use_near, s_near, s_far)\n",
" hit = (disc >= 0) & (use_near | use_far) & valid\n",
" impulse = np.where(hit, p0 / np.maximum(s, 1e-6) ** 2 * k_const, 0.0)\n",
" out[e] = np.sum(impulse * cos_el * res * res)\n",
" return out\n",
"\n",
"\n",
"reference = no_armor_reference(cfg, sensor, grid, k_const)\n",
"diff = np.abs(results[\"no_armor\"] - reference)\n",
"scale = reference.max()\n",
"print(f\"max |warp - numpy| = {diff.max():.3e} \"\n",
" f\"({diff.max() / scale:.2e} of the data scale)\")\n",
"assert diff.max() <= 1e-2 * scale\n",
"print(\"Cross-check passed: Warp's no_armor column matches the independent \"\n",
" \"NumPy reference.\")\n",
"\n",
"# Aside: the total captured solid angle should be close to the spherical cap\n",
"# 2*pi*(1 - cos(theta)). Implied here from the impulse at the cube center.\n",
"cap = 2.0 * np.pi * (1.0 - np.cos(theta0))\n",
"implied = results[\"no_armor\"][e0] / (p0 / d0 ** 2 * k_const)\n",
"print(f\"solid angle at the cube center: cap formula {cap:.3e} sr, \"\n",
" f\"implied by the simulation {implied:.3e} sr\")"
]
},
{
"id": "diffraction-md",
"cell_type": "markdown",
"metadata": {},
"source": [
"## 8. Optional: edge diffraction (filling the shadow)\n",
"\n",
"Everything so far is pure ray tracing, which casts *sharp* shadows — a blast the\n",
"armor occludes contributes exactly zero. But the pulse is acoustically **large**:\n",
"its positive phase spans $c\\,t_d \\approx 3.4$ m, roughly ten times the helmet, so\n",
"real blast energy diffracts *around* the shell into those shadows. Setting\n",
"`diffraction.enabled: true` adds **first-order edge diffraction** so the geometric\n",
"shadows fill in physically.\n",
"\n",
"The method is **Biot–Tolstoy–Medwin (BTM)**, not the more familiar GTD/UTD.\n",
"Impulse is the zero-frequency (DC) component of pressure, and the high-frequency\n",
"GTD/UTD coefficient *diverges* at DC (it scales $1/\\sqrt{k}$); BTM's time-domain\n",
"response is integrable, so its DC edge gain is finite — the right tool for an\n",
"impulse metric. The diffracted impulse folds into the same accumulator as the\n",
"direct rays, matched so geometric-optics + diffraction is continuous across the\n",
"shadow boundary. The coefficient is validated against the exact rigid-wedge\n",
"Green's function in `tests/test_diffraction.py`.\n",
"\n",
"Because the MICH mesh is a smooth watertight shell with **no rim edge**, the\n",
"helmet's diffractor is sourced procedurally — the shell's intersection with a\n",
"horizontal (constant-Z) plane. Treat diffraction-on numbers as a *comparative\n",
"robustness bound* (does the armor ranking survive the model's biggest omission?),\n",
"not absolute dosimetry.\n",
"\n",
"On this small **frontal** demo grid the magnitude change is mild — the takeaway\n",
"below is that the per-point armor ranking shifts for almost every blast point,\n",
"so the geometric-only ordering is not robust to diffraction. The stronger\n",
"shadow-fill and amplification show up in a full centered-cube sweep (see the\n",
"README's *Edge diffraction* section). The diffraction term is signed, so the\n",
"engine clamps the (physical, nonnegative) impulse at zero and logs how many\n",
"transition cells it clamped."
]
},
{
"id": "diffraction-code",
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"import diffraction\n",
"\n",
"# The helmet is a smooth closed shell, so the diffracting rim is the shell's\n",
"# intersection with a horizontal plane (a clean closed curve lying on the shell).\n",
"rim = diffraction.rim_from_plane(helmet, normal=(0, 0, 1), offset_m=-0.09,\n",
" seg_len_m=0.02, nu=0.5)\n",
"print(f\"procedural helmet rim: {rim.apex.shape[0]} segments\")\n",
"\n",
"cfg_d = config_mod.load_config(\"config.yaml\", DEMO_OVERRIDES + [\n",
" \"diffraction.enabled=true\", \"diffraction.segment_length_m=0.02\"])\n",
"res_d = warp_core.simulate(\n",
" cfg_d, helmet, vest, back, sensor, grid, k_const, device,\n",
" diffraction.build_sensor_coupling(grid),\n",
" diffraction.build_edge_sets(cfg_d, helmet, vest, back))\n",
"\n",
"print(f\"{'config':>20} {'shadow cells':>13} {'filled':>7} {'lit mean':>10}\")\n",
"for name in CONFIG_ORDER:\n",
" go, di = results[name], res_d[name]\n",
" shadow = go == 0.0\n",
" filled = int(((di > 0.0) & shadow).sum())\n",
" lit = go > 0.0\n",
" pct = 100.0 * np.mean((di[lit] - go[lit]) / go[lit]) if lit.any() else 0.0\n",
" print(f\"{name:>20} {int(shadow.sum()):>13} {filled:>7} {pct:>+9.1f}%\")\n",
"\n",
"# does enabling diffraction change the per-point armor ranking?\n",
"rank_go = np.argsort([-results[n] for n in CONFIG_ORDER], axis=0)\n",
"rank_di = np.argsort([-res_d[n] for n in CONFIG_ORDER], axis=0)\n",
"changed = int(np.any(rank_go != rank_di, axis=0).sum())\n",
"print(f\"\\nblast points whose armor ranking changed: {changed}/{grid.blast_os.shape[0]}\")"
]
},
{
"id": "next-md",
"cell_type": "markdown",
"metadata": {},
"source": [
"## 9. Where to go from here\n",
"\n",
"- **Real sweeps** run through the CLI, which writes a CSV\n",
" (`bp_x, bp_y, bp_z, full_armor, front_plate_helmet, helmet_only, vest_only,\n",
" no_armor`):\n",
"\n",
" ```sh\n",
" uv run simulate_blasts.py --config config.yaml \\\n",
" --set blast.cube_segments=15 --set rays.resolution_arcmin=2\n",
" ```\n",
"\n",
" Any override you experimented with above works there unchanged. Result CSVs\n",
" from earlier runs are committed under `data/` (pre-backplate files have four\n",
" armor columns; their `full_armor` is today's `front_plate_helmet`).\n",
"- **Where the plates come alive**: this notebook's demo casts rays only within\n",
" each sensor cone, which structurally cannot see paths that *start* aimed at\n",
" armor and reflect onto the sensor. Re-run with a centered cube and the\n",
" geometry-wide fan and the plates turn into small reflection amplifiers:\n",
"\n",
" ```sh\n",
" uv run simulate_blasts.py --config config.yaml \\\n",
" --set blast.cube_center_m=[0,0,0] --set blast.cube_extents_m=24 \\\n",
" --set rays.fan_mode=geometry --set rays.resolution_arcmin=1\n",
" ```\n",
"- **Statistics**: `analysis.R` runs an ANOVA and Holm-adjusted paired t-tests\n",
" across the armor configurations; `visualize.R` has ggplot/rgl snippets\n",
" (source `analysis.R` first).\n",
"- **Tests**: `uv run pytest` checks the physical invariants and performs the\n",
" same NumPy cross-check as §7.\n",
"- **Performance**: see the README's *Precision and performance notes* — Warp's\n",
" mesh queries are float32 on a software BVH; the next big lever at very fine\n",
" resolutions is an RT-core path (Mitsuba 3 / Dr.Jit)."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"name": "python"
}
},
"nbformat": 4,
"nbformat_minor": 5
}