Adaptive Volume Estimators

[jtramm]

Overview

In the random ray solver there are a variety of ways to numerically treat the volume term that shows up in the MOC equations. You can read more background info in the "Volume Dilemma" section of the theory docs.

The crux of the issue is that the most natural treatment (known as the "naive" volume estimator) generally guarantees positivity and stability and has a low variance, though at the expense of being a biased ratio estimator. As such, eigenvalue solves using this naive setting might experience hundreds or even thousands of pcm error in k-eff, which is unacceptably high. Early in the random ray development the "simulation averaged" estimator was developed that accumulates volume from rays run over the course of the entire simulation, rather than only the last batch as in the naive approach. The simulation averaged estimator is thus unbiased (as it is effectively no longer a ratio estimator) but comes at the cost of 1) higher variance, and 2) potential for instability and negative fluxes under certain conditions.

Subsequent work introduced the "hybrid" estimator (merged in #3060) that took a "best of both worlds" approach. Generally the stability/negativity problems with the simulation averaged estimator are worst when there is an external source present in a region with a very low Sigma_t. Thus, the hybrid estimator tries to use the unbiased simulation averaged estimator everywhere it can, but then demoted regions that had an external source present.

To date the hybrid estimator has been the default and works well in most cases. However, in working on the JET problem, it became apparent that some fixed source problems were still becoming unstable and producing significant numbers of source regions with negative (or extremely large) fluxes. As such, JET simulations needed to be performed with the naive estimator. Given that the naive estimator is biased, I wondered if we might be able to make more efficient weight windows if we were able to come up with a better volume estimator that allowed more cells to use the simulation averaged estimator while ensuring stability.

This PR makes more progress to a "best of both worlds" approach between the naive and simulation averaged estimators by introducing more intelligent criteria for when to demote regions to naive that strike at the root issues. It adds two new volume estimators to the random ray solver (adaptive and strict_adaptive), and also introduces an auto setting that allows for OpenMC to make an intelligent choice about which estimator to use given the type of problem being simulated.

As a result of this PR, the user can now simply select "auto" for this setting (or leave it at default), and shouldn't need to think about it anymore:

settings.random_ray['volume_estimator'] = 'auto'

The Adaptive Estimator

adaptive uses simulation_averaged by default and "demotes" individual source regions to use the naive estimator where simulation_averaged is unsafe. Regions are selected in three ways:

1. Hit-starved regions (averaging fewer than 1.5 ray crossings per batch) receive the naive treatment in every batch, as under hybrid.
2. While the source is converging (the inactive batches), a per-batch guard applies the naive treatment to any region whose ratio of reduced source to flux estimate exceeds a fixed threshold (4×), or whose reduced source is negative while its flux is not (possible only with transport-corrected cross sections). This guard is what keeps the iteration finite on problems where simulation_averaged diverges.
3. From each region's running accumulated flux, first at the end of the inactive batches and then re-evaluated every active batch as the accumulation grows, two demote-only decisions are made. Regions whose accumulated flux is negative are demoted to the naive treatment (in the active phase this also watches the active-only accumulation that tally means are computed from), as are regions whose reduced source is dominated by contributions that do not derive from their own flux (scattering from other groups, fission, or an external source) in excess of 4x their accumulated flux. Once demoted, a region is never returned.

The division of labor matters. Per-batch tests act on noisy single-batch values, so they are confined to the inactive phase, where they serve only stability. In the active batches, the ones tallies accumulate, the treatment of every region changes only by demote-only decisions made from accumulated estimates, so the choice never churns with the sign or size of single-batch noise. Continuing the accumulated decisions through the active phase also matters at scale. On problems too large to afford an inactive phase that converges every deep region (verified on a 105M-region shielding model), the transition-time decisions alone miss regions whose instability only surfaces after tallies begin, and a single unprotected region can corrupt the solution far beyond its own boundary through scattering feedback.

The Strict Adaptive Estimator

adaptive selects estimators but never modifies a computed flux value. That is what keeps it unbiased, yet in some cases very slightly negative fluxes can still occur. Adaptive is configured to ensure simulation averaged is used maximally while preserving stability and still ensuring minimal bias from use of the naive estimator. This ensures good results on almost all problems I tested. However, on weight window generation problems that computed an adjoint source based on the forward flux, I did notice some problems popping up still. On the 105M source region JET model, a 0.67% negative fraction in the forward flux grew to ~3% in the adjoint and produced visibly noisy weight windows. This is a major improvement over 'hybrid', but is worse in practice than windows generated from the crude naive estimator. As my goal for this project was to improve performance over naive, this was quite disappointing. Thus, I developed the "strict adaptive estimator"

Fixing the remaining rare negative source regions via more demotion cannot close this gap. Negativity is contagious (a region on a safe estimator can still inherit it through in-scatter from neighbors that are not), so no selection rule alone controls the sign. strict_adaptive therefore runs the adaptive machinery unchanged and adds a per-batch fixup applied to the flux values themselves:

1. Rescue. A region whose batch flux comes out negative is first recomputed with that batch's own volume (algebraically, the naive update), whose term-by-term consistency removes the volume-mismatch noise responsible for most negative excursions.
2. Floor. If the flux is still negative, it is replaced with the previous iterate, keeping the flat flux moments non-negative batch over batch. (In every test in this PR, including the 105M-region model below, the resulting solutions contain no negative fluxes. With a linear source shape the reconstructed in-region flux can still locally dip below zero, so this is enforced behavior on the flux iterates rather than a theoretical guarantee on all outputs.)
3. Chronic demotion. A region that keeps going negative is demoted outright to the naive treatment. This channel is essential for bias rather than coverage. The floor prevents the accumulated flux from ever showing the negative sign that adaptive's own demotion watches for, so without an exit a chronically noisy region would be clipped upward every batch, biasing its flux. Demotion instead moves such regions onto an estimator that needs no clipping.

The auto default

The default volume_estimator is now auto, resolved once at the start of the run. Solves whose results feed variance reduction (weight window generation, and any adjoint workflow, including the forward solve an adjoint source is derived from, since that is where the negativity originates) receive strict_adaptive, and everything else receives adaptive. The report prints the resolution (e.g. Volume Estimator Type = Strict Adaptive (auto)), and setting any concrete estimator disables the routing.

Validation

Five problems, each comparing several estimators against a converged multigroup Monte Carlo or published reference. Every table reports the count of negative flux tally bins. The "miss rate" is the average fraction of source regions not crossed by any ray in a batch. OpenMC recommends staying below 1%, and several tests deliberately run far above it as a stress test. (A streaming-dominated problem, a one-group point source in a near-void room where simulation_averaged's failure is driven by ray-count variance rather than by the source-to-flux ratio, was also checked. adaptive matches hybrid to within 0.25% on all region tallies there while leaving 2 negative tally bins to hybrid's 26, and simulation_averaged fails catastrophically.)

1. C5G7 eigenvalue benchmark: no change where the default was already good

Pin-resolved 2D C5G7 (143k source regions), scored against the published eigenvalue (1.18655) and pin powers (AAPE = average absolute pin power error, MAX = worst pin):

miss rate	estimator	k dev (pcm)	AAPE	MAX pin	neg bins
1.0%	naive	−3,118	1.58%	7.32%	0
1.0%	simulation_averaged	+40	0.61%	2.14%	0
1.0%	hybrid	+40	0.61%	2.14%	0
1.0%	adaptive	+30	0.37%	1.66%	0
1.0%	strict_adaptive	+365	0.79%	4.36%	0
10.2%	naive	−5,662	1.29%	6.33%	0
10.2%	simulation_averaged	+61	0.67%	2.96%	0
10.2%	hybrid	+61	0.67%	2.96%	0
10.2%	adaptive	−1	0.58%	2.80%	0
10.2%	strict_adaptive	−671	4.65%	11.03%	0
20.4%	naive	−6,451	1.68%	7.53%	0
20.4%	simulation_averaged	+92	0.67%	3.04%	0
20.4%	hybrid	−79	2.15%	11.47%	0
20.4%	adaptive	−20	1.03%	6.95%	0
20.4%	strict_adaptive	−2,924	4.91%	10.30%	0

adaptive matches the accurate estimators from well-resolved conditions through heavy starvation, while naive shows the bias that motivates not simply using it everywhere. The strict_adaptive rows show the cost of its repair machinery growing with starvation, which is part of why it is not used for standard solves.

2. Shielding problem with scatter-fed regions: the failure the default has

The "irradiation vault" problem consists of a 12 cm cavity emitting 2 MeV neutrons at one side of a large near-void air hall (Σt ≈ 3×10⁻³ cm⁻¹ fast, 3×10⁻⁵ cm⁻¹ slow), several absorber blocks, and a thick concrete shell, discretized on a 1.3 cm source-region mesh (1.24M regions). The hall's slow group is fed almost entirely by down-scatter from the fast group, so its reduced source exceeds its flux by orders of magnitude with no external source anywhere, which is exactly the configuration hybrid's trigger cannot see. Slow-group fluxes at 12,800 rays (2.8% miss) vs multigroup Monte Carlo:

estimator	total	cavity	hall	shield	neg bins (of 85,750)
naive	+486% ± 0.3%	−0.1%	+508%	+469%	0
simulation_averaged	crashes (flux not finite)	—	—	—	—
hybrid	−22% ± 32%	−0.4%	−37%	+117%	5,468
adaptive	−0.5% ± 0.1%	−0.3%	−0.6%	+0.1%	0
strict_adaptive	0.0% ± 0.1%	−0.2%	−0.2%	+0.9%	0

Geometry and converged adaptive-estimator flux (z = 0 slice, log scale). The fast group beams from the source cavity and is shadowed by the absorber blocks, while the slow group is the scatter-fed distribution across the hall.

3. Published stress test: the three-region cube

The cube problem from here, as built by openmc.examples.random_ray_three_region_cube with the source region at the severe end of the paper's cross-section sweep, discretized on a 0.6 cm source-region mesh (216k regions). Region-integrated fluxes vs multigroup Monte Carlo at two ray densities:

miss rate	estimator	source	void	absorber	neg bins
18.4%	naive	−1.13%	+1.26%	−8.67%	0
18.4%	simulation_averaged	crashes (flux not finite)	—	—	—
18.4%	hybrid	−1.12%	−0.89%	−2.71%	0
18.4%	adaptive	−1.09%	+1.03%	+2.88%	0
18.4%	strict_adaptive	−1.01%	+2.36%	+2.78%	0
1.7%	naive	−0.94%	+0.48%	+0.93%	0
1.7%	simulation_averaged	crashes (flux not finite)	—	—	—
1.7%	hybrid	−0.95%	−0.41%	−0.20%	0
1.7%	adaptive	−0.96%	−0.29%	−0.04%	0
1.7%	strict_adaptive	−0.94%	+0.55%	+1.46%	0

4. Variations of the hall

This test isolates the case that motivated the strict adaptive estimator. It uses small variants of the irradiation vault from problem 2 (35^3 cells rather than 1.24M) run on a short schedule of 100 inactive and 100 active batches, which reproduces at small scale the conditions under which adaptive leaves a few negative bins on large problems. Two versions of the hall are used. The thin hall keeps the near-void air of the original problem, while the collisional variant raises the slow-group cross section of the hall gas by two orders of magnitude (keeping the scattering ratio near one), so the slow flux in the hall is built by many small collisions rather than by streaming. The collisional variant is also run at a quarter of the ray density (25% miss vs 4.3%) as a stress case. Slow-group region errors vs multigroup Monte Carlo, with negative bins counted over the 85,750 mesh tally bins:

configuration	estimator	cavity	hall	shield	neg bins
thin hall (Σt = 3×10⁻⁵), 2000 rays	naive	−4.9%	+63%	+874%	0
	adaptive	−4.9%	−5.3%	−2.5%	0
	strict_adaptive	−4.8%	−4.8%	−1.8%	0
collisional hall (Σt = 3×10⁻³), 2000 rays	naive	−4.9%	−4.4%	+7.5%	0
	adaptive	−4.9%	−5.3%	−2.6%	81
	strict_adaptive	−4.8%	−4.7%	−1.8%	0
collisional hall, 500 rays	naive	−3.1%	+134%	+66%	0
	adaptive	−5.9%	−4.6%	+3.0%	531
	strict_adaptive	−5.4%	+1.7%	+6.8%	0

strict_adaptive matches adaptive's accuracy where both are stable, produces zero negative fluxes where adaptive leaves dozens to hundreds, and degrades gracefully where naive blows up. The cost of the one-sided fixup is a conservative bias. On the 1%-miss C5G7 benchmark above, k rises to +365 pcm (vs +30 pcm for adaptive) and pin-power AAPE to 0.79% (vs 0.37%), which is why it is not used as the standard default. For weight window generation, there is an extreme preference for guarantees on positivity

5. JET

This was the original test case that motivated this project. The hybrid estimator resulted in instability. The naive estimator worked but I had wondered if we were leaving WW efficiency on the table. I compared the two estimators on a full CADIS test of the JET model, where weight windows are generated with the given estimator and then an MC solve is run. The maximum relative standard deviation across the 20 TLD tallies is then reported in each run. The baseline naive run resulted in a maximum rel std dev of 4.54%. The new default auto setting (which routes to the strict_adaptive estimator when WW's are being generated) resulted in a max rel std dev of 3.24%, both cases using essentially the same exact runtime. This PR is therefore responsible for a nearly 2x improvement in overall MC FOM for the JET test case

Compatibility and testing

• The default volume_estimator changes from hybrid to auto; hybrid, naive, and simulation_averaged are unchanged and remain selectable, and any explicit selection disables the automatic routing.
• Existing random-ray regression tests that relied on the implicit default are pinned to hybrid, so their reference results are unchanged. The adaptive and strict_adaptive estimators are covered by the volume-estimator regression tests (flat and linear) and by two deliberately ray-starved test variants that are sensitive to the demotion decisions. A new three-case test pins the auto routing itself, with forward, adjoint, and weight-window-generator runs that set no estimator, so a regression in either routing trigger shifts a reference result.
• At the default verbosity the end-of-run report adds a single line, the number of regions receiving the naive treatment in the final batch, with a per-cause breakdown (including the strict rescue/floor/chronic counts) available at verbosity 8. This may be helpful for understanding any problems that may arise with the new estimators in the future.
• The methods and user's guide documentation describe both estimators, the enforcement mechanics, and the auto default.

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)

[paulromano]

Thanks @jtramm! Doing a review in full right now but here is some immediate feedback from GPT-6-Astra:

Detailed findings

Blocking issues

[Major] Apply diagonal stabilization before flooring negative flux

In src/random_ray/flat_source_domain.cpp, the strict estimator replaces a negative raw flux with the previous iterate before apply_transport_stabilization() runs.

With transport-corrected cross sections, a negative diagonal (within-group) scattering term can produce a negative intermediate source-iteration value. This is expected: OpenMC's existing diagonal stabilization transforms that raw value into the appropriate stabilized value. Flooring the value first prevents stabilization from seeing the intermediate result and can freeze the iteration at an incorrect solution.

For example, if

$$ D = -2,\qquad \phi^{(n)} = 1,\qquad \phi_{\mathrm{raw}}^{(n+1)} = -1, $$

the existing stabilization gives

$$ \phi^{(n+1)} = \frac{\phi_{\mathrm{raw}}^{(n+1)}-D\phi^{(n)}}{1-D} = \frac{-1-(-2)(1)}{3} = \frac{1}{3}. $$

Strict adaptive instead floors the raw value to the previous iterate, producing

$$ \frac{1-(-2)(1)}{3}=1, $$

so the iterate does not move.

I confirmed this with a homogeneous two-group eigenvalue problem using a negative diagonal scattering term:

• Expected k-effective: $1/6$
• adaptive: 0.1666666642
• strict_adaptive: 2.5, unchanged through 200 batches

The relevant ordering is:

1. add_source_to_scalar_flux() performs the strict rescue/floor.
2. apply_transport_stabilization() runs afterward.

The positivity fix should therefore be applied after diagonal stabilization. If the naive-volume rescue is retained, the rescued candidate should also receive the corresponding stabilization before its sign is assessed.

The transport-corrected regression test currently exercises only adaptive. Please add a strict_adaptive case that verifies convergence to the expected result.

Non-blocking suggestions

[Moderate] Use the void-specific additive term when rescuing void flux

The rescue calculation assumes the non-void update

$$ \phi = \frac{T}{V} + \frac{q}{\Sigma_t} $$

and reconstructs the transport contribution by subtracting source(sr, g).

True void regions use a different update:

$$ \phi = \frac{T}{V} + \frac{1}{2}q_{\mathrm{ext}},\overline{\ell^2}. $$

Consequently, subtracting and restoring source(sr, g) does not algebraically reproduce the naive-volume update for a sourced void region.

Using the implemented routines with the same transport contribution, I obtained:

• Naive-volume update: 4.5
• Strict rescue: 0

The rescue should preserve the original transport contribution directly or subtract the actual additive term used by the void update.

[jtramm]

Thanks @paulromano! I agree with both issues and was able to reproduce the diagonal stabilization issue on my end. I've added fixes for both.

[paulromano]

Good to go now; I made one small update to avoid allocations when they weren't needed but everything else looked fine. Thanks for the improvement @jtramm!