# SPDX-License-Identifier: AGPL-3.0-or-later
# Copyright (C) 2025 SWGY, Inc
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
"""
Correctness checks for the Warp blast engine.
Specular-reflection model with the sensor (a brain) inside the helmet shell:
armor can amplify impulse (helmet reflecting rays onto the sensor) as well as
attenuate it, so there is NO simple ordering between armor configs. The
invariants asserted here hold regardless. The Warp mesh query is float32, so
comparisons use float32-scale tolerances.
"""
import math
import numpy as np
import pytest
import config as config_mod
import sampling
import simulate_blasts
# Small, fast grid; 30 arcmin resolves the sensor cone at these distances.
BASE_OVERRIDES = [
"blast.cube_segments=4",
"blast.cube_extents_m=6.0",
"blast.cube_center_m=[0,-4,0]",
"blast.min_dist_m=1.0",
"rays.resolution_arcmin=30",
"rays.max_bounces=2",
]
def _run(extra=None):
cfg = config_mod.load_config("config.yaml", BASE_OVERRIDES + (extra or []))
return cfg, simulate_blasts.run(cfg)
@pytest.fixture(scope="module")
def result():
return _run()
def _key(x, y, z):
return (round(x, 4), round(y, 4), round(z, 4))
def _direct_no_armor_reference(cfg):
"""Independent, geometry-free NumPy computation of the no_armor column.
Casts each cone ray from the blast, finds the first sensor-sphere hit along
the segment, and accumulates p0/eff**2 * K weighted by the differential
solid angle -- sharing no code with the Warp kernel.
"""
sensor = sampling.Sensor(np.zeros(3), cfg.sensor.radius_m)
grid = sampling.build_grid(cfg, sensor)
res = grid.res_rad
r = sensor.radius
p0 = cfg.physics.initial_peak_pressure_kpa
ml = cfg.rays.max_length_m
k_const = sampling.precompute_waveform_constant(
cfg.physics.positive_phase_duration_s, cfg.physics.decay_alpha)
out = {}
for e in range(grid.blast_os.shape[0]):
bp = grid.blast_os[e]
fwd, right, up = grid.frames[e]
el = grid.el_off[e][:, None]
az = grid.az_off[e][None, :]
vmask = grid.valid_el[e][:, None] & grid.valid_az[e][None, :]
cos_el = np.cos(el)
dirs = (fwd * (cos_el * np.cos(az))[..., None]
+ right * (cos_el * np.sin(az))[..., None]
+ up * np.sin(el)[..., None])
b = dirs @ bp
cq = bp @ bp - r * r
disc = b * b - cq
sq = np.sqrt(np.maximum(disc, 0.0))
s_lo, s_hi = -b - sq, -b + sq
use_lo = (s_lo >= 0) & (s_lo <= ml)
use_hi = (s_hi >= 0) & (s_hi <= ml)
s = np.where(use_lo, s_lo, s_hi)
hit = (disc >= 0) & (use_lo | use_hi)
eff = np.maximum(s, 1e-6)
imp = np.where(hit, p0 / (eff * eff) * k_const, 0.0)
out[_key(*bp)] = float(np.sum(imp * cos_el * res * res * vmask))
return out
def test_row_shape_matches_columns(result):
_, rows = result
assert len(rows[0]) == len(simulate_blasts.COLUMNS)
def test_nonnegative_and_finite(result):
_, rows = result
for row in rows:
for val in row[3:]:
assert np.isfinite(val) and val >= 0.0
def test_no_armor_matches_independent_reference(result):
cfg, rows = result
ref = _direct_no_armor_reference(cfg)
w = np.array([na for *_, na in rows])
r = np.array([ref[_key(x, y, z)] for x, y, z, *_ in rows])
scale = max(float(r.max()), 1e-12)
# float32 mesh-free direct path vs the float64 reference: bulk agreement to
# float32, with the occasional cone-edge ray flipping hit/miss (bounded by
# the data scale, not the tiny per-point value it lands on).
assert np.max(np.abs(w - r)) <= 1e-2 * scale, \
"no_armor differs from the independent reference beyond float32 scale"
def test_no_armor_positive_somewhere(result):
_, rows = result
assert max(na for *_, na in rows) > 0.0, \
"direct path never reached the sensor -- ray grid likely under-resolved"
def test_geometry_participates(result):
_, rows = result
ho = simulate_blasts.COLUMNS.index("helmet_only")
na = simulate_blasts.COLUMNS.index("no_armor")
# The helmet surrounds the sensor, so at least one point must differ from
# the geometry-free result (confirms the mesh query is active).
assert any(abs(r[ho] - r[na]) > 1e-9 for r in rows)
def test_min_distance_filter(result):
cfg, rows = result
for x, y, z, *_ in rows:
assert math.dist((x, y, z), (0, 0, 0)) >= cfg.blast.min_dist_m - 1e-9
def test_determinism(result):
# Warp accumulates via float64 atomic_add, whose order is not fixed, so allow
# tiny reassociation differences rather than requiring bit-equality.
_, rows1 = result
_, rows2 = _run()
a = np.array([r[3:] for r in rows1])
b = np.array([r[3:] for r in rows2])
assert np.allclose(a, b, rtol=1e-6, atol=1e-9), "same config must reproduce"