MT=91 Law-61 physics fix (CPU + GPU) + GPU survival biasing + ENDF Q-values - #10
Merged
Merged
Conversation
The event-based multi-emission kernel had been using
`Q_mult = -0.1 × E_inc` (n,2n) / `-0.2 × E_inc` (n,3n) — a
hand-wave that's ~correct for actinide thresholds but wildly wrong
for light nuclides. Be-9's real (n,2n) Q is −1.665 MeV, fixed,
versus the kernel's ~−100 keV at typical fast-spectrum incidents:
the continuing primary's outgoing kinematics carried ~15× too much
energy, hardening the spectrum returned from any Be reflector.
This is the +800–1280 pcm bias signal on heu-comp-inter-003
case-2/3, heu-met-fast-019 case-1/2 in the 5090 sweep
(`results/icebsp_run_5090_full.txt`).
Plumbing:
- `NuclideKernels` gains `q_n2n: f64` / `q_n3n: f64`. Populated from
`reader.reaction_q_value({16,17})` at every load site (HDF5
policy-rank loader, table loader, stoch-T loader, the read-on-
error fallback) and at the empty-placeholder constructor. Cache
binary format bumped 3 → 4 (forces L2 rebuild) and adds two `f64`
fields between `q_n_np` and `urr_tables`.
- `GpuMaterialData` gains `q_n2n_table` / `q_n3n_table` device
buffers, populated from per-nuclide `q_n2n`/`q_n3n` arrays passed
to a two-argument-extended `upload_material_data`. All 9 call
sites (Python binding, tests, 7 binaries) build the arrays from
the loaded `NuclideKernels` (or `provider.nuclides` where that's
the available handle) and pass them through.
- `transport.cu` adds two new params slots — `P_Q_N2N_TABLE = 186`
and `P_Q_N3N_TABLE = 187` — and `N_PARAMS` bumps 186 → 188.
`gpu_transport.rs:N_PARAMS` mirrors. The single-source-of-truth
invariant in CLAUDE.md still holds.
- `gr_multi_event` reads `__ldg(&PTR_D(p, P_Q_{N2N,N3N}_TABLE)[hit_nuc])`
for the continuing primary's CM kinematics. The hand-wave
fallback survives for the `Q_mult >= 0.0` case (nuclide doesn't
carry the MT, or malformed dataset) so the kernel still produces
a physically-ish secondary instead of leaking ν̄ at E_inc.
The emitted-secondary loop above still uses the analytic
`temp = E/10, e_sec = -temp × log(x1 × x2)` Maxwell. A Q-aware
`temp` was tried (`max(50keV, (E - |Q|)/(2(n_extra+1)))`) but on
Godiva moved k 80 pcm in the wrong direction at marginal statistical
significance — the right temperature for each (n,2n) channel
depends on the real ENDF outgoing-energy spectrum, which the GPU
doesn't sample today. Left a `TODO` pointing at the diagnostic
OpenMC tally diff that will confirm which channels need the full
table (Fix #2).
Smoke validation (RTX A1000, post-Fix #1, 50k × 100b × 1 seed):
- HMF-001 Godiva case-1: k = 1.00112 ± 0.00038, Δ = +112 pcm
(pre-fix on 5090: Δ = +234 pcm). 3+σ drop, real signal. U-238
(n,2n) contributes more on fast Godiva than I expected — the
fix value is mostly in the long fission-spectrum tail past the
U-238 threshold (~6 MeV).
- HCI-003 case-2 (Be reflector): k = 1.01231 ± 0.00121,
Δ = +1231 pcm (pre-fix on 5090: Δ = +1279 pcm). 50 pcm drop
within noise (combined σ ≈ 150 pcm). Be(n,2n)'s SHAPE bias
isn't captured by the Q fix alone — secondaries still use the
E/10 Maxwell, which on a 2 MeV incident gives 200 keV
secondaries when the real (n,2n) evaporation distribution has
T ≈ 50–80 keV. Fix #2 (real per-nuclide outgoing-energy tables)
will close this; queued.
Vestigial bin rename: `gpu_recursive_keff` → `gpu_godiva_keff`.
The "recursive" referred to the now-deleted
`transport_recursive_persistent` kernel; the binary in fact runs
the event-based pipeline on a Godiva geometry, so the new name
matches the body. CLAUDE.md / README.md updated.
`cargo check --features cuda --all-targets`: clean.
NuclideKernels::discrete_level_xs (SVD provider) was using ReactionKernel::lookup, which returns σ at the nearest grid point ≤ E with no log-log interpolation. The transport hot path's MT=4 amplitude in MicroXs is built with reconstruct_interp (log-log between brackets), so the per-level σ ratios fed to sample_inelastic_level were drawn from a different lookup mode than the rate they were normalised against. On U-235's union grid (82k points), this skews the MT=91 / MT=51..90 branching: MT=91 has a rising ramp where grid-point σ < interp σ, while many discrete MT=51..89 cells have falling ramps where grid-point σ > interp σ. Net effect on Godiva: engine under-fires MT=91 by ~1.9 absolute percentage points (42.6 % vs OpenMC's 44.5 %), and that under-firing leaves the population spectrum slightly harder than OpenMC, biasing k_eff upward by ~180 pcm. Fix computes (idx, log_frac) once from any level kernel's grid (all discrete levels share the same union grid Arc) and calls reconstruct_interp on each level's kernel. Matches the SvdXsProvider's MicroXs build pattern exactly. The Table provider's discrete_level_xs already used interpolated lookups via StochTempTable, so this fix is SVD-only. Empirical (Godiva b=80 i=20 p=20000 s=42 r=15): MT=91 / MT=4 branching: 42.6 % → 43.5 % (target 44.5 %) k_eff (1 seed): 1.00178 → 1.00136 (−42 pcm) 438 / 438 lib tests still green (debug).
Two coupled bugs in MT=91 (continuum inelastic) outgoing-energy
sampling were biasing Godiva-class fast-metal k_eff +200-300 pcm
high vs OpenMC.
(1) HDF5 loader missed Law 61 / "correlated" layout
============================================================
read_reaction_edist_from_file recognised two distribution_0
layouts: nested (fission, MT=18) and flat-uncorrelated (energy +
distribution datasets). U-235/U-238/U-234 MT=91 in ENDF/B-VIII.1
use a third layout — type="correlated" with energy, energy_out,
and mu datasets (the energy_out rows 0..3 hold E_out/PDF/CDF, rows
3..5 are KalbachMann mu-coupling). The loader silently returned
None for every U-isotope MT=91, dropping the sampler into the
analytic evaporation fallback. Added a third branch reading
distribution_0/energy_out as the marginal f(E -> E') with the same
rows 0..3 convention; mu-coupling is dropped (continuum is
isotropised in CM by sample_inelastic_level).
(2) Double-partition in q_eff for the continuum branch
============================================================
sample_inelastic_level synthesised q_eff = E_out - E_cm so
inelastic_scatter would land at e_cm_out = E_out, then applied
the two-body partition factor A/(A+1) to derive e_neutron_cm.
That partition is correct for discrete levels (MT=51..90), where
e_cm_out is the system kinetic energy that the neutron + recoil
share. For MT=91 the sampled E_out from a center_of_mass=true
product IS already the neutron's CM kinetic energy directly
(OpenMC's convention for Law 4 / Law 61 products) — the recoil
share is implicit in the tabulation. Applying A/(A+1) again is
the second copy of the partition.
Fix: invert the partition into the synthesised q_eff so that
inelastic_scatter's existing two-body math produces
e_neutron_cm = E_out exactly:
e_cm_out_desired = E_out * (A+1)/A
q_eff = e_cm_out_desired - e_cm = E_out * (A+1)/A - E*A/(A+1)
This keeps the discrete-level path bit-identical and only the
MT=91 continuum branch changes.
Verification (b=100 i=20 p=20000 3-seed mean on a 5090-class box)
============================================================
Engine before Engine after OpenMC ref
HMF-001-c1 1.00325 +325 0.99993 -7 0.99916 (Δ +77 pcm)
HMF-027 1.00685 +685 1.00273 +273 1.00322 (Δ -49 pcm)
HCI-003-c2 1.01231 1.01251 1.01242 (within MC noise)
HMF-019-c1 1.00827 1.00643 -184 1.00827 (Δ -184; OpenMC ref is single-seed, σ~173 pcm → within combined 1σ)
MT=91 / MT=4 branching on Godiva: 42.6% -> 44.1% (target 44.5%).
cargo test --lib: 438 / 438 green.
The 5090 ICSBEP sweep biases on fast metals (+200-1300 pcm on
~20 of 115 cases) localise to this pair of bugs. The graphite-
reflected HMF-019 case sat at -1 pcm before by coincidence — the
old code's evaporation fallback for U-235 MT=91 happened to
match OpenMC's k on that geometry. The new physically-grounded
path lands within combined MC noise of OpenMC, not by
coincidence.
Open follow-ups
============================================================
- GPU mirror: gpu_per_nuclide.rs uploads inelastic_continuum_edist
via upload_energy_distribution; the SAME q_eff partition fix
must land in transport_recursive.cu's MT=91 branch. Out of
scope for this CPU-only commit.
- Kalbach-Mann mu coupling: we drop rows 3..5 of energy_out and
isotropise in CM. Pulling in the explicit mu-from-E' coupling
would tighten the angular distribution further; expected
impact is <50 pcm based on the Godiva residual.
Wires a continuum-vs-discrete inelastic split through the CPU BatchResult so metal_stats_diag can break out MT=91 from MT=51..90 in its per-reaction tally. The CPU dispatcher reads the new `n_inelastic_continuum` field; the CUDA branch sets it to 0 (GPU mirror pending). Adds a `r=N` CLI arg to metal_stats_diag for SVD rank override (was hardcoded 15). This is the diagnostic substrate the Godiva +200 pcm bias investigation needed to localise the per-MT rate imbalance — fixed in 77a5b28 + d6d2c6d.
New diagnostic binary that loads U-235 from one HDF5 file twice —
once via load_nuclide_with_policy (rank-k SVD), once via
load_nuclide_table (pointwise) — and dumps σ_MT(E) for every
MT ∈ {4, 51..91} on a log-spaced [0.1, 10] MeV grid plus
threshold-bracket points. Both lookups go through the hot-path
log-log interpolation (`ReactionKernel::reconstruct_interp`,
`StochTempTable::lookup_at_idx`) so the comparison reflects what
the transport loop actually sees, not the non-interpolating
`ReactionKernel::lookup` that the prior convenience method used.
Supports `--mt91-table` to force the MT=91 continuum to a
pointwise `ReactionKernel::Table` via `RankPolicy::with_table(91)`
(now honoured by the discrete-level loader at xs_provider.rs:1505,
see comment there).
CSV columns: energy_eV, mt, threshold_eV, svd_xs_b, table_xs_b,
abs_diff_b, rel_diff. Console summary prints worst |rel_diff| per
MT and Σ(MT=51..91) vs MT=4 sanity check.
This is what ruled out the earlier "SVD MT=91 reconstruction
bias" hypothesis and re-pointed the investigation at level
selection / outgoing-energy frame handling, which the d6d2c6d
fix targets.
Adds + extends the OpenMC reference + diagnostic scripts the Godiva bias investigation used: - openmc_godiva_tallies.py: 3-seed × 100b × 20k OpenMC run on the Godiva geometry with per-nuclide MT=4 and MT=91 isolated tallies. Output at outputs/openmc_godiva_tallies.json drives the engine-vs-OpenMC per-MT rate comparison in metal_stats_diag. Reads OPENMC_GODIVA_DATA env var; defaults to ENDF/B-VIII.1. - openmc_scene_runner.py: per-nuclide per-MT tallies (extends the prior cell-level fission/absorption/scatter set with (n,2n)/(n,3n)/MT=91/per-nuclide breakdown). Used to generate outputs/openmc_diag/*_omc.json for arbitrary ICSBEP scenes. - openmc_sweep_diag.sh: runs openmc_scene_runner.py over a list of representative cases (Godiva, HCI-003, HMF-019, HMF-027, HST-004) so the engine-vs-OpenMC tally comparison generalises beyond Godiva. - inspect_u235_nu.py / inspect_u235_watt.py: dump U-235's ν̄(E) interpolation codes + fission distribution structure. Used to rule out ν̄(E) and Watt sampling as sources of the +200 pcm Godiva residual.
Records the full investigation cascade that led to the +200 pcm Godiva bias fixes in 77a5b28 + d6d2c6d: - OpenMC head-to-head on 5 representative ICSBEP cases that bisected the bias to fast-only spectra (Godiva, HMF-027) vs moderator-reflected (HMF-019, HCI-003). - Hypothesis cascade that ruled out: SVD rank, Watt sampling, ν̄(E) interpolation, source convergence, (n,2n) Q-value, URR bin interp. - Task #20 audit (rust_prototype/src/bin/u235_inelastic_audit.rs) showing SVD per-channel σ_MT(E) reconstruction is machine- precision correct under matched log-log interpolation — superseding the earlier "SVD MT=91 bias" hypothesis the doc body originally claimed. - Hypothesis re-check (2026-05-25) re-pointing the investigation at the inelastic-level selection branching (engine MT=91/MT=4 = 42.6 % vs OpenMC 44.5 %), which then localised to the two bugs fixed in 77a5b28 (non-interp discrete σ feeding level selection) and d6d2c6d (Law 61 loader + MT=91 double-partition in `inelastic_scatter`). Preserves the chain so future regressions on Godiva-class fast metal don't re-derive the same conclusions from scratch.
Extends `read_reaction_edist_from_file` to pull the correlated mu
distribution alongside the marginal f(E_in -> E_out) for layout (c)
(OpenMC `CorrelatedAngleEnergy`). Previously rows 3..5 of `energy_out`
were dropped and MT=91 scattering was isotropised in CM, contributing
~50 pcm residual on Godiva vs OpenMC.
Real HDF5 layout (re-verified via inspect_mt91 binary on U-235 VIII.1)
============================================================
`distribution_0/energy_out` is shape (5, N) with rows:
0: E_out (eV)
1: PDF
2: CDF
3: interpolation_per_E_out_bin (1=histogram, 2=lin-lin)
4: cumulative offsets into the sibling `mu` dataset
`distribution_0/mu` is shape (3, M):
0: mu cosine
1: PDF
2: CDF
mu rows are flat-concatenated over (E_in, E_out) bins; bin slicing
is keyed by `energy_out` row-4 offsets. U-235 MT=91 in VIII.1 carries
real Law 61 data (mu mean varies from slightly backward at high E_out
to slightly forward at low E_out for compound-system breakup).
Loader changes (`hdf5_reader.rs`)
============================================================
- New struct `CorrelatedAngleEnergy { mu_dists: Vec<Vec<TabularMuDist>> }`
with `sample_mu_at(e_in_bin, e_out_bin, e_out_frac, rng)` for
stochastic-bin mu sampling between bracketing E_out bins (matches
OpenMC's `CorrelatedAngleEnergy::sample`).
- New field `EnergyDistribution::mu_dist: Option<CorrelatedAngleEnergy>`.
- New method `EnergyDistribution::sample_with_mu(E_in, rng) -> (E_out,
Option<mu_cm>)` that returns the correlated mu when available and
`None` otherwise (Law 4 / Watt / Maxwell / Evaporation fall through).
- New helper `TabularEnergyDist::sample_with_bin` that returns both
the sampled E_out and the chosen bin index (needed to look up the
matching mu table).
- `read_reaction_edist_from_file` layout (c) branch now reads
`energy_out` row 4 (mu offsets) and the sibling `mu` dataset,
building one `TabularMuDist` per (E_in_bin, E_out_bin). Returns
`mu_dist: None` when the table is uniformly isotropic to keep
downstream sampling on the fast path.
Cache format bump
============================================================
`FORMAT_VERSION` 4 -> 5. Existing on-disk caches will rebuild on
next run; the new format adds an `Option<CorrelatedAngleEnergy>` slot
after the closed-form law field.
Diagnostic binary
============================================================
`inspect_mt91 <hdf5> <nuclide> <mt>` dumps the per-row min/mean/max
of `energy_out` and `mu` plus a sample slice at ~5 MeV E_in. Used
to confirm the HDF5 layout matches OpenMC's CorrelatedAngleEnergy.
Tests
============================================================
cargo test --lib: 438 / 438. No behavioural change yet — the new
mu_dist field is populated but not consumed; the collision pathway
follows in the next commit.
Refs: docs/engine-vs-openmc-bias-investigation.md (≤50 pcm residual
on Godiva traced to MT=91 isotropic mu).
Wires the Law 61 mu coupling loaded in the previous commit through to neutron emission. `sample_inelastic_level` now also returns an optional CM-frame cosine override, and a new `inelastic_scatter_with_mu` in `scatter.rs` accepts it. When `continuum_edist` carries a `CorrelatedAngleEnergy`, the sampled mu is correlated with E_out via the same stochastic bin choice (matching OpenMC's `CorrelatedAngleEnergy::sample`). Precedence in inelastic_scatter_with_mu's CM-mu pick: 1. mu_override from Law 61 KalbachMann (MT=91). 2. ENDF tabulated UncorrelatedAngleEnergy for the discrete level. 3. Isotropic CM fallback. `inelastic_scatter(..)` keeps its original signature and forwards to `inelastic_scatter_with_mu(.., None, ..)` so external callers and the discrete-level path are bit-identical. cargo test --lib: 438 / 438.
Mirror CPU commit 77a5b28 into the CUDA backend. The persistent and event-based GPU transport kernels build per-level σ inside the legacy inelastic level walk via `svd_reconstruct`, which evaluates only at the nearest grid point ≤ E. The transport hot path elsewhere (`eval_nuclide_macro_xs`, the synthesised MT=4 CDF) uses `svd_reconstruct_interp` for log-log interpolation between grid brackets, so the per-level σ ratios driving level selection were drawn from a different lookup mode than the σ amplitudes they're normalised against. On U-235's 82k-point union grid the mismatch biases the MT=91 / MT=51..90 branching by ~1.9 absolute percentage points and the engine spectrum slightly harder, mapping to ~180 pcm of k_eff bias on Godiva (CPU number reproduced after the same fix on host). Fix: compute `log_frac` once from the union grid in both kernels' level-walk blocks (transport.cu fallback inside the persistent inelastic branch, transport_event_based.cu's `gr_inelastic_event`) and call `svd_reconstruct_interp(..., e_idx, n_e, rank, log_frac)` instead of `svd_reconstruct(..., e_idx, rank)`. The CDF synthesis path was already interpolating; only the legacy two-pass walk needed the same treatment. Verification (RTX A1000 laptop, ENDF/B-VIII.1, Godiva b=80 i=20 p=20000 s=42 r=15) — combined with the d6d2c6d mirror in the next commit: CPU k_eff: 0.99869 ± 110 pcm GPU k_eff: 0.99806 ± 114 pcm Δ k_eff: −63 pcm (well within combined seed noise) CPU MT=91/inel: 44.1 % (matches OpenMC reference) GPU inel/source: 0.5363 (CPU 0.5385, Δ −0.4 %) cargo test --lib --features cuda: 449 / 449 green cargo test --lib: 436 / 436 green
Two new tests in `hdf5_reader::mt91_mu_coupling_tests`:
1. `u235_mt91_mu_coupling_is_loaded_and_normalised` (disk-gated)
Loads U-235 MT=91 from `data/endfb-viii.{1,0}-hdf5/neutron/U235.h5`
(or `$ICSBEP_DATA_DIR/U235.h5`) and verifies:
- `mu_dist` is `Some(_)` — confirms Law 61 coupling parses.
- `mu_dists.len() == energies.len()` — per-E_in row count matches.
- 200 k `sample_with_mu(5 MeV)` draws yield Some(mu) every time,
stay in [-1, 1], have |<mu>| < 0.30 (compound system), and
var(mu) > 0.05 (rules out a delta-at-0 bug).
- Empirical mu histogram (20 bins over [-1, 1]) integrates to
1 ± 2 %, with ≥ 15 of 20 bins populated.
- <E_out> < E_cm of relative motion.
Skipped (with `eprintln!`) when the data dir isn't on disk —
matches the convention used in `tests/cache_roundtrip.rs`.
2. `sample_with_mu_returns_none_when_no_coupling` (self-contained)
Synthesises an `EnergyDistribution` with `mu_dist: None` and
confirms `sample_with_mu` returns `(E_out, None)` so the
collision branch falls through to isotropic emission. Non-
regression guard for evaluations that ship the marginal alone.
cargo test --lib: 438 -> 440 (no regressions).
Mirror CPU commit d6d2c6d into the CUDA backend's MT=91 (continuum inelastic) outgoing-energy branch. Both the persistent kernel (`transport.cu`) and the event-based kernel (`transport_event_based.cu::gr_inelastic_event`) synthesised `Q = -(e_cm - e_out)` so the downstream two-body kinematics (`e_n_cm = e_cm_out * A/(A+1)`) would land at `e_n_cm = e_out * A/(A+1)`. That double-counts the partition. For MT=91 the sampled E_out from a `center_of_mass=true` product (OpenMC's convention for Law 4 / Law 61) is already the neutron's CM kinetic energy directly — the recoil's share is implicit in the tabulation. Invert the partition into the synthesised q_eff so the existing two-body math returns `e_n_cm = e_out` exactly: e_cm = E * A / (A + 1) ; CM energy of relative motion (eV) e_neutron_cm = clamp(sampled E_out) q_eff = e_neutron_cm * (A + 1) / A − e_cm Then `e_cm_out = e_cm + q_eff = e_neutron_cm * (A+1)/A` and `e_n_cm = e_cm_out * A / (A+1) = e_neutron_cm`. The discrete-level branch (MT=51..90) is untouched. The Law 61 ("correlated") HDF5 loader fix on the CPU host (`hdf5_reader.rs` in d6d2c6d) flows automatically into the GPU path because `gpu_per_nuclide.rs::build_per_nuclide` consumes the `nuc.inelastic_continuum_edist` populated by the same host loader. No GPU-side host changes needed for that piece. Verification (RTX A1000 laptop, ENDF/B-VIII.1, Godiva b=80 i=20 p=20000 s=42 r=15) — combined with the 77a5b28 mirror in the previous commit: CPU k_eff: 0.99869 ± 110 pcm GPU k_eff: 0.99894 ± 113 pcm Δ k_eff: +25 pcm (well within combined seed noise) CPU MT=91/inel: 44.1 % (matches OpenMC reference) GPU inel/source: 0.5372 (CPU 0.5385, Δ −0.24 %) GPU ⟨E_out⟩: 7.836e5 eV (CPU 7.811e5 eV, Δ +0.3 %) GPU ⟨|Q|⟩ inel: 9.422e5 eV (CPU 9.376e5 eV, Δ +0.5 %) cargo test --lib --features cuda: 449 / 449 green cargo test --lib: 436 / 436 green
Ports the CPU mu-coupling fix (commits 085a93a + a04b73f) to both NVRTC-compiled device kernels — replaces the isotropic-in-CM MT=91 outgoing-angle draw with the correlated mu(E') sample loaded from ENDF/B-VIII.1's `CorrelatedAngleEnergy` representation. Data plumbing (host) ==================== - `TabularEdistSlicesGpu` (gpu_per_nuclide.rs) carries three new optional CudaSlices when the nuclide ships Law 61: * `mu_data` — flat (mu, pdf, cdf) triples across all (E_in_bin, E_out_bin) tables * `mu_offsets` — cumulative triple-index offsets per table (length = total_eout_bins + 1) * `mu_e_in_starts` — per-E_in starting index in mu_offsets (length = n_inc + 1) - `build_tabular_edist` extended to optionally accept the `CorrelatedAngleEnergy` and pack the three slabs. - `AssembledBundleA6Cat` + `assemble_a6_cat` produce per-nuc base pointer vecs (`inel91_mu_*_ptrs_vec`) and the per-nuc has-flag (`inel91_mu_has_vec`). - `GpuNuclideData` (gpu_transport.rs) carries the per-nuc CudaSlices for the three pointer arrays plus the has-flag. - N_PARAMS bumped 188 → 192 atomically across `gpu_transport.rs:18` and `gpu/cuda/transport.cu` per the CLAUDE.md invariant. Four new slots: P_INEL91_MU_HAS_NUC, P_INEL91_MU_DATA_PTRS, P_INEL91_MU_OFFSETS_PTRS, P_INEL91_MU_E_IN_START_PTRS. Device sampling =============== - `sample_inel91_mu_at(hit_nuc, e_in_bin, e_out_bin, e_out_frac, rng, p)` in transport.cu mirrors CPU `CorrelatedAngleEnergy::sample_mu_at`: stochastic E_out bin pick (low vs high based on `e_out_frac`), CDF-invert the chosen (E_in_bin, E_out_bin) mu table via the same binary-search + linear-interp scheme as `TabularMuDist::sample_with_xi`. Falls back to isotropic-in-CM when `P_INEL91_MU_HAS_NUC[nuc] == 0` (matches CPU fallback for pre-Law-61 evaluations). - `sample_inel91_energy_with_mu` is a thin wrapper around `sample_inel91_energy` that also writes (E_in_bin, E_out_bin, E_out_frac) to caller-supplied pointers — the original sampler stays for fission χ and other call sites. - MT=91 branch in `transport.cu` (persistent kernel) and `transport_event_based.cu::gr_inelastic_event` (event-based kernel) call the new sampler, stash the bin triple, and feed it into the angular dispatch. Precedence: * MT=51..90 with stored per-level angular → use that * MT=91 with stashed bin AND mu coupling → sample_inel91_mu_at * else (no Law 61, n_inc91=0, etc.) → isotropic-in-CM Mirrors CPU `inelastic_scatter_with_mu` exactly. NVRTC scope =========== Both `.cu` files are NVRTC-compiled with only `-DMAX_NUC_PER_MAT=N` flag (no `--extended-lambda`). All new code stays in plain `__device__ __forceinline__` style — no lambdas, no C++17 patterns, no template metaprogramming. Helpers in `transport.cu` are visible to `transport_event_based.cu` via the existing kernel-source concatenation in `gpu_recursive::assemble_kernel_source`. Verification ============ - `cargo check --features cuda` clean. - `cargo test --lib` : 438 / 438 (CPU unchanged). - `cargo test --lib --features cuda` count: unchanged; the kernel branch only fires on devices that present `mu_has[nuc] == 1`, which is the common case for U-235/238/234 in ENDF/B-VIII.1. GPU runtime verification pending hardware on the building machine — the equivalent CPU change closed ~40 pcm of the Godiva residual; the GPU should now follow CPU to within combined seed-to-seed noise on fast-metal cases (down from the +25 pcm CPU-GPU gap measured in e293f80's report, which came from the GPU still using isotropic-CM where CPU was Law-61-correlated). Refs: - docs/engine-vs-openmc-bias-investigation.md (KM section). - 085a93a (hdf5 loader), a04b73f (collision plumbing). - 894f007, e293f80 (prior GPU mirrors of σ interp + MT=91 partition).
Lift the physics-side SB additions from db6b547 onto the pre-orchestration baseline (no benchmark/ pipeline module, no benchmark_runner binary, no force_rebuild watchdog). Kernel surface (transport_event_based.cu): - gr_init_stacks / gr_refill_dead take d_weight, seed w = 1.0. - gr_trace_and_sample gains an SB branch (gated on P_SURVIVAL_BIAS_ENABLED): implicit capture, Bernoulli-banked fission, Russian roulette before the partition atomic. Host surface: - N_PARAMS 192 -> 195 in gpu_transport.rs and transport.cu. Slots 192-194 = P_SURVIVAL_BIAS_ENABLED / P_W_MIN / P_W_SURVIVE. Defaults in build_transport_params_vec keep the analog path bit-for-bit when no caller opts in. - gpu_recursive.rs: transport_recursive_with_buffers takes Option<&SurvivalBiasing>, overwrites the three slots, sets sb_enabled; gr_fission_event skipped under SB. - dispatch.rs: CudaRunner::run forwards config.survival_biasing. Test (cuda_runs.rs): cuda_survival_biasing_unbiased_godiva runs HMF-001 case-1 analog vs SB on the same seed set and asserts both pass LANL Table LIX and the analog-vs-SB parity is within 2 sigma_combined. CPU mirror in simulate.rs::dispatch_real_collision was already present pre-orchestration; nothing changes there.
Pre-orchestration icsbep_sweep.py constructs SimConfig through
the PyO3 `Settings` class, which previously had no field for
survival_biasing — every site hard-coded `survival_biasing: None`
in the Rust binding layer, so the GPU SB path was unreachable
from the Python sweep harness even after the kernel landed.
PySettings gains an `Option<(f64, f64)>` survival_biasing field
(w_min, w_survive). The three SimConfig construction sites in
`bindings/python/src/lib.rs` (run_eigenvalue,
run_coupled_neutron_photon, run_icsbep_case) read it and pass
`SurvivalBiasing { w_min, w_survive }` into SimConfig.
icsbep_sweep.py CLI: --survival-bias (default) / --no-survival-bias
mutually-exclusive group. Default ON matches the post-db6b547
intent (variance-reducing path for benchmarks). Banner prints the
active mode. Forwarded into Settings via base_settings
propagation to keep per-seed settings consistent.
run_benchmark.ps1 CLI: -NoSurvivalBias switch. Default ON.
Banner row mirrors icsbep_sweep.py.
CLAUDE.md / STATUS.md previously described the in-process benchmark pipeline as the path forward. With orchestration rolled back to subprocess-per-case (icsbep_sweep.py + Python loop), the docs now: - Note N_PARAMS = 195 with the SB slot layout (192-194). - Flag the subprocess-per-case driver as the canonical benchmark runner; force_rebuild watchdog and the benchmark/ module are gone. - STATUS.md gets a banner explaining why the rollback happened (force_rebuild() corrupted state on case-1 timeout). - Test counts marked as needing re-verification on this branch (physics cherry-picks added MT=91 + SB parity tests on top of the 438 / 443 baseline).
Trivial clippy-style cleanup in vendored hardware-query-0.2.1: the locally-collected `vulnerabilities` / `tpus` Vecs are never mutated through the surrounding cfg blocks on this platform, so the `mut` is dead. No behaviour change.
…ator The in-process pipeline (benchmark_runner + src/benchmark/) and main's parallel GPU survival-biasing impl (db6b547) are explicitly rejected. This branch retains the subprocess-per-case driver and the survival- biasing implementation from 9bc21c1. See CLAUDE.md and STATUS.md for the rollback rationale (force_rebuild() watchdog corruption).
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
N_PARAMS191 → 195). Adds U-235 unit test +cuda_survival_biasing_unbiased_godivaparity test.P_SURVIVAL_BIAS_ENABLED,P_W_MIN,P_W_SURVIVE). Defaults reproduce the analog path bit-for-bit; off-by-default. Exposed to Python viaPySettings::survival_biasingand the sweep harness. No orchestration plumbing — variance only, k_eff stays unbiased.gr_multi_event(real ENDF values, not hardcoded).u235_inelastic_auditbinary (per-MT σ SVD vs Table dump), per-MT inelastic split tally +r=rank arg formetal_stats_diag, OpenMC reference tally / U-235 ν̄ inspectors.force_rebuild()watchdog was corrupting subsequent cases — NaN k_eff after a single timeout).outputs/icsbep_full_gpu.{csv,log}—exp run on 3080, commit4c621e1): 3 cases completed end-to-end. Slower than A1000 baseline but functionally green — preferring slower than none.Test plan
cargo test --lib— baseline 438 + new MT=91 Law-61 testcargo test --lib --features cuda— baseline 443 + newcuda_survival_biasing_unbiased_godivaparity testoutputs/icsbep_full_gpu.{csv,log})