Analog delayed-nu-fission tallies with an EnergyoutFilter and a DelayedGroupFilter apply weighted filters twice #4181

Bug Description

When score_fission_eout() (src/tallies/tally_scoring.cpp) scores a delayed fission neutron on a tally that has both an EnergyoutFilter and a DelayedGroupFilter, it applies the filter weights twice. At a5bc348 it multiplies the score by the product of the filter weights (lines 482-489) and passes score * filter_weight to score_fission_delayed_dg() (lines 491-492), which multiplies by the same product again (lines 162-174). A bin that should score w therefore scores w^2: for example, an EnergyFunctionFilter with a constant value of 2 multiplies such a tally by 4 instead of 2.

Affected are delayed-nu-fission tallies (an EnergyoutFilter always selects the analog estimator) that combine an EnergyoutFilter, a DelayedGroupFilter and a filter with non-unit weights, such as EnergyFunctionFilter (continuous energy only) or one of the expansion filters (LegendreFilter, SphericalHarmonicsFilter, SpatialLegendreFilter, ZernikeFilter, ZernikeRadialFilter). Both continuous-energy and multigroup runs take this path. Not affected: tallies whose filters all have unit weights (cell, material, mesh, energy and so on), nu-fission and prompt-nu-fission with an EnergyoutFilter, and delayed-nu-fission without an EnergyoutFilter or without a DelayedGroupFilter.

The call site's multiplication was correct until d5a2340 ("Consistent use of filter_weight", #1462), which moved the weight product into score_fission_delayed_dg() but left it at this call site as well. Every release since v0.12.0 has the double weighting.

I will open a pull request with a fix and a regression test.

Steps to Reproduce

The script below needs no nuclear data: it writes its own two-group library with two delayed groups and runs the multigroup slab of openmc.examples.slab_mg() (0 <= x <= 929.45 cm, 1000 particles, 10 batches, 5 inactive, seed 1). It tallies analog delayed-nu-fission three ways, [EnergyoutFilter, DelayedGroupFilter, SpatialLegendreFilter], [DelayedGroupFilter, SpatialLegendreFilter] and [EnergyoutFilter, SpatialLegendreFilter], so all three tallies score the same delayed neutrons. Summed over its outgoing-energy bins, the first tally must equal the second; summed over its delayed groups, it must equal the third. The P0 weight is 1, so the P0 bins agree whatever the bug does. The P1 weight is x/929.45, between 0 and 1 in the slab, so the P1 bins of the first tally come out too small when the weight is applied twice. Run it in an empty directory with openmc on the PATH:

"""Delayed-nu-fission tallied with an EnergyoutFilter, a DelayedGroupFilter and
a SpatialLegendreFilter: summing over the outgoing-energy bins should give the
same tally without the EnergyoutFilter. Multigroup, no nuclear data needed."""
import numpy as np
import openmc
from openmc.examples import slab_mg

# Two-group library with two delayed groups
groups = openmc.mgxs.EnergyGroups(group_edges=[0.0, 0.625, 20.0e6])
beta = np.array([0.1, 0.2])
fission = np.array([0.002817, 0.097])
nu_fission = 2.5 * fission
capture = np.array([0.008708, 0.02518])
scatter = np.array(
    [[[0.31980, 0.06694], [0.004555, -0.0003972]],
     [[0.00000, 0.00000], [0.424100, 0.05439000]]])
xsdata = openmc.XSdata('mat_1', groups, num_delayed_groups=2)
xsdata.order = 1
xsdata.set_total([0.33588, 0.54628])
xsdata.set_absorption(capture + fission)
xsdata.set_scatter_matrix(scatter)
xsdata.set_fission(fission)
xsdata.set_prompt_nu_fission((1.0 - beta.sum()) * nu_fission)
xsdata.set_delayed_nu_fission(np.outer(beta, nu_fission))
xsdata.set_chi_prompt([1.0, 0.0])
xsdata.set_chi_delayed([[0.8, 0.2], [0.6, 0.4]])
xsdata.set_decay_rate([0.0127, 0.0317])
library = openmc.MGXSLibrary(groups, num_delayed_groups=2)
library.add_xsdata(xsdata)
library.export_to_hdf5('2g.h5')

# Slab 0 <= x <= 929.45 cm; 1000 particles, 10 batches (5 inactive)
model = slab_mg()
model.settings.seed = 1

# P1 weight x/929.45 lies in [0, 1] inside the slab
legendre = openmc.SpatialLegendreFilter(1, 'x', -929.45, 929.45)
energyout = openmc.EnergyoutFilter([0.0, 0.625, 20.0e6])
delayed_groups = openmc.DelayedGroupFilter([1, 2])
for name, filters in [
    ('energyout-delayedgroup', [energyout, delayed_groups, legendre]),
    ('delayedgroup', [delayed_groups, legendre]),
    ('energyout', [energyout, legendre]),
]:
    tally = openmc.Tally(name=name)
    tally.filters = filters
    tally.scores = ['delayed-nu-fission']
    tally.estimator = 'analog'
    model.tallies.append(tally)

sp_path = model.run()
with openmc.StatePoint(sp_path) as sp:
    a = sp.get_tally(name='energyout-delayedgroup').get_reshaped_data(value='mean')
    b = sp.get_tally(name='delayedgroup').get_reshaped_data(value='mean')
    c = sp.get_tally(name='energyout').get_reshaped_data(value='mean')
# a[energyout, delayed group, Legendre order], b[delayed group, Legendre order],
# c[energyout, Legendre order]
ratio = a[..., 0, 0].sum(axis=0) / b[..., 0, 0]
print('sum over energyout bins / tally without EnergyoutFilter')
print('(rows: delayed groups 1, 2; columns: P0, P1):')
print(ratio)
print(f'max |ratio - 1| over the P1 bins: {np.abs(ratio[:, 1] - 1.0).max():.3e}')
ratio_dg = a[..., 0, 0].sum(axis=1) / c[..., 0, 0]
print('sum over delayed groups / tally without DelayedGroupFilter')
print('(rows: outgoing energy groups; columns: P0, P1):')
print(ratio_dg)

Output on develop (a5bc348), after OpenMC's own output:

sum over energyout bins / tally without EnergyoutFilter
(rows: delayed groups 1, 2; columns: P0, P1):
[[1.         0.70032602]
 [1.         0.68763328]]
max |ratio - 1| over the P1 bins: 3.124e-01
sum over delayed groups / tally without DelayedGroupFilter
(rows: outgoing energy groups; columns: P0, P1):
[[1.         0.70122105]
 [1.         0.68720811]]

Expected: every ratio is 1 to round-off. With the multiplication at the call site removed, the same script prints 1 for every ratio, and a P1 deviation at round-off (about 1e-16).

Environment

• OpenMC develop at a5bc348 (openmc --version: 0.16.1-dev57), built from source with CMake: RelWithDebInfo, -DOPENMC_ENABLE_STRICT_FP=on, OpenMP on, MPI off, HDF5 1.10.10 (serial).
• Ubuntu 24.04 (x86_64), GCC 13.3.0, Python 3.12.3 with numpy 2.5.3 and h5py 3.16.0; OMP_NUM_THREADS=2.
• Nuclear data: none; the script writes its own multigroup library.