From cfca8bee2d272de55e707560d1d6fa1fa2b7419b Mon Sep 17 00:00:00 2001 From: SHI JINKANG <139621604+kongshifer@users.noreply.github.com> Date: Fri, 23 Oct 2026 18:40:22 +0800 Subject: [PATCH] Fix multigroup scattering weights for outgoing-group-dependent multiplicity Use an incoming-group multiplier with the production-weighted outgoing-group distribution so that the expected scattering production matches the input. Keep channel data and angular distributions intact, correct the total event cross-section query, and cover the change with deterministic and analytical tests. Refs #4205 --- docs/source/methods/cross_sections.rst | 38 +++- include/openmc/scattdata.h | 4 + src/scattdata.cpp | 20 +- tests/cpp_unit_tests/CMakeLists.txt | 1 + tests/cpp_unit_tests/test_scattdata.cpp | 187 ++++++++++++++++++ .../mg_scattering_multiplicity/__init__.py | 0 .../mg_scattering_multiplicity/test.py | 99 ++++++++++ 7 files changed, 337 insertions(+), 12 deletions(-) create mode 100644 tests/cpp_unit_tests/test_scattdata.cpp create mode 100644 tests/regression_tests/mg_scattering_multiplicity/__init__.py create mode 100644 tests/regression_tests/mg_scattering_multiplicity/test.py diff --git a/docs/source/methods/cross_sections.rst b/docs/source/methods/cross_sections.rst index 764c4c62808..d088a46bc76 100644 --- a/docs/source/methods/cross_sections.rst +++ b/docs/source/methods/cross_sections.rst @@ -236,19 +236,43 @@ Unlike the continuous-energy mode, the multi-group mode does not explicitly track particles produced from scattering multiplication (i.e., :math:`(n,xn)`) reactions. These are instead accounted for by adjusting the weight of the particle after the collision such that the correct total weight is maintained. -The weight adjustment factor is optionally provided by the `multiplicity` data -which is required to be provided in the form of a group-wise matrix. This data -is provided as a group-wise matrix since the probability of producing multiple -particles in a scattering reaction depends on both the incoming energy, `g`, and -the sampled outgoing energy, `g'`. This data represents the average number of -particles emitted from a scattering reaction, given a scattering reaction has -occurred: +The optional `multiplicity` data is provided as a group-wise matrix, since the +average number of particles produced in a scattering reaction can depend on both +the incoming energy group, `g`, and the outgoing energy group, `g'`. Each entry +is the ratio of the neutron-production cross section to the scattering-event +cross section for that channel: .. math:: multiplicity_{g \rightarrow g'} = \frac{\nu_{scatter}\sigma_{s,g \rightarrow g'}}{ \sigma_{s,g \rightarrow g'}} +When multiplicity data are provided, the scattering matrix contains neutron +production cross sections. Let :math:`\nu_{g \rightarrow h}` denote the +multiplicity for an incoming group `g` and outgoing group `h`. The outgoing group +is sampled from the normalized production row: + +.. math:: + + P(g' \mid g) = + \frac{\nu_{g \rightarrow g'}\sigma_{s,g \rightarrow g'}} + {\sum_h \nu_{g \rightarrow h}\sigma_{s,g \rightarrow h}}. + +The particle weight is therefore multiplied by the ratio of total neutron +production to total scattering events for the incoming group: + +.. math:: + + w' = w\frac{\sum_h \nu_{g \rightarrow h}\sigma_{s,g \rightarrow h}} + {\sum_h \sigma_{s,g \rightarrow h}}. + +For consistent total and absorption cross sections, the denominator equals +:math:`\sigma_{t,g} - \sigma_{a,g}`. This multiplier is independent of the sampled +outgoing group. Multiplying by that channel's multiplicity instead would count +its production bias a second time. If multiplicity is constant across the +nonzero channels in a row, the weight multiplier reduces to that constant. +The conditional angular distribution for each channel is unchanged. + If this scattering multiplication information is not provided in the library then no weight adjustment will be performed. This is equivalent to neglecting any additional particles produced in scattering multiplication reactions. diff --git a/include/openmc/scattdata.h b/include/openmc/scattdata.h index bea881402ac..21fe7f15ecb 100644 --- a/include/openmc/scattdata.h +++ b/include/openmc/scattdata.h @@ -25,6 +25,10 @@ class ScattData { virtual ~ScattData() = default; protected: + // B_g / C_g: average production per event when sampling the production + // matrix. + vector weight_mult_; + //! \brief Initializes the attributes of the base class. void base_init(int order, const tensor::Tensor& in_gmin, const tensor::Tensor& in_gmax, const double_2dvec& in_energy, diff --git a/src/scattdata.cpp b/src/scattdata.cpp index 7dd53817ea4..55b995d7730 100644 --- a/src/scattdata.cpp +++ b/src/scattdata.cpp @@ -29,6 +29,7 @@ void ScattData::base_init(int order, const tensor::Tensor& in_gmin, energy.resize(groups); mult.resize(groups); dist.resize(groups); + weight_mult_.assign(groups, 1.); for (int gin = 0; gin < groups; gin++) { // Store the inputted data @@ -46,6 +47,16 @@ void ScattData::base_init(int order, const tensor::Tensor& in_gmin, double norm = std::accumulate(energy[gin].begin(), energy[gin].end(), 0.); if (norm != 0.) { + // The outgoing group is sampled from neutron production, B_gh / B_g. + // Use the row multiplier B_g / C_g, not the sampled channel multiplicity, + // so the expected production remains C_g * (B_gh / B_g) * (B_g / C_g). + double scatter_norm = 0.; + for (int go = 0; go < energy[gin].size(); go++) { + if (energy[gin][go] != 0.) { + scatter_norm += energy[gin][go] / mult[gin][go]; + } + } + weight_mult_[gin] = norm / scatter_norm; for (auto& n : energy[gin]) n /= norm; } @@ -197,8 +208,7 @@ double ScattData::get_xs( if (gout != nullptr) { val *= energy[gin][i_gout] / mult[gin][i_gout]; } else { - val /= std::inner_product( - mult[gin].begin(), mult[gin].end(), energy[gin].begin(), 0.0); + val /= weight_mult_[gin]; } break; case MgxsType::NU_SCATTER_FMU: @@ -365,7 +375,7 @@ void ScattDataLegendre::sample( } // Update the weight to reflect neutron multiplicity - wgt *= mult[gin][i_gout]; + wgt *= weight_mult_[gin]; } //============================================================================== @@ -553,7 +563,7 @@ void ScattDataHistogram::sample( mu = std::clamp(mu, -1., 1.); // Update the weight to reflect neutron multiplicity - wgt *= mult[gin][i_gout]; + wgt *= weight_mult_[gin]; } //============================================================================== @@ -766,7 +776,7 @@ void ScattDataTabular::sample( mu = std::clamp(mu, -1., 1.); // Update the weight to reflect neutron multiplicity - wgt *= mult[gin][i_gout]; + wgt *= weight_mult_[gin]; } //============================================================================== diff --git a/tests/cpp_unit_tests/CMakeLists.txt b/tests/cpp_unit_tests/CMakeLists.txt index 82bcf319ca0..61fb22be9da 100644 --- a/tests/cpp_unit_tests/CMakeLists.txt +++ b/tests/cpp_unit_tests/CMakeLists.txt @@ -10,6 +10,7 @@ set(TEST_NAMES test_photon test_ray test_region + test_scattdata test_surface test_tensor test_geometry diff --git a/tests/cpp_unit_tests/test_scattdata.cpp b/tests/cpp_unit_tests/test_scattdata.cpp new file mode 100644 index 00000000000..283f9c4d357 --- /dev/null +++ b/tests/cpp_unit_tests/test_scattdata.cpp @@ -0,0 +1,187 @@ +#include "openmc/scattdata.h" + +#include +#include + +#include +#include +#include + +#include "openmc/random_lcg.h" + +namespace { + +using Catch::Matchers::WithinAbs; +using openmc::double_2dvec; +using openmc::double_3dvec; +using openmc::MgxsType; + +// Use the same isotropic production cross sections in all representations. +// With two histogram bins, each bin has width one. +template +T make_scatter(const double_2dvec& production, const double_2dvec& multiplicity) +{ + openmc::tensor::Tensor gmin({2}, 0); + openmc::tensor::Tensor gmax({2}, 1); + gmin[1] = 1; + + double_3dvec coefficients(2); + for (int gin = 0; gin < 2; ++gin) { + for (double value : production[gin]) { + if constexpr (std::is_same_v) { + coefficients[gin].push_back({value}); + } else { + coefficients[gin].push_back({value / 2.0, value / 2.0}); + } + } + } + + T data; + data.init(gmin, gmax, multiplicity, coefficients); + return data; +} + +void check_cross_sections(openmc::ScattData& data, + const double_2dvec& production, const double_2dvec& scattering) +{ + for (int gin = 0; gin < 2; ++gin) { + CAPTURE(gin); + double sum_production = 0.0; + double sum_scattering = 0.0; + for (int i = 0; i < production[gin].size(); ++i) { + sum_production += production[gin][i]; + sum_scattering += scattering[gin][i]; + int gout = data.gmin[gin] + i; + CHECK_THAT(data.get_xs(MgxsType::NU_SCATTER, gin, &gout, nullptr), + WithinAbs(production[gin][i], 1e-12)); + CHECK_THAT(data.get_xs(MgxsType::SCATTER, gin, &gout, nullptr), + WithinAbs(scattering[gin][i], 1e-12)); + } + CHECK_THAT(data.get_xs(MgxsType::NU_SCATTER, gin, nullptr, nullptr), + WithinAbs(sum_production, 1e-12)); + CHECK_THAT(data.get_xs(MgxsType::SCATTER, gin, nullptr, nullptr), + WithinAbs(sum_scattering, 1e-12)); + } +} + +void check_sample_weights( + openmc::ScattData& data, int gin, double expected_multiplicity) +{ + uint64_t seed = openmc::init_seed(1, 0); + constexpr double initial_weight = 2.5; + for (int i = 0; i < 32; ++i) { + CAPTURE(gin, i); + int gout; + double mu; + double weight = initial_weight; + data.sample(gin, gout, mu, weight, &seed); + REQUIRE(gout >= data.gmin[gin]); + REQUIRE(gout <= data.gmax[gin]); + REQUIRE(mu >= -1.0); + REQUIRE(mu <= 1.0); + // The weight multiplier is fixed by the incoming group, regardless of the + // outgoing group drawn from the production-weighted distribution. + REQUIRE_THAT( + weight, WithinAbs(initial_weight * expected_multiplicity, 1e-12)); + } +} + +} // namespace + +TEMPLATE_TEST_CASE( + "Multigroup scattering cross sections with variable multiplicity", + "[scattdata]", openmc::ScattDataLegendre, openmc::ScattDataHistogram, + openmc::ScattDataTabular) +{ + auto data = make_scatter({{0.2, 0.4}, {0.2}}, {{1.0, 2.0}, {1.0}}); + CHECK_THAT(data.energy[0][0], WithinAbs(1.0 / 3.0, 1e-12)); + CHECK_THAT(data.energy[0][1], WithinAbs(2.0 / 3.0, 1e-12)); + CHECK_THAT(data.energy[1][0], WithinAbs(1.0, 1e-12)); + check_cross_sections(data, {{0.2, 0.4}, {0.2}}, {{0.2, 0.2}, {0.2}}); +} + +TEMPLATE_TEST_CASE("Multigroup sampling with variable multiplicity", + "[scattdata]", openmc::ScattDataLegendre, openmc::ScattDataHistogram, + openmc::ScattDataTabular) +{ + auto data = make_scatter({{0.2, 0.4}, {0.2}}, {{1.0, 2.0}, {1.0}}); + check_sample_weights(data, 0, 1.5); + check_sample_weights(data, 1, 1.0); +} + +TEMPLATE_TEST_CASE("Multigroup sampling with constant multiplicity", + "[scattdata]", openmc::ScattDataLegendre, openmc::ScattDataHistogram, + openmc::ScattDataTabular) +{ + double multiplicity = 1.0; + SECTION("Unit multiplicity") {} + SECTION("Nonunit multiplicity") + { + multiplicity = 1.5; + } + + auto data = make_scatter( + {{0.2, 0.4}, {0.2}}, {{multiplicity, multiplicity}, {multiplicity}}); + CHECK_THAT(data.energy[0][0], WithinAbs(1.0 / 3.0, 1e-12)); + CHECK_THAT(data.energy[0][1], WithinAbs(2.0 / 3.0, 1e-12)); + check_cross_sections(data, {{0.2, 0.4}, {0.2}}, + {{0.2 / multiplicity, 0.4 / multiplicity}, {0.2 / multiplicity}}); + check_sample_weights(data, 0, multiplicity); + check_sample_weights(data, 1, multiplicity); +} + +TEMPLATE_TEST_CASE("Combining multigroup data with different multiplicities", + "[scattdata]", openmc::ScattDataLegendre, openmc::ScattDataHistogram, + openmc::ScattDataTabular) +{ + auto first = make_scatter({{0.2, 0.4}, {0.2}}, {{1.0, 2.0}, {1.0}}); + auto second = + make_scatter({{0.3, 0.1}, {0.3}}, {{3.0, 1.0}, {1.5}}); + TestType combined; + combined.combine({&first, &second}, {0.5, 2.0}); + + // Cross sections, rather than multiplicities, add with number density. + check_cross_sections(combined, {{0.7, 0.4}, {0.7}}, {{0.3, 0.3}, {0.5}}); + CHECK_THAT(combined.energy[0][0], WithinAbs(7.0 / 11.0, 1e-12)); + CHECK_THAT(combined.energy[0][1], WithinAbs(4.0 / 11.0, 1e-12)); + check_sample_weights(combined, 0, 11.0 / 6.0); + check_sample_weights(combined, 1, 1.4); +} + +TEST_CASE("Legendre to tabular conversion preserves variable multiplicity", + "[scattdata]") +{ + openmc::tensor::Tensor gmin({2}, 0); + openmc::tensor::Tensor gmax({2}, 1); + gmin[1] = 1; + openmc::ScattDataLegendre legendre; + // A positive first-order angular distribution exercises the conversion + // without requiring the negative-probability correction. + legendre.init(gmin, gmax, {{1.0, 2.0}, {1.0}}, + {{{0.2, 0.04}, {0.4, 0.08}}, {{0.2, 0.04}}}); + + openmc::ScattDataTabular tabular; + convert_legendre_to_tabular(legendre, tabular); + check_cross_sections(tabular, {{0.2, 0.4}, {0.2}}, {{0.2, 0.2}, {0.2}}); + CHECK_THAT(tabular.energy[0][0], WithinAbs(1.0 / 3.0, 1e-12)); + CHECK_THAT(tabular.energy[0][1], WithinAbs(2.0 / 3.0, 1e-12)); + for (double mu : {-1.0, -0.5, 0.0, 0.5, 1.0}) { + CHECK_THAT( + tabular.calc_f(0, 0, mu), WithinAbs(legendre.calc_f(0, 0, mu), 1e-12)); + CHECK_THAT( + tabular.calc_f(0, 1, mu), WithinAbs(legendre.calc_f(0, 1, mu), 1e-12)); + } + check_sample_weights(tabular, 0, 1.5); + check_sample_weights(tabular, 1, 1.0); +} + +TEMPLATE_TEST_CASE("Zero multigroup scattering row has finite cross sections", + "[scattdata]", openmc::ScattDataLegendre, openmc::ScattDataHistogram, + openmc::ScattDataTabular) +{ + auto data = make_scatter({{0.0, 0.0}, {0.2}}, {{0.0, 0.0}, {1.0}}); + check_cross_sections(data, {{0.0, 0.0}, {0.2}}, {{0.0, 0.0}, {0.2}}); + CHECK(std::isfinite(data.get_xs(MgxsType::SCATTER, 0, nullptr, nullptr))); + // This incoming group has no scattering events and must never be sampled. + check_sample_weights(data, 1, 1.0); +} diff --git a/tests/regression_tests/mg_scattering_multiplicity/__init__.py b/tests/regression_tests/mg_scattering_multiplicity/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/mg_scattering_multiplicity/test.py b/tests/regression_tests/mg_scattering_multiplicity/test.py new file mode 100644 index 00000000000..72421fc5067 --- /dev/null +++ b/tests/regression_tests/mg_scattering_multiplicity/test.py @@ -0,0 +1,99 @@ +"""Analytical fixed-source checks for multigroup scattering multiplicity.""" + +import numpy as np +import pytest + +import openmc +from openmc.examples import slab_mg +from tests.regression_tests import config + + +@pytest.mark.parametrize('scatter_format', ['legendre', 'histogram', 'tabular']) +@pytest.mark.parametrize('multiplicity', [ + pytest.param(np.array([[1.0, 2.0], [1.0, 1.0]]), id='variable'), + pytest.param(np.ones((2, 2)), id='unit'), + pytest.param(np.full((2, 2), 1.5), id='constant-1.5'), +]) +@pytest.mark.parametrize('survival_biasing', [False, True], + ids=['analog', 'survival']) +def test_scattering_multiplicity(run_in_tmpdir, scatter_format, multiplicity, + survival_biasing): + # Matrices use incoming groups as rows, ordered from high to low energy. + total = np.ones(2) + absorption = np.array([0.6, 0.8]) + scattering = np.array([[0.2, 0.2], [0.0, 0.2]]) + production = scattering * multiplicity + groups = openmc.mgxs.EnergyGroups([0.0, 1.0, 2.0]) + data = openmc.XSdata('mat_1', groups) + data.order = 0 + data.set_total(total) + data.set_absorption(absorption) + data.set_scatter_matrix(production[:, :, np.newaxis]) + data.set_multiplicity_matrix(multiplicity) + if scatter_format != 'legendre': + data = data.convert_scatter_format(scatter_format, 2) + library = openmc.MGXSLibrary(groups) + library.add_xsdata(data) + library.export_to_hdf5('mgxs.h5') + + # Reuse the multigroup example, replacing its leaking slab by a unit-volume + # reflected box so that the spatially integrated balance has no leakage. + model = slab_mg(mgxslib_name='mgxs.h5') + box = openmc.model.RectangularParallelepiped( + 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, boundary_type='reflective') + cell = next(iter(model.geometry.root_universe.cells.values())) + cell.region = -box + model.geometry = openmc.Geometry([cell]) + settings = model.settings + settings.run_mode = 'fixed source' + settings.inactive = 0 + settings.batches = 20 + settings.particles = 5000 + settings.seed = 12345 + settings.survival_biasing = survival_biasing + settings.tabular_legendre = {'enable': False} + settings.source = openmc.IndependentSource( + space=openmc.stats.Box([0.0, 0.0, 0.0], [1.0, 1.0, 1.0]), + energy=openmc.stats.Discrete([1.5], [1.0])) + settings.output = {'summary': False, 'tallies': False} + + tally = openmc.Tally(name='group balance') + tally.filters = [openmc.EnergyFilter(groups.group_edges)] + tally.scores = ['flux', 'absorption'] + tally.estimator = 'tracklength' + model.tallies = [tally] + model.export_to_model_xml() + if config['build_inputs']: + return + + run_args = {'openmc_exec': config['exe'], 'event_based': config['event']} + if config['mpi']: + run_args['mpi_args'] = [config['mpiexec'], '-n', config['mpi_np']] + openmc.run(**run_args) + + # Unit source in the high-energy group: (diag(total) - B.T) Phi = Q. + # In the variable-multiplicity case, Phi = (1.25, 0.625), and the + # integrated absorption rate is 0.6*1.25 + 0.8*0.625 = 1.25. + expected_flux = np.linalg.solve(np.diag(total) - production.T, [1.0, 0.0]) + if multiplicity[0, 1] == 2.0: + np.testing.assert_allclose(expected_flux, [1.25, 0.625]) + expected_absorption = absorption * expected_flux + with openmc.StatePoint('statepoint.20.h5', autolink=False) as statepoint: + result = statepoint.get_tally(name='group balance') + for score, expected in [('flux', expected_flux), + ('absorption', expected_absorption)]: + # EnergyFilter bins are in increasing energy; reverse them to the + # high-to-low convention used by the input cross sections. + mean = result.get_values(scores=[score]).ravel()[::-1] + std_dev = result.get_values(scores=[score], + value='std_dev').ravel()[::-1] + assert np.all(np.isfinite(mean)) + assert np.all(np.isfinite(std_dev)) + # Require enough precision that a large error cannot be hidden by + # a broad confidence interval. The original variable-multiplicity + # defect biases the low-energy flux by about 23%, versus at most + # 9% allowed here. The 0.5% floor accommodates near-zero variance. + assert np.all(std_dev < 0.015 * expected), (score, std_dev, expected) + tolerance = np.maximum(6.0 * std_dev, 0.005 * expected) + assert np.all(np.abs(mean - expected) <= tolerance), ( + score, mean, expected, std_dev)