FW/CADIS-Ω angle-informed weight windows for fusion VR: extend random ray (not Sₙ), staged & gated

[jon-proximafusion]

Summary

Investigate adding angle-informed weight windows (FW/CADIS-Ω) to OpenMC's random-ray FW-CADIS pipeline for fusion fixed-source VR (deep penetration + streaming through ports/ducts/gaps), and do it as a staged, gated plan that measures benefit before building anything large. The conclusion of the design analysis below is: extend random ray — do not build a from-scratch Denovo-like Sₙ solver — and prove the benefit over a converged ℓ=0 FW-CADIS map before writing any angular-moment machinery.

Related: #106 (RR FW-CADIS distance auto-tune) and the convert_to_multigroup convergence auto-tune. Tier 0 of this plan is that work.

Background: the angular-importance gap

CADIS/FW-CADIS build weight windows (WW) from the scalar adjoint flux φ†. In streaming geometry the importance is strongly direction-dependent (a neutron headed into a duct mouth is worth orders of magnitude more than one headed at the wall beside it), and the scalar adjoint averages that away.

FW/CADIS-Ω (Munk & Slaybaugh, ANS M&C 2017, arXiv:1612.00793) replaces φ† with a forward-weighted adjoint scalar flux

φ†_Ω(r,E) = [ ∫_4π ψ(r,E,Ω) ψ†(r,E,Ω) dΩ ] / [ ∫_4π ψ(r,E,Ω) dΩ ]

and feeds it into the same FW-CADIS WW formula. The output is still ordinary space/energy WW (angle-informed, not angle-dependent), so the WW file and the CE transport step are unchanged. By spherical-harmonic orthogonality the numerator is a contraction of the angular flux moments of the forward and adjoint solves; the dominant correction is the ℓ=1 (current) term J·J†, which is large and positive where the forward and adjoint currents align along a streaming channel.

Key finding: today this is a no-op in OpenMC

OpenMC's random-ray solver is strictly P0 / isotropic-scatter and stores only the scalar (ℓ=0) flux:

• scatter source is a pure G×G sigma_s × scalar-flux contraction, no Legendre index — src/random_ray/flat_source_domain.cpp:121-139
• anisotropic MGXS is hard-rejected — src/random_ray/random_ray_simulation.cpp:79 (fatal_error("Anisotropic MGXS detected..."))
• the only stored "moments" are spatial linear-source gradients (MomentArray = Position), not angular — include/openmc/random_ray/moment_matrix.h:10-14

With only ℓ=0 available, φ†_Ω collapses exactly to φ† → CADIS-Ω is bit-identical to FW-CADIS. The missing input is angular flux moments (ℓ≥1) for both the forward and adjoint solves.

Decision: extend random ray, not build Sₙ

A 4-lens design review (performance / engineering-cost / physics-accuracy / OpenMC-strategy) was unanimous for extending RR. Rationale:

1. Ray effects cut for RR. Sₙ ray effects are worst in optically-thin ducts/voids/streaming with localized sources — the exact fusion regime — and the standard mitigations (first-collision source, quadrature rotation, finer order) are problem-specific knobs. RR/MOC samples direction continuously and is ray-effect-free by construction. (The CADIS-Ω paper itself flagged ray effects on its labyrinth case.)
2. Low-knob. Sₙ irreducibly adds quadrature order/type, spatial differencing, DSA/upscatter acceleration, negative-flux fixup, and a CAD→Cartesian voxel mesh. RR adds ~one internal knob (moment order, default P1) and reuses the existing CSG/DAGMC geometry and CE→MGXS pipeline.
3. Engineering / strategy. The FW-CADIS pipeline already exists and is tested on RR (forward+adjoint solves, set_fw_adjoint_sources at flat_source_domain.cpp:1234, WW export, 26 random_ray_* regression dirs incl. weightwindows_fw_cadis*); CADIS-Ω reuses it unchanged. A Denovo-like Sₙ is person-decades of zero-overlap code that doubles OpenMC's transport surface area and breaks the self-contained-MC ethos RR was built to honor.

Sₙ would only win if angular fidelity (high moments) dominated over ray-effect-freedom and the importance solve became the memory/wall-clock wall and the low-knob constraint were relaxed — a conjunction that does not hold here.

Important caveats / risks

