Fix double filter weighting of delayed-nu-fission with EnergyoutFilter and DelayedGroupFilter

src/tallies/tally_scoring.cpp

@@ -476,20 +476,13 @@ void score_fission_eout(Particle& p, int i_tally, int i_score, int score_bin)
         const DelayedGroupFilter& dg_filt {*dynamic_cast<DelayedGroupFilter*>(
           model::tally_filters[i_dg_filt].get())};
 
-        // Loop over delayed group bins until the corresponding bin is found
+        // Loop over delayed group bins until the corresponding bin is found.
+        // The helper applies the filter weights itself, so the score passed
+        // to it must not already include them.
         for (auto d_bin = 0; d_bin < dg_filt.n_bins(); ++d_bin) {
           if (dg_filt.groups()[d_bin] == g) {
-            // Find the filter index and weight for this filter combination
-            double filter_weight = 1.;
-            for (auto j = 0; j < tally.filters().size(); ++j) {
-              auto i_filt = tally.filters(j);
-              auto& match {p.filter_matches(i_filt)};
-              auto i_bin = match.i_bin_;
-              filter_weight *= match.weights_[i_bin];
-            }
-
-            score_fission_delayed_dg(i_tally, d_bin, score * filter_weight,
-              i_score, p.filter_matches());
+            score_fission_delayed_dg(
+              i_tally, d_bin, score, i_score, p.filter_matches());
           }
         }
 

tests/regression_tests/mg_delayed_eout_weights/inputs_true.dat

@@ -0,0 +1,62 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <cross_sections>2g.h5</cross_sections>
+    <material id="1" name="mat_1">
+      <density value="1.0" units="macro"/>
+      <macroscopic name="mat_1"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="1 -2" universe="0"/>
+    <surface id="1" type="x-plane" boundary="reflective" coeffs="0.0"/>
+    <surface id="2" type="x-plane" boundary="vacuum" coeffs="929.45"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>1000</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <source type="independent" strength="1.0" particle="neutron">
+      <space type="box">
+        <parameters>0.0 -1000.0 -1000.0 929.45 1000.0 1000.0</parameters>
+      </space>
+    </source>
+    <output>
+      <summary>false</summary>
+    </output>
+    <energy_mode>multi-group</energy_mode>
+    <tabular_legendre>
+      <enable>false</enable>
+    </tabular_legendre>
+  </settings>
+  <tallies>
+    <filter id="2" type="energyout">
+      <bins>0.0 0.625 20000000.0</bins>
+    </filter>
+    <filter id="3" type="delayedgroup">
+      <bins>1 2</bins>
+    </filter>
+    <filter id="1" type="spatiallegendre">
+      <order>1</order>
+      <axis>x</axis>
+      <min>-929.45</min>
+      <max>929.45</max>
+    </filter>
+    <tally id="1" name="energyout-delayedgroup">
+      <filters>2 3 1</filters>
+      <scores>delayed-nu-fission</scores>
+      <estimator>analog</estimator>
+    </tally>
+    <tally id="2" name="delayedgroup">
+      <filters>3 1</filters>
+      <scores>delayed-nu-fission</scores>
+      <estimator>analog</estimator>
+    </tally>
+    <tally id="3" name="energyout">
+      <filters>2 1</filters>
+      <scores>delayed-nu-fission</scores>
+      <estimator>analog</estimator>
+    </tally>
+  </tallies>
+</model>

tests/regression_tests/mg_delayed_eout_weights/results_true.dat

@@ -0,0 +1,37 @@
+k-combined:
+1.117031E+00 1.558955E-02
+tally 1:
+1.309236E-01
+3.592460E-03
+6.736655E-02
+9.501908E-04
+4.519839E-01
+4.128356E-02
+2.120194E-01
+9.242207E-03
+4.456540E-01
+3.981255E-02
+2.141939E-01
+9.262961E-03
+6.551057E-01
+8.622179E-02
+3.090335E-01
+1.917417E-02
+tally 2:
+5.765776E-01
+6.658661E-02
+2.815605E-01
+1.597098E-02
+1.107090E+00
+2.465786E-01
+5.210529E-01
+5.486740E-02
+tally 3:
+5.829075E-01
+6.879483E-02
+2.793860E-01
+1.598254E-02
+1.100760E+00
+2.428386E-01
+5.232274E-01
+5.502966E-02

tests/regression_tests/mg_delayed_eout_weights/test.py

@@ -0,0 +1,109 @@
+"""Filter weights on delayed-nu-fission tallies with an outgoing energy filter.
+
+With an EnergyoutFilter, delayed-nu-fission is scored once per banked fission
+neutron, and a DelayedGroupFilter sends each delayed neutron to its group bin.
+The weight of any other filter on the tally, here a SpatialLegendreFilter,
+must be applied exactly once. The same delayed neutrons are also tallied
+without the EnergyoutFilter and without the DelayedGroupFilter, which are
+scored along other paths, so summing the first tally over either of those
+filters must reproduce the corresponding tally to round-off.
+
+"""
+
+import os
+
+import numpy as np
+import openmc
+from openmc.examples import slab_mg
+
+from tests.testing_harness import PyAPITestHarness
+
+
+def create_library():
+    # Two-group cross sections from the mg_tallies test with two delayed
+    # groups added. The delayed neutron fractions are unphysically large so
+    # that a short run banks many delayed neutrons, and delayed neutrons are
+    # born in both groups so that every outgoing energy bin is scored.
+    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')
+
+
+class MGDelayedEoutHarness(PyAPITestHarness):
+    def _compare_results(self):
+        # The tallies must agree with each other whatever the reference
+        # results are, so check that first
+        with openmc.StatePoint(self._sp_name) as sp:
+            eout_dg = sp.get_tally(name='energyout-delayedgroup')
+            dg = sp.get_tally(name='delayedgroup')
+            eout = sp.get_tally(name='energyout')
+            eout_dg = eout_dg.get_reshaped_data(value='mean')[..., 0, 0]
+            dg = dg.get_reshaped_data(value='mean')[..., 0, 0]
+            eout = eout.get_reshaped_data(value='mean')[..., 0, 0]
+
+        # Every bin must be scored for the comparisons to mean anything
+        assert np.all(dg > 0.0)
+        assert np.all(eout > 0.0)
+
+        # Indices of eout_dg are (energyout, delayed group, Legendre order)
+        np.testing.assert_allclose(eout_dg.sum(axis=0), dg, rtol=1e-10)
+        np.testing.assert_allclose(eout_dg.sum(axis=1), eout, rtol=1e-10)
+
+        super()._compare_results()
+
+    def _cleanup(self):
+        super()._cleanup()
+        f = '2g.h5'
+        if os.path.exists(f):
+            os.remove(f)
+
+
+def test_mg_delayed_eout_weights():
+    create_library()
+    model = slab_mg()
+
+    # The slab spans 0 <= x <= 929.45 cm. With the filter defined over twice
+    # that width, the P1 weight is x/929.45, which lies between 0 and 1 in the
+    # slab, so the P1 results are sums of non-negative terms that can be
+    # compared with a relative tolerance.
+    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']
+        # Required: the SpatialLegendreFilter would otherwise give the
+        # delayedgroup tally a collision estimator, and the checks rely on all
+        # three tallies scoring the same analog events
+        tally.estimator = 'analog'
+        model.tallies.append(tally)
+
+    harness = MGDelayedEoutHarness('statepoint.10.h5', model)
+    harness.main()