Random Ray Linear Source Gradient Limiter

[jtramm]

Overview

This PR adds an optional limiter for linear source gradients. It is off by default, but can be enabled via:

settings.random_ray['source_gradient_limiter'] = True

Background

The fitted gradient inverts the region's spatial moment matrix, which amplifies noise along any thin extent of the region, so a poorly sampled region can carry a spuriously steep gradient. Such a region emits a negative source over part of its extent, and rays crossing the negative lobe carry negative angular flux downstream.

When enabled, the limiter rescales the gradient so the modeled source stays non-negative over the region's bounding box, preserving mean emission exactly. Linear source mode only, flat source solves are unaffected.

A similar gradient limiting treatment exists in production MOC as MPACT's limited linear source approximation, developed by Sooyoung Choi. MPACT computes each region's minimum source exactly, as the smaller of the modeled source at the entrance and exit of every segment crossing the region, over all rays and angles, and when that minimum is negative it rescales the gradient by the mean-preserving factor q/(q − q_min).

The reduction factor being introduced here follows similar logic. What differs is how the minimum is found. MPACT's exact minimum needs the whole track set before the sweep, which deterministic MOC has because its rays are laid down once and reused every iteration. Random ray samples new rays every batch, so no segment set exists when the source is built. Instead, this PR bounds the modeled source over the region's axis-aligned bounding box. Ray segment endpoints lie on the region boundary (or inside it, where a ray starts or ends), so accumulating their componentwise minimum and maximum over the batches gives a box that, once the region's extreme points along each axis have been hit, contains the region. The lowest value the linear term can take over the box is at a corner, a sum computed from the gradient components and the centroid's offsets to the box faces. When the flat source does not cover that dip, the gradient is scaled by their ratio, which is MPACT's factor with the box minimum standing in for the segment minimum.

The box only grows, since the geometry is static, and after the first few batches it covers the region. The bound is exact for axis-aligned box regions and conservative for others: a sphere is limited by up to sqrt(3) more than necessary, and a thin region lying diagonally to the axes by much more, since its bounding box is far larger than the region. It costs two positions per source region, accumulated only when the limiter is enabled.

The limiter is off by default because limiting cuts both ways. A steep fit can be noise, but it can also be physical. The clearest case is deep penetration, where a linear fit to an exponentially attenuated flux legitimately crosses zero near the far edge of a source region, and clipping it produces an error that compounds with depth. The results below quantify both sides.

I experimented with a number of methods for attempting to dynamically limit gradients in specific source regions when certain problematic conditions were detected. However, I was generally not able to find a strategy that worked in all types of cases, so am leaving that as a potential for future work. Users should therefore enable the limiter when they observe negative sources destabilizing a linear source simulation.

This PR is part of a set of recent numerical improvements to the random ray solver, e.g., the adaptive volume estimator PR (#4110), and the linear source consistency PR (#4111). These PRs are reviewable in any order. All results below are measured against a "baseline" of the combination of #4110 and #4111 together with the limiter disabled, so the deltas shown are what enabling the limiter does on top of everything else. Both arms of every comparison run the same binary and differ only in the setting.

Results

Stability was measured across a range of ray densities, three estimators, and four problems, all with linear sources. The test problems:

• Thin hall: the irradiation vault hall of Adaptive Volume Estimators #4110 (a near-void air hall fed by down-scatter, with slow-group scattering ratios near one). Negativity is counted in slow-group bins, as that is the group that amplifies it. Accuracy is the fast-group total flux error vs Monte Carlo, since at high miss rates the slow group's percentage errors are inflated by the near-unity scattering ratio for every method alike.
• Collisional hall: the Adaptive Volume Estimators #4110 variant with the hall cross section raised by two orders of magnitude, so the slow flux is built by many small collisions rather than by streaming. Accuracy is the slow-group total flux error vs Monte Carlo.
• TLD: a point source at the center of a large near-void room wrapped in a thin wall, in a single energy group. Distant cells accumulate their flux from a few highly variable ray crossings per batch, so their fitted gradients are noise dominated. Accuracy is the total flux error vs Monte Carlo.
• Three-region cube: the published source/void/absorber stress problem used in Adaptive Volume Estimators #4110 and Naive Flux Linear Source Consistency Fix #4111. Accuracy is the total flux error vs Monte Carlo.

Thin hall with the naive estimator, varying the ray density (gradient limiter off -> on):

rays	miss rate	neg slow bins	most negative bin	fast flux err
500	25%	668 -> 0	-4.8e-2 -> none	-9.9% -> -10.6%
1,000	9.5%	340 -> 0	-1.1e-2 -> none	-3.0% -> -4.6%
2,000	4.3%	25 -> 0	-2.7e-3 -> none	-1.4% -> -2.5%
4,000	2.8%	2 -> 0	-2.0e-4 -> none	+1.0% -> -0.3%

As we can see above, limiting gradients does greatly reduce negative flux occurrences. However, at lower ray densities, this can also come with a small loss of accuracy.

Now we look at a wider variety of problems and volume estimators:

case	estimator	neg bins	flux err
thin hall	auto	2 -> 0	-1.3% -> -2.5%
thin hall	hybrid	5,645 -> 5,598	-1.2% -> -2.4%
collisional hall	auto	120 -> 80	-2.3% -> -4.1%
collisional hall, 500 rays	auto	576 -> 454	-3.4% -> -3.2%
TLD	naive	3,250 -> 0	-5.8% -> -7.1%
TLD	auto	1,474 -> 456	-6.0% -> -7.1%
three-region cube	auto	0 -> 0	-0.4% -> -0.4%

The limiter removes all of the negativity in the problems where noisy gradients cause it. The TLD naive case drops from 3,250 negative bins to zero, and the thin hall reaches zero at every density and seed. The hybrid row is unaffected because its instability comes from the volume estimator issue addressed in #4110 rather than from source gradients. The collisional hall and the TLD auto case retain a residual because part of their negativity arises in the flat source and flux themselves, which the limiter leaves alone by design.

The flux error column shows the cost. Enabling the limiter shifts the flux in these problems by up to about two percentage points, nearly always in the negative direction. This is the clipping bias. The walls of these problems carry physically steep gradients, and limiting them reduces transmission. In problems where steep gradients dominate the cost is larger, as the slab results below show. This is the tradeoff that makes the limiter an option rather than a default.

The cost where linear sources are working well

Coarse-mesh C5G7 (102x102 overlay, ~15k source regions) is the case where linear sources deliver a large accuracy win over flat. Here we vary the ray density and observe the eigenvalue and pin power (AAPE) errors.

rays	miss rate	flat AAPE	LS AAPE, off	LS AAPE, on	k dev, off	k dev, on
200	0.16%	1.51%	0.65%	0.86%	+45 pcm	+163 pcm
100	3.7%	1.82%	1.23%	1.43%	-131	+201
70	9.8%	1.76%	2.22%	2.05%	-2,341	-590
50	19%	2.36%	2.92%	2.78%	-3,721	-1,114

At fine ray counts the limiter costs about 0.2 AAPE and about 100 pcm, and limited linear sources still beat flat (by 1.75x at 200 rays and 1.3x at 100). Under heavy ray starvation it removes most of the large k bias that noisy gradients otherwise develop and pulls the pin powers back toward the flat-source level.

Runtime

The limiter's cost is the bounding-box update, two endpoint comparisons per ray segment inside the source region lock, plus a few operations per region and group in the source update. Measured on the 200-ray C5G7 case above (3,000 batches, 4 threads, three repeats per arm, run back to back on an otherwise idle machine):

	limiter off	limiter on	change
total simulation time	70.3 s (69.7--70.7)	73.4 s (73.1--74.0)	+4.4%
transport sweep	52.2 s (51.4--52.6)	54.8 s (54.4--55.2)	+5.0%

MPACT reports 5--6% for its limited linear source approximation, most of it in the segment search. Since the box is not needed when the limiter is off, the bookkeeping exists only when it is on, and flat source and unlimited linear source runs carry no cost. Enabling the limiter on a well-resolved problem is mildly harmful and unnecessary, which is the fine-ray half of the default-off rationale.

Where steep gradients are physical

The other half of the default-off rationale is that whole problem classes offer the limiter nothing to fix. The test problem is a slab 24 mean free paths deep, with a source at one end, mild scattering (c = 0.3), and source regions 3 mean free paths thick. The reference is a discretization-converged fine-mesh flat-source run (0.1 MFP regions). Flux errors at depth:

depth (MFP)	flat sources	linear, off	linear, on
3--6	+32%	-4.7%	+4.6%
9--12	+174%	-16.0%	+18.0%
15--18	+350%	-19.4%	+25.3%
21--24	+594%	-21.9%	+31.8%

Linear sources are the whole game on this problem, cutting the deep-flux error from +594% to -22% at this coarse region size. Every fitted gradient here is physical. The flux decays exponentially, so the linear fit legitimately crosses zero near each region's far edge, and the limiter clips every one of them. The result is not a breakdown. Enabling the limiter trades the -22% under-prediction for a +32% over-prediction at every depth. Both errors shrink as regions are refined: at 1.5 mean free paths the deepest bin is -1.9% with or without the limiter, which by then fires only at the 1e-7 level, and at 1.0 it is -0.6% with the limiter fully quiet. But at the coarse size the limited answer comes from discarding correct shape information, so on a problem like this the limiter changes every deep result while having nothing to offer.

Compatibility and testing

• A new random_ray_linear_source_stability regression test pins the limiter on in a naive linear source configuration whose spherical region boundaries cut the overlay mesh into pieces with off-center centroids, the geometry where the bound's asymmetry matters.
• The methods documentation gives the bounding-box bound and the mean-emission-preserving rescale, and relates the bound to MPACT's exact segment minimum. The user's guide shows how to enable the limiter and when to want it.

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)