• RR moment noise is not benign. φ†_Ω's numerator is a product of two noisy estimators, and RR has documented structured biases (large overestimates, negative fluxes) in exactly the low-flux void/duct regions of interest (docs/source/methods/random_ray.rst:533-561). RR's real edge over Sₙ is ray-effect-freedom + geometry/pipeline reuse — not that its error averages away. Needs denominator flooring, moment smoothing, and an auto-fallback to ℓ=0.
• Self-consistency. A physically consistent ℓ=1 current moment requires anisotropic scatter in the solve (lift the fatal_error, generate anisotropic MGXS — which grows the MGXS burden adding trigger to mgxs convert openmc-dev/openmc#3925/auto-tune is shrinking, and adds per-region angular source state). Prior art exists: anisotropic scattering (and linear sources) in the random ray method (Neame & Cosgrove, NSE 2024) — but validated on a thick reactor benchmark (C5G7), not sparse-hit fusion streaming.
• No MPI domain decomposition in RR (random_ray_simulation.cpp:269, rank-0 only). Adding forward+adjoint × P1/P3 moments multiplies the single-node footprint and could hit the memory ceiling on the first real fusion CAD model. This is the one factor that would flip the answer toward Sₙ at extreme scale — gate on a memory/runtime check.
• WW are applied angle-blind (space/energy bins only). CADIS-Ω is a scalar reweighting of the same map: in a straight duct it largely rescales importance ℓ=0 already places along the channel; the genuine redistribution is at bends/corners/junctions. So the benefit over a converged ℓ=0 map may be modest, and the published ~28× (Denovo vs coarse CADIS, one labyrinth) does not transfer.

Proposed staged plan

Tier 0 — ship auto-tuned ℓ=0 FW-CADIS (this largely exists).
Close #106 (runtime miss-rate-driven distance_inactive/distance_active) and fold in the convert_to_multigroup convergence auto-tune. This is today's genuine lowest-knob / highest-performance baseline — and possibly the whole answer.

GO / NO-GO gate (before any moment plumbing).
On a representative fusion single-bend duct problem, measure the FOM delta of an angle-informed importance map over the converged Tier-0 map. If the delta is small → stop; Tier 0 is the workflow.

Tier 1 — cheap post-hoc P1 current tally (only if Tier 0 underperforms).
Accumulate a per-region/per-group current moment J = ∫ ψ Ω dΩ on the existing P0 solve (accept P0-transport / P1-diagnostic inconsistency), for forward and adjoint, and form φ†_Ω from the ℓ=0 + ℓ=1 (J·J†) terms. Avoids lifting the fatal_error and avoids anisotropic MGXS. Reuses the existing directional-accumulation + locking + reduction pattern already proven for the linear-source spatial moments (src/random_ray/random_ray.cpp ~446-461, ~583-639). Substitute φ†_Ω at the single WW scoring point flat_source_domain.cpp:640-641 (scoped to the FW-CADIS flux path — that line also feeds total/fission scores). Add an auto-fallback to ℓ=0 when the J·J† moment is too noisy.

Tier 2 — full self-consistent CADIS-Ω (only if Tier 1 insufficient and measured FOM gap justifies).
Lift random_ray_simulation.cpp:79, add anisotropic-scatter transport + anisotropic MGXS (fold into the auto-tune, not a user knob), leverage the NSE-2024 anisotropic-RR machinery, extend to P3.

Parallel watch-item. Add MPI domain decomposition (or a hierarchical/coarsened importance mesh) to RR if the importance solve becomes the bottleneck.

The φ†_Ω contraction (P1 form to implement)

Expanding to first order, ψ(Ω) ≈ (1/4π)(φ + 3 J·Ω) with scalar flux φ = ∫ψ dΩ and current J = ∫ψΩ dΩ, and using ∫Ω dΩ = 0 and ∫(J·Ω)(J†·Ω) dΩ = (4π/3) J·J†, the CADIS-Ω quantity reduces (up to an overall constant absorbed by the WW 1/(2·max) normalization) to:

φ†_Ω(r,E)  ∝  φ†_adj(r,E)  +  3 · [ J_fwd(r,E) · J†_adj(r,E) ] / φ_fwd(r,E)

• First term = the ordinary scalar adjoint flux → recovers plain FW-CADIS exactly (the N=0 baseline / regression check).
• Second term = the angular correction: forward current · adjoint current, normalized by forward scalar flux; large and positive where both currents align along a streaming channel.
• Inputs per (region, group): φ_fwd, φ†_adj (both already computed) plus the NEW current vectors J_fwd, J†_adj.
• Watch: φ_fwd → 0 in deep voids (needs a floor); J·J† can be negative when currents anti-align (clamp φ†_Ω ≥ the ℓ=0 term, or floor at 0).
• P3 would add Σ_{ℓ=2,3} φ_ℓm φ†_ℓm; defer unless P1 proves insufficient.

Implementer's map (code pointers, branch making-wheel-3)

Task	Location	Note
Accumulate new current moment J (fwd + adj)	src/random_ray/random_ray.cpp:440-462 (flat-source scalar accum) + void/linear variants	Mirror the linear-source spatial-moment accumulation at random_ray.cpp:600-646 (MomentArray, locking, MPI reduce). Ray direction u() and angular_flux_[g] are already in scope.
New per-region/group storage for J	include/openmc/random_ray/source_region.h:367-385 (beside scalar_flux_*)	A 3-vector per (region, group), × {fwd, adj}.
Snapshot forward J before the adjoint overwrites it	flux_swap (src/random_ray/source_region.cpp:175-181)	Mirror how set_fw_adjoint_sources (flat_source_domain.cpp:1234-1273) snapshots the forward scalar flux.
Adjoint construction (exists, no change)	transpose_scattering_matrix (flat_source_domain.cpp:1395-1414)
Form φ†_Ω & substitute	flat_source_domain.cpp:640-649 (random_ray_tally, SCORE_FLUX)	Scope to the FW-CADIS flux path — that line also feeds total/fission scores.
CADIS inversion / WW formula (unchanged)	src/weight_windows.cpp:624-658	WW = 1/φ†, normalized by 1/(2·max).
New method string fw_cadis_omega	Python openmc/weight_windows.py:614; C++ src/weight_windows.cpp:790-814 (branch after :798)	Or a moment-order option defaulting to 0 (= FW-CADIS).
Tier 2 only: allow anisotropic scatter	lift fatal_error at random_ray_simulation.cpp:79	Needs anisotropic MGXS + per-region angular source state.
Scaling watch-item	no MPI domain decomposition: random_ray_simulation.cpp:268-271	fwd + adj × moments live on rank 0 only.
WW export/import (unchanged)	to_hdf5 weight_windows.cpp:739-747; import simulation.cpp:162

GO/NO-GO benchmark & metric

• No duct/streaming geometry exists in-tree. The RR FW-CADIS regression tests tests/regression_tests/weightwindows_fw_cadis{,_local,_mesh}/ all use a simple three-region void-layered cube (random_ray_three_region_cube() in openmc/examples.py ~:875), not a true duct. The gate needs a new single-bend duct benchmark — extend the random_ray_* model factories in openmc/examples.py:875+, modeled on the labyrinth / single-bend duct from the references.
• The metric is built in: tally.figure_of_merit = 1/(rel_err² · simulation_time) (openmc/tallies.py:955-962), with tally.std_dev (:589-609) and tally.mean. Gate = FOM(angle-informed WW) vs FOM(converged ℓ=0 FW-CADIS WW) at the deepest tally, both consuming their WW in the same CE transport.

Acceptance criteria

• Tier-0 auto-tuned ℓ=0 FW-CADIS baseline runs knob-free on a fusion duct problem (ties to Auto-tune random ray for FW-CADIS: physics-aware distance_inactive/distance_active (toward zero-knob weight-window generation) #106)
• Documented FOM comparison: angle-informed vs converged ℓ=0 on that problem (GO/NO-GO)
• If GO: Tier-1 P1 current tally with N=0 reproducing FW-CADIS bit-for-bit as a regression baseline, plus auto-fallback to ℓ=0
• Noise/robustness check at the deepest tally site / duct exit

Open questions / decisions to make

• Moment order: start at P1 (current only). Is P3 ever worth it, given RR is ray-effect-free?
• Consistency (the Tier-1→Tier-2 gate): does the cheap post-hoc P1 tally on a P0 solve yield a usable φ†_Ω, or is self-consistent anisotropic transport required? (Empirical.)
• Denominator flooring: what floor for φ_fwd → 0 in deep voids? (reuse the ZERO_FLUX_CUTOFF pattern, or a relative floor.)
• Negative correction: clamp φ†_Ω ≥ the ℓ=0 term, or floor at 0, when J·J† < 0?
• Moment denoising: is coarse-mesh averaging of the currents needed, and at what length scale?
• Auto-fallback trigger: what rel-err threshold on the J·J† moment falls back to ℓ=0 FW-CADIS?
• Memory: does fwd+adj × (scalar + 3-vector current) per region·group fit a single node (no MPI DD) on a real fusion CAD model?

Side discussions & dead ends (for the record)

Approaches considered and rejected:

• Build an in-tree Denovo-like Sₙ solver — rejected (see Decision): ray effects are worst in the fusion ducts/voids we target, it explodes user-facing knobs, it's person-decades of zero-overlap code, and it breaks the self-contained-MC ethos RR exists to honor.
• Couple to an external Sₙ code (Denovo/ADVANTG) — rejected: reintroduces exactly the external deterministic dependency RR was built to eliminate; Denovo is SCALE-gated and drives MCNP, not OpenMC.
• "Just swap φ†_Ω in" as a quick patch — dead end: there are no ℓ>0 angular moments in RR to swap, and the single WW scoring line (flat_source_domain.cpp:640-641) also feeds total/fission scores, so a blind variable replacement would corrupt other RR tallies. The substitution must be scoped to the FW-CADIS flux path.

Useful side insights that shaped the plan:

• Backend-invariance collapses the RR-vs-Sₙ question. CADIS-Ω emits ordinary space/energy WW, so the FOM-dominating CE transport step is identical regardless of backend. The decision reduces to "which backend gives the best importance map per unit of one-time preprocessing effort" — which is why solver raw speed matters far less than map quality in the ducts.
• The angular correction is dominated by ℓ=1. By spherical-harmonic orthogonality ∫ ψ ψ† dΩ = Σ_ℓm φ_ℓm φ†_ℓm; in a streaming duct the dominant correction is the current term J·J† (both currents align with the duct axis). So low order (P1) likely captures most of the benefit — this is what makes the cheap Tier-1 path viable and means high-order Sₙ fidelity is not where the win is.
• "Add P1 to RR" has two distinct meanings (often conflated): a cheap post-hoc current tally on the existing P0 solve (Tier 1, inconsistent but possibly sufficient) vs self-consistent anisotropic P1 transport (Tier 2, needs anisotropic MGXS + lifting the fatal_error). Keeping these separate is what keeps Tier 1 cheap.
• Course-correction during analysis: RR moment noise was initially framed as "benign zero-mean that averages out" — that's wrong. RR has structured biases (huge overestimates, negative fluxes) in low-flux void/duct regions, so the correct justification for RR is ray-effect-freedom + pipeline reuse, and the plan must include moment smoothing + ℓ=0 auto-fallback.
• WW are applied angle-blind, so CADIS-Ω only reshapes the scalar map — it rescales importance ℓ=0 already places in a straight duct and redistributes meaningfully only at bends/corners. Hence the benefit-vs-converged-ℓ=0 question, and the GO/NO-GO gate, rather than assuming the paper's ~28× (which was Denovo vs coarse CADIS on one labyrinth and does not transfer).

Related in-flight work

• Auto-tune random ray for FW-CADIS: physics-aware distance_inactive/distance_active (toward zero-knob weight-window generation) #106 — RR FW-CADIS distance auto-tune; this is Tier 0.
• MGXS auto-tune — branch trigger-mgxs-generation (commit 9f1dea8f0): convert_to_multigroup(nparticles=None) converges via a rel-err trigger (knob-free redo of the abandoned PR adding trigger to mgxs convert openmc-dev/openmc#3925). Tier 2's anisotropic-MGXS requirement must fold into this, not become a new knob.

References

• Munk & Slaybaugh, "FW/CADIS-Ω: An Angle-Informed Hybrid Method for Deep-Penetration Radiation Transport," ANS M&C 2017 — https://arxiv.org/pdf/1612.00793
• Munk & Slaybaugh, "Review of Hybrid Methods for Deep-Penetration Neutron Transport," Nucl. Sci. Eng. 193(10), 2019 — https://doi.org/10.1080/00295639.2019.1586273
• Neame & Cosgrove, "Linear Sources and Anisotropic Scattering in the Random Ray Method," Nucl. Sci. Eng. 2024 — https://doi.org/10.1080/00295639.2024.2394729
• Auto-tune random ray for FW-CADIS: physics-aware distance_inactive/distance_active (toward zero-knob weight-window generation) #106 — RR FW-CADIS distance auto-tune (Tier 0)

[jon-proximafusion]

Tier 1 design verified against develop @ f018524 (post openmc-dev#3962)

A full code sweep confirms the Tier 1 slot-in is feasible with no changes to the WW file format, the CE transport step, or the statepoint format. Findings that update or sharpen the implementer's map above, verified June 2026:

openmc-dev#3962 landed and helps. There is now a three-state solve enum RandomRaySolve {FORWARD, FORWARD_FOR_ADJOINT, ADJOINT} (constants.h:371). During an FW-CADIS run the forward solve runs as FORWARD_FOR_ADJOINT, its statepoint/tallies/plots are tagged forward.*, and the weight-window export is skipped until the adjoint solve completes (simulation.cpp:203-210). One RandomRaySimulation/domain object persists across both solves, which is exactly the lifetime needed for a forward-moment snapshot.

Where the forward flux dies today. set_fw_adjoint_sources() (flat_source_domain.cpp:1242) builds the adjoint source as 1/phi_fwd from scalar_flux_final_ and then zeroes it in place. adjoint_reset() (source_region.cpp:231) zeroes everything else but deliberately leaves scalar_flux_final_ alone. So the snapshot point for phi_fwd and J_fwd is inside set_fw_adjoint_sources, before the zeroing.

Current accumulation is a four-kernel change, not one. random_ray.cpp has four leaf kernels (flat, flat void, linear, linear void). All accumulate a per-group quantity into scalar_flux_new under a per-region lock with the ray direction u() in scope, so the J estimator is current_new(g) += u() * delta_psi[g] (or u() * psi * L in void kernels) added inside the existing lock. Normalization differs by source shape: flat divides by sigma_t*V, linear divides by V only, both already routed through the virtual set_flux_to_flux_plus_source, so the J normalization hooks in the same two overrides. The linear-source spatial-moment machinery (delta_moments_, MomentArray = Position) is a line-for-line template.

Storage is SoA plus dynamic discovery. Per-(region,group) arrays live in SourceRegionContainer and must also be threaded through SourceRegion (AoS used for mesh-subdivided regions discovered mid-ray), SourceRegionHandle, push_back, assign, get_source_region_handle, and adjoint_reset (source_region.cpp). Missing any of these compiles fine and silently drops J for subdivided regions. New arrays: current_new_, current_t_ (working), current_fwd_, scalar_flux_fwd_ (forward snapshots), gated behind a static flag so the default path allocates nothing.

The substitution point moved (and must be scoped by tally index, not mesh). There is no dedicated FW-CADIS scoring path. The WWG creates a plain flux tally (weight_windows.cpp:852) scored generically by random_ray_tally (flat_source_domain.cpp:640-651). Two traps: (1) create_tally reuses an existing user MeshFilter when one matches, so scoping the omega substitution by mesh would corrupt ordinary user tallies on the same mesh; scope by task.tally_idx matched against omega-enabled generators. (2) The tally scores the per-iteration scalar_flux_new each active batch and the WW update (last batch only in RR mode) uses the realization mean, so omega must be applied per batch as a multiplicative factor built from running-accumulated moments (accumulate_iteration_flux runs before random_ray_tally, so totals include the current batch).

Adjoint sign convention. The adjoint solve reuses the forward kernel with a transposed scattering matrix, so it solves for chi(omega) = psi_adj(-omega). The accumulated adjoint current is therefore the reversed physical adjoint current, J_adj_phys = -J_chi. The P1 contraction to implement is

omega(r,g) = 1 - 3 * (J_fwd . J_chi) / (phi_fwd * phi_adj),   phi_omega = phi_adj * omega

clamped to a positive range (about [0.05, 4]; the P1 theoretical range of the bracket is [-2, 4]) with omega = 1 fallback wherever phi_fwd is below the existing ZERO_FLUX_CUTOFF-scaled floor, phi_adj is non-positive, or the region is flagged is_small. Sanity check for deep penetration: both physical currents point outward (toward low-phi regions where the 1/phi adjoint source is largest), so J_fwd . J_adj_phys > 0 and omega > 1, growing with depth as the angular flux becomes more forward-peaked.

Python/XML surface. The method whitelist is ('magic', 'fw_cadis') at openmc/weight_windows.py:614, and the targets gate at :733-737 is an exact string compare that would reject a new method string until widened. The method string never crosses the C API (openmc.lib only exposes update_magic), so no ctypes changes. update_weights treats every non-MAGIC method as CADIS-family, so a new FW_CADIS_OMEGA enum value reuses the existing invert-and-normalize update unchanged.

No serialization changes. No SourceRegion data is written to statepoint or HDF5 anywhere, so J arrays need zero I/O changes. The WW file format is untouched.

Next step (in progress)

Implementing the above on branch fw-cadis-omega as a new WeightWindowGenerator(method='fw_cadis_omega'), with the default fw_cadis path bit-for-bit unchanged. First benchmark will be a deliberately adversarial control case for the GO/NO-GO gate: a 2 m thick ordinary-concrete box bioshield with a 14.1 MeV ring source in the central cavity and no penetrations, comparing analog vs fw_cadis vs fw_cadis_omega figure of merit on an outer-surface tally. Expectation from the caveats above: in a homogeneous shield with no streaming paths, both currents align with the shield normal everywhere, so omega is a smooth monotone steepening of the same importance map and the FOM delta should be modest. The interesting result is measuring how modest; the streaming/duct case where the correction genuinely redistributes importance comes after.

[jon-proximafusion]

Tier 1 implemented and benchmarked: 2 m bioshield control case

The Tier 1 plan above is now implemented on branch fw-cadis-omega (commits 4b381c13d + 3273fe8be, based on develop @ f018524) as WeightWindowGenerator(method='fw_cadis_omega'):

• P1 net current J accumulated in all four ray kernels (flat/linear, void/non-void) under the existing per-region lock, in both the forward and adjoint solves. Storage gated behind a static flag so the default path allocates nothing.
• Forward flux and current snapshotted inside set_fw_adjoint_sources() before the working arrays are consumed.
• The WW tally score is multiplied by omega = clamp(1 - 3 (J_fwd . J_chi) / (phi_fwd phi_adj), 0.05, 4) (the minus sign because the kernel-accumulated adjoint current is the reversed physical adjoint current), scoped by tally index to omega generator tallies only, adjoint solve only, with fallback to 1 for small regions, sub-cutoff forward flux, or non-positive adjoint flux.
• All 5 existing FW-CADIS regression tests pass bit-for-bit on the modified build; a new weightwindows_fw_cadis_omega regression test and a unit test for the method string are included. RR solve cost: about +11 percent; the CE transport step is untouched.

Sign convention verified empirically on a 1-group 1D absorber slab: the omega/fw_cadis WW shape ratio is 1.0 at the source, drops to ~0.42 through the bulk (constant amplification ~2.4, as expected for the asymptotic angular shape of a homogeneous 1-group slab), and relaxes to ~0.94 at the deep vacuum boundary where the adjoint current reverses at the adjoint source peak.

Benchmark: 2 m ordinary-concrete box bioshield, 14.06 MeV ring source, no penetrations

Setup: cubic void cavity (half-width 2 m), 2 m NIST ordinary concrete shield (outer half-width 4 m, vacuum), DT ring source (r = 1.5 m, isotropic, 14.06 MeV) in a thin annular void cell. Pipeline: convert_to_multigroup(material_wise, 10 groups) then RR forward+adjoint on a 20 cm source-region mesh (flat source, naive volume estimator, 30k rays x 150 batches) then CE transport with the generated WW. Deep tally: flux in the outermost 10 cm concrete shell. Global tally: 40^3 mesh flux.

The omega map itself is physically textbook. The correction concentrates exactly where angular information lives: WW lower bounds in the 6 to 15 MeV groups (which carry the deep penetration) are lowered by up to ~30 percent through the shield bulk relative to plain FW-CADIS, while thermal groups are essentially untouched (ratio ~0.99, the flux there is isotropic so J.J_adj ~ 0). Sign correct everywhere.

But the FOM verdict on this geometry is a wash. Four seeds per method, 40M CE particles each (seed pairs run under identical conditions):

case	det tally rel err (seeds 2-4 avg)	det FOM geo-mean (paired seeds 2-4)	mesh bins < 10 percent err
analog	no scores at all (12M particles, zero hits in outer shell)	0	26 percent (cavity region only)
fw_cadis	4.2 percent	0.878	34.8 percent
fw_cadis_omega	4.1 percent	0.874	35.1 percent

Paired per-seed FOM ratios (omega/fw_cadis): 0.84, 0.87, 1.36; geometric mean 1.00. Depth-resolved mesh rel err profiles are identical (ratio 1.00 to 1.01 in every 40 cm shell out to the surface). Detector means agree across all runs (1.09e-6 vs 1.11e-6), confirming both estimators are unbiased.

A methodological caution worth recording: the first single-run comparison showed an apparent 3x FOM gain for omega. That was entirely a heavy-tail fluctuation in the fw_cadis variance estimate (its seed-1 rel err was 9.5 percent where the other three seeds gave ~4.2 percent at identical particle counts). Single-run FOM comparisons in tail-dominated deep-penetration problems can easily produce spurious 2 to 3x swings in either direction; seed replication is mandatory for the GO/NO-GO gate.

Interpretation and answer to the gate question for this geometry

For a thick homogeneous shield with no penetrations, angle-informed weight windows do not make the simulation converge measurably faster than converged scalar FW-CADIS. This is the outcome the caveats above predicted: with no streaming paths, the forward and adjoint currents align with the shield normal everywhere, so the P1 correction is a smooth monotone reshaping (bounded at ~30 percent on the WW bounds, largest in the hardest groups) that the CADIS normalization mostly absorbs. The scalar adjoint already encodes essentially all of the importance structure this problem has.

Equally important: plain FW-CADIS is doing the heavy lifting. Analog CE cannot score the outer shell at all, while either WW set reaches ~4 percent on the deep tally in ~10 minutes.

Status of the staged plan: Tier 1 machinery is built, verified, and cheap, but the unpenetrated bioshield is the control case where the null result is expected. The GO/NO-GO measurement that can actually distinguish the methods is the streaming case: the same bioshield with a penetration (straight duct first, then a single-bend duct), where the correction genuinely redistributes importance at the duct mouth and bend rather than rescaling it. That benchmark is the natural next step on this branch.

[jon-proximafusion]

GO/NO-GO gate result: Protvino-style labyrinth (streaming with two bends)

Second benchmark on branch fw-cadis-omega, modeled on the ALARM-CF-AIR-LAB-001 Protvino labyrinth (exact drawings are NEA-DVD-only, so this is faithful in spirit): a U-shaped 3-section corridor (2 m wide x 2.6 m high, 6 m legs, two 90 degree bends) through a concrete monolith with >= 1 m walls, Cf-252 Watt point source at the entrance of leg 1, five 10 cm volumetric detectors down the legs (D1 end of leg 1 through D5 at the leg 3 dead end, ~6 decades of importance along the path). Same pipeline as before: 10-group material-wise MGXS, RR forward+adjoint on a 20 cm source-region mesh, CE transport, 3 paired seeds per method plus analog.

The omega map now shows genuine redistribution, exactly as the theory promises. Along the corridor centerline the omega/fw_cadis WW lower-bound ratio is non-monotone: 0.56 to 0.94 in the straight streaming runs (more splitting where forward and adjoint currents align down the leg) and 2 to 4 in the bend pockets and at the dead end (less splitting where the forward current points into a wall). The structure grows with depth and energy, spanning ratio 0.59 to 4.2 in leg 3.

But the redistribution is counterproductive. FOM results (geo-means over seeds, paired ratios):

detector	analog	fw_cadis	fw_cadis_omega	paired ratio omega/fw_cadis
D1 (end leg 1)	1893	7.5	5.8	0.77
D2 (mid leg 2)	84	4.2	2.8	0.66
D3 (start leg 3)	10.2	6.5	2.1	0.32
D4 (mid leg 3)	5.2	5.6	1.9	0.35
D5 (dead end)	3.0	5.1	1.8	0.34

All nine paired seed comparisons beyond the second bend favor plain fw_cadis (ratios 0.11 to 0.79). Detector means agree across all three methods to within 2 percent everywhere, so both estimators are unbiased; this is purely an efficiency loss. Note also that global FW-CADIS itself only beats analog by ~1.7x at the deepest detector and loses to analog at D1-D3: at ~3 decades of streaming attenuation this maze is still largely analog-friendly.

Interpretation: the P1 contraction misranks importance exactly at the feature it targets. Transport around a bend proceeds through wall scattering in the bend pocket: neutrons reach leg 3 by hitting the back wall of the junction and re-emerging. There the angular flux is bimodal (incident stream into the wall plus albedo around the corner), the net currents anti-align, and the P1 correction J.J_adj < 0 suppresses importance in precisely the cells whose wall-albedo events feed the next leg. The particles that matter get rouletted or under-split at the bends, and the variance moves downstream. This is the failure mode flagged in the plan's caveats: the published ~28x CADIS-Omega gain came from a deterministic S16 Denovo adjoint with full angular fidelity versus a coarse CADIS baseline; a P1 net-current moment from random ray does not carry enough angular information at bimodal-flux locations, and at bends it gets the sign of the redistribution wrong.

Combined verdict for the staged plan

• Unpenetrated 2 m bioshield: omega = wash (paired FOM ratio 1.00).
• 3-section labyrinth: omega = consistent ~3x FOM loss beyond the second bend.

Tier 1 (post-hoc P1 current on the P0 solve) is a NO-GO as an FOM improver in both regimes tested. The machinery is correct, verified, and cheap, and the maps behave exactly as P1 theory says they should; the problem is that P1 is the wrong truncation for the streaming features that motivate angle-informed VR. Options this leaves, in rough order of my read on value:

1. Accept Tier 0 as the answer (converged scalar FW-CADIS plus auto-tuning, Auto-tune random ray for FW-CADIS: physics-aware distance_inactive/distance_active (toward zero-knob weight-window generation) #106): it is unbeaten in both benchmarks here.
2. Before abandoning omega entirely, two cheap experiments on this branch could isolate the failure: (a) clamp the correction to amplification-only (omega >= 1, i.e. never suppress below the scalar map) which keeps the streaming-run boost but removes the bend-pocket penalty; (b) tighten the clamp range (e.g. [0.5, 2]). Both are one-line constant changes.
3. Tier 2 (self-consistent higher-order moments, P3+, anisotropic scatter in the solve) is what the theory actually requires at bends, with the cost and knob growth the plan already documents. The Tier 1 result strengthens the case that nothing short of it will beat the scalar map, and weakens the case that it is worth building.

Benchmark model and scripts: omega_bench/labyrinth.py, gen_ww_lab.py, run_ce_lab.py, analyze_lab.py (untracked, repo root).

[jon-proximafusion]

Tiers 2 and 3 implemented and benchmarked: higher-order moments and true angle-dependent windows

Following the Tier 1 NO-GO above, two escalations were implemented on branch fw-cadis-omega (now 4 commits, pushed; defaults remain bit-identical to develop) and benchmarked on the same labyrinth with the same paired seeds.

Tier 2: P2/P3 contraction plus configurable clamp (commit c42cd33be)

The random ray solver now optionally accumulates the l = 2 and l = 3 orthonormal real spherical harmonic angular flux moments (12 components per source element, forward and adjoint), and the omega factor generalizes to the full contraction 1 + 4pi sum_(l>=1,m) (-1)^l psi_lm chi_lm / (phi_fwd phi_adj), with the parity factor from the adjoint solve computing psi_adj(-omega). New update_parameters on the generator: order (1 to 3), clamp_min, clamp_max. Setting clamp_min: 1.0 gives an amplification-only variant that can never suppress importance below the scalar map.

Result: P3 made it worse, and the clamp experiments prove the mechanism. Paired FOM ratios vs scalar fw_cadis (geo-mean over 3 matched seeds):

detector	omega P1	omega P3	omega P3 amp-only	omega P1 amp-only
D1 (end leg 1)	0.77	0.39	0.39	0.54
D2 (mid leg 2)	0.66	0.28	0.63	0.66
D3 (after bend 2)	0.32	0.18	0.24	0.29
D4 (mid leg 3)	0.35	0.22	0.27	0.35
D5 (dead end)	0.34	0.28	0.37	0.48

The WW map probe explains it: where P1 suppressed the bend pockets ~2.6x relative to the legs, P3 suppresses them 10 to 16x. Higher moments represent the anti-alignment of the forward and adjoint angular fluxes at the pockets more faithfully, and that is precisely the wrong thing for an angle-blind window to act on: the particle flying into the bend wall is about to become important via albedo, but the scalar window can only kill it. Forbidding suppression (amp-only) recovers most of the P3 penalty, confirming the diagnosis, but never beats scalar.

Tier 3: true angle-dependent (octant) weight windows (commit f0fb3fb07)

Since the ceiling was the angle-blind window application itself, the stronger thing was built: WeightWindows gains a num_angle dimension (8 direction octants), the transport lookup selects the octant matching the particle's direction of flight at every WW checkpoint, and the octant-resolved importance is reconstructed from the adjoint scalar flux plus the P1 (and optionally P3) moments evaluated at reversed octant centers. num_angle round-trips through weight_windows.h5, XML, and the Python API. Enabled with update_parameters={'angular_bins': 8}. This is the AVATAR idea (Van Riper et al., LANL 1997) driven by random ray moments instead of a deterministic Sn adjoint.

The generated maps finally express the right physics. At the leg 2 / leg 3 junction the most-split octants are +x,-y, i.e. around the bend, with 57x anisotropy; mid-leg-2 favors +x; leg 3 favors -y. This is exactly the information every scalar variant could not act on.

And the transport result still loses. Seed 3, same conditions as all other variants:

method	D3 FOM	D5 FOM	active time
analog	10.2	3.0	700 s
fw_cadis (scalar)	6.5	5.1	~1500 s
angular octant P1	1.1	1.3	3462 s
angular octant P3	2.6	1.4	1829 s

The angular P1 run took 2.3x longer at the same source particle count. The cause is split/roulette churn: octant windows place 10x to 57x weight-bound discontinuities between directions in the same cell, so every scattering event that carries a particle across an octant boundary triggers a large split or a roulette. The churn costs both time and variance; even the timing-free variance-per-particle comparison loses to scalar. This is the classic reason angle-dependent weight windows never became a production default. The standard mitigation, damping the octant anisotropy (capping the inter-octant ratio at ~2), is a one-line change in the new angular update branch and remains the one untested idea.

Full comparison

All detector means agree across all eight methods to within ~3 percent, so every estimator is unbiased; the differences are purely efficiency. Benchmark scripts and all figures live on branch omega-bench-results.

On analog winning at D1 to D3

Worth recording so the ranking is not misread: this maze only attenuates ~2.6 decades end to end (flux 2.76 to 0.0062 per source particle), the corridors are open air, and the detectors are 10 cm slabs spanning the full corridor cross section. Analog reaches D5 by ordinary random walk with cheap histories, while global FW-CADIS pays WW lookups, splitting trees, and effort converging the concrete island that no detector cares about. Analog falls ~2.5 orders of magnitude over the path while fw_cadis stays flat at FOM ~5; the curves cross between D3 and D5, and by the dead end scalar FW-CADIS leads ~1.7x. The published hard cases (5A with polyethylene plates in the second section, 6B with covered source and dead-end) add the decades that make analog collapse; reproducing 5A is a ~20 line change to labyrinth.py if a deep-regime separation is ever needed.

Closing verdict for this issue

The hypothesis chain is now fully tested, in code, on a streaming benchmark:

1. Angle-informed P1 (the paper's FW/CADIS-Omega): loses ~3x at bends. Suppresses wall-bound particles that carry the bend transport.
2. Angle-informed P3: loses ~4 to 5x. Better angular fidelity in the contraction sharpens a correction the angle-blind window cannot use.
3. Amplification-only clamps: recover to ~half of scalar. Confirms the suppression mechanism; still never a win.
4. True angle-dependent octant windows: express the correct directional physics and lose ~3 to 4x anyway, to inter-octant split/roulette churn.

Converged scalar FW-CADIS (Tier 0, with the #106 auto-tuning) survived every angular scheme implemented here and remains the recommendation for the fusion VR stack. The branch is kept as a complete, tested reference implementation (moment accumulation, omega contraction to P3, angular window storage/lookup/IO) should the damped-anisotropy variant or a different application of the moments ever be worth a look.

[jon-proximafusion]

Status update and the benchmark still needed

Where things stand (2026-07-03)

Everything proposed in this issue up to and beyond the original Tier 2 has now been implemented and measured. Inventory:

Code (branch fw-cadis-omega, 4 commits off develop f01852411, defaults bit-identical to develop, all existing FW-CADIS regression tests passing plus new omega regression and unit tests):

commit	capability
4b381c13d	Tier 1: P1 current accumulation in all 4 RR kernels (fwd+adj), forward snapshot, omega factor on the WW tally, method string fw_cadis_omega
3273fe8be	Review fixes, tests, docs
c42cd33be	Tier 2: l = 2, 3 spherical harmonic moments, full P3 contraction, order / clamp_min / clamp_max update parameters
f0fb3fb07	Tier 3: angle-dependent octant weight windows (angular_bins: 8), direction-of-flight lookup in transport, octant importance from adjoint moments, num_angle through HDF5/XML/Python

Benchmarks run (scripts and figures on omega-bench-results):

1. 2 m concrete box bioshield, 14.06 MeV ring source, no penetrations (~7 decades, pure deep penetration). Analog scores nothing; scalar FW-CADIS and omega P1 statistically tied (paired FOM ratio 1.00, 3 seeds).
2. Protvino-style 3-section labyrinth, Cf-252, two 90 degree bends (~2.6 decades, pure streaming). Results across all eight methods, paired seeds, all estimators unbiased (means agree to ~3 percent):

method	verdict at D3 to D5 (vs scalar fw_cadis)
analog	wins at D1 to D3, loses ~1.7x at D5 (problem too shallow for VR)
fw_cadis (scalar, Tier 0)	reference, FOM ~5 flat along the whole path
omega P1 (Tier 1)	0.32 to 0.35
omega P3 (Tier 2)	0.18 to 0.28 (worse: sharpens bend-pocket suppression 10 to 16x)
omega P1/P3 amp-only (clamp_min 1)	0.24 to 0.48 (recovers half; confirms the suppression mechanism)
angular octant P1/P3 (Tier 3)	0.2 to 0.5, with 1.2 to 2.3x longer runtimes from inter-octant split/roulette churn, despite directionally correct maps (57x anisotropy pointing around the bend)

Conclusion so far: converged scalar FW-CADIS is unbeaten by every angular scheme implemented. But there is a gap in the evidence, and it motivates one more benchmark before this issue can close cleanly.

The benchmark we still need: deep streaming (Protvino Case 5A analog)

Every angular-vs-scalar comparison so far sits in a regime where analog transport still works (the labyrinth is only ~2.6 decades, open air, large detectors). The regime this issue was opened for is the one where VR is mandatory: fusion bioshields with penetrations, 6+ decades, where the question "does angular information help" has never actually been measured by us. The churn cost and the directional gain of octant windows scale differently when the analog contribution is nil, so the labyrinth ranking does not automatically transfer.

The original Protvino benchmark anticipates exactly this: its hard configurations block the corridor. Proposed spec, a ~20 line change to omega_bench/labyrinth.py:

Geometry. Same 3-section labyrinth (2 x 2.6 m corridor, 6 m legs, 1 m concrete walls) plus two polyethylene plates (0.92 g/cm3 CH2, 10 cm thick) spanning the full leg 2 cross section at x = 400 cm and x = 600 cm (Case 5A style). Optionally a third configuration with a 10 cm poly cap over the source (Case 6B style). Fast neutrons must now penetrate hydrogenous slabs twice on top of the two bends, pushing D5 to a target of >= 6 decades below D1.

Detectors. The existing D1 to D5 plus one detector immediately downstream of each plate (D2a, D2b), so the attenuation budget is resolved plate by plate.

Protocol (encoding the lessons from the runs above):

1. Verify the attenuation with the cheap RR forward solve before any CE compute; if D5 is under 6 decades, thicken the plates to 15 cm.
2. Analog control first: expected to produce zero scores at D4/D5 at feasible particle counts. That is the entry condition for the deep regime; if analog still scores, the benchmark is not hard enough.
3. Method matrix: analog, scalar fw_cadis (Tier 0), omega P1, angular octant P1, angular octant P3, and the one variant not yet tested anywhere: damped-anisotropy octant windows (cap the inter-octant bound ratio at ~2, a one line change in the angular branch of update_weights), which targets the churn failure directly.
4. Three paired seeds per method, runs executed sequentially at a fixed thread count (the 3-way concurrent runs above inflated wall-clock unevenly), FOM on active time, plus the timing-free variance-at-equal-particles comparison as a cross check, plus the means-agreement check (<= 3 percent).

GO/NO-GO. Angular (any variant) is declared useful if its paired FOM ratio vs converged scalar FW-CADIS is >= 1.5 at D4/D5, consistent in sign across all three seeds. Anything less closes this issue with a final verdict: scalar FW-CADIS plus the #106 auto-tuning is the production answer, and the branch remains a reference implementation.

Estimated cost: one RR verification solve (~5 min), then roughly 6 methods x 3 seeds of CE at the labyrinth run size, executable overnight on the current machine.