Naive Flux Linear Source Consistency Fix

[jtramm]

Overview

This PR fixes a statistical inconsistency in the linear source flux update. Basically, the naive estimator currently uses per-batch "naive" volumes, but simulation-averaged centroids.

The naive volume treatment updates the flux from a single batch's rays, but the linear source is anchored to the simulation-averaged centroid, so the two halves of the update describe different sets of tracks. The difference (the source gradient dotted with the batch's centroid fluctuation) was silently absorbed by the flux estimate as noise, which breaks the exactness the naive treatment promises and can ignite negative fluxes in optically thin scatter-fed regions. About 15 lines of logic that affect linear source mode only.

Independent of the adaptive volume estimator PR (#4110), and reviewable in either order. That said, as I developed PR #4110 first and it also improves performance, I measure results against a "baseline" value generated by the PR #4110 branch, which is where the auto results below come from.

The fix

For source regions using the naive flux estimator, the flux update now adds back the source as the batch's tracks actually integrated it, evaluated at the batch centroid rather than the accumulated one. The term is applied only when a region updates with its own batch (naive) volume, which is exactly where it makes the update an exact per-batch track identity (the flux inherits the sign of the angular fluxes it averages). Regions updating with the simulation-averaged volume keep the original treatment, which already correctly matches the simulation-averaged volume with a simulation-averaged centroid.

Results

Stability on the shielding test problem (the "irradiation vault hall" problem from #4110, naive estimator, linear sources, 2,000 rays, 100 inactive and 100 active batches), negative slow-group bins of 85,750 and worst amplitude:

	negative bins	worst amplitude
before	5,724	−3.8e-1
after	72	−2.7e-3

Accuracy on the C5G7 overlay-mesh eigenvalue case (linear_xy, 200 rays, 102x102 mesh), where the naive estimator's linear-mode noise was previously severe:

	k dev (pcm)	AAPE
naive, before	−5,514	4.09%
naive, after	−3,524	1.58%
auto (default), before	+45	0.65%
auto (default), after	+45	0.65%

The high error on C5G7 is expected as the naive estimator is being used. For accurate results, the simulation averaged, hybrid, or adaptive estimators are far better. However, these results do demonstrate that the naive estimator becomes much more accurate, which is the goal, and does allow for hybrid/adaptive estimators to maintain accuracy in cases where more cells are demoted under low ray density conditions etc.

The auto rows were measured with this fix combined with the adaptive estimator PR: the default path is untouched at measurement precision, since under the adaptive family the term only fires in the rare batch where a transiently hit-starved region still carries its fitted gradients. Pinned naive and hybrid linear runs see the full correction.

Compatibility and testing

• Flat source mode is bit-for-bit unaffected (the term is identically zero), and the simulation-averaged linear reference results are untouched byte for byte, confirming the original arithmetic survives where the term does not apply.
• Reference results are regenerated for the linear-mode tests whose estimators receive the term (naive, and hybrid through its external-source and hit-starved regions). A new random_ray_linear_source_stability regression test pins the fix on a starved naive linear configuration.
• The methods documentation gains one paragraph deriving the batch-sampled source and stating the pairing rule.

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[jtramm]

Merged develop, which now includes the adaptive volume estimators (#4110). One interaction with #4110 needed a fix. The strict adaptive estimator rescues a negative flux iterate by rebuilding the batch-volume update from an additive term. With this PR, that term depends on which volume the update used. The additive term now takes an explicit batch-volume flag, so the flux update and the rescue always agree.

I also made two other small changes: 1) the regression test now uses 20 rays and no overlay mesh, since the old mesh matched the lattice and no region was ever missed, and 2) the methods note now says the term omitted on the simulation-averaged path is zero-mean and shrinks with the batch count.

[paulromano]

All looks good; thanks @jtramm!