From 781dbf9c12b70937e9d3dc10c2ba29228e996b5c Mon Sep 17 00:00:00 2001 From: Alex Nellis <79525520+realmisch@users.noreply.github.com> Date: Fri, 26 Sep 2025 12:27:05 -0400 Subject: [PATCH] Multi-group capability for kinetics parameter calculations with Iterated Fission Probability (#3425) Co-authored-by: GuySten Co-authored-by: Paul Romano --- docs/source/usersguide/kinetics.rst | 25 +++++++- openmc/model/model.py | 41 +++++++++++-- openmc/statepoint.py | 58 +++++++++++++++++++ src/tallies/tally.cpp | 3 +- src/tallies/tally_scoring.cpp | 9 +++ .../ifp/groupwise/__init__.py | 0 .../ifp/groupwise/inputs_true.dat | 43 ++++++++++++++ .../ifp/groupwise/results_true.dat | 21 +++++++ tests/regression_tests/ifp/groupwise/test.py | 40 +++++++++++++ tests/regression_tests/ifp/total/__init__.py | 0 .../ifp/{ => total}/inputs_true.dat | 0 .../ifp/{ => total}/results_true.dat | 0 .../regression_tests/ifp/{ => total}/test.py | 1 - tests/unit_tests/test_ifp.py | 39 +++++++++++++ 14 files changed, 272 insertions(+), 8 deletions(-) create mode 100644 tests/regression_tests/ifp/groupwise/__init__.py create mode 100644 tests/regression_tests/ifp/groupwise/inputs_true.dat create mode 100644 tests/regression_tests/ifp/groupwise/results_true.dat create mode 100644 tests/regression_tests/ifp/groupwise/test.py create mode 100644 tests/regression_tests/ifp/total/__init__.py rename tests/regression_tests/ifp/{ => total}/inputs_true.dat (100%) rename tests/regression_tests/ifp/{ => total}/results_true.dat (100%) rename tests/regression_tests/ifp/{ => total}/test.py (99%) diff --git a/docs/source/usersguide/kinetics.rst b/docs/source/usersguide/kinetics.rst index bdf26d341b..9024ff8227 100644 --- a/docs/source/usersguide/kinetics.rst +++ b/docs/source/usersguide/kinetics.rst @@ -67,6 +67,23 @@ are needed to compute kinetics parameters in OpenMC: Obtaining kinetics parameters ----------------------------- +The ``Model`` class can be used to automatically generate all IFP tallies using +the Python API with :attr:`openmc.Settings.ifp_n_generation` greater than 0 and +the :meth:`openmc.Model.add_ifp_kinetics_tallies` method:: + + model = openmc.Model(geometry, settings=settings) + model.add_kinetics_parameters_tallies(num_groups=6) # Add 6 precursor groups + +Alternatively, each of the tallies can be manually defined using group-wise or +total :math:`\beta_{\text{eff}}` specified by providing a 6-group +:class:`openmc.DelayedGroupFilter`:: + + beta_tally = openmc.Tally(name="group-beta-score") + beta_tally.scores = ["ifp-beta-numerator"] + + # Add DelayedGroupFilter to enable group-wise tallies + beta_tally.filters = [openmc.DelayedGroupFilter(list(range(1, 7)))] + Here is an example showing how to declare the three available IFP scores in a single tally:: @@ -95,6 +112,12 @@ for ``ifp-denominator``: \beta_{\text{eff}} = \frac{S_{\text{ifp-beta-numerator}}}{S_{\text{ifp-denominator}}} +The kinetics parameters can be retrieved directly from a statepoint file using +the :meth:`openmc.StatePoint.ifp_results` method:: + + with openmc.StatePoint(output_path) as sp: + generation_time, beta_eff = sp.get_kinetics_parameters() + .. only:: html .. rubric:: References @@ -107,4 +130,4 @@ for ``ifp-denominator``: of the Iterated Fission Probability Method in OpenMC to Compute Adjoint-Weighted Kinetics Parameters", International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering (M&C 2025), Denver, April 27-30, - 2025 (to be presented). + 2025. diff --git a/openmc/model/model.py b/openmc/model/model.py index 59e6fa511e..c1ffafafd3 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -1,7 +1,7 @@ from __future__ import annotations from collections.abc import Iterable, Sequence import copy -from functools import lru_cache +from functools import cache from pathlib import Path import math from numbers import Integral, Real @@ -160,7 +160,7 @@ class Model: return False @property - @lru_cache(maxsize=None) + @cache def _materials_by_id(self) -> dict: """Dictionary mapping material ID --> material""" if self.materials: @@ -170,14 +170,14 @@ class Model: return {mat.id: mat for mat in mats} @property - @lru_cache(maxsize=None) + @cache def _cells_by_id(self) -> dict: """Dictionary mapping cell ID --> cell""" cells = self.geometry.get_all_cells() return {cell.id: cell for cell in cells.values()} @property - @lru_cache(maxsize=None) + @cache def _cells_by_name(self) -> dict[int, openmc.Cell]: # Get the names maps, but since names are not unique, store a set for # each name key. In this way when the user requests a change by a name, @@ -190,7 +190,7 @@ class Model: return result @property - @lru_cache(maxsize=None) + @cache def _materials_by_name(self) -> dict[int, openmc.Material]: if self.materials is None: mats = self.geometry.get_all_materials().values() @@ -203,6 +203,37 @@ class Model: result[mat.name].add(mat) return result + def add_kinetics_parameters_tallies(self, num_groups: int | None = None): + """Add tallies for calculating kinetics parameters using the IFP method. + + This method adds tallies to the model for calculating two kinetics + parameters, the generation time and the effective delayed neutron + fraction (beta effective). After a model is run, these parameters can be + determined through the :meth:`openmc.StatePoint.ifp_results` method. + + Parameters + ---------- + num_groups : int, optional + Number of precursor groups to filter the delayed neutron fraction. + If None, only the total effective delayed neutron fraction is + tallied. + + """ + if not any('ifp-time-numerator' in t.scores for t in self.tallies): + gen_time_tally = openmc.Tally(name='IFP time numerator') + gen_time_tally.scores = ['ifp-time-numerator'] + self.tallies.append(gen_time_tally) + if not any('ifp-beta-numerator' in t.scores for t in self.tallies): + beta_tally = openmc.Tally(name='IFP beta numerator') + beta_tally.scores = ['ifp-beta-numerator'] + if num_groups is not None: + beta_tally.filters = [openmc.DelayedGroupFilter(list(range(1, num_groups + 1)))] + self.tallies.append(beta_tally) + if not any('ifp-denominator' in t.scores for t in self.tallies): + denom_tally = openmc.Tally(name='IFP denominator') + denom_tally.scores = ['ifp-denominator'] + self.tallies.append(denom_tally) + @classmethod def from_xml( cls, diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 715becf488..29c11921cb 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -1,4 +1,5 @@ from datetime import datetime +from collections import namedtuple import glob import re import os @@ -8,6 +9,7 @@ import h5py import numpy as np from pathlib import Path from uncertainties import ufloat +from uncertainties.unumpy import uarray import openmc import openmc.checkvalue as cv @@ -15,6 +17,9 @@ import openmc.checkvalue as cv _VERSION_STATEPOINT = 18 +KineticsParameters = namedtuple("KineticsParameters", ["generation_time", "beta_effective"]) + + class StatePoint: """State information on a simulation at a certain point in time (at the end of a given batch). Statepoints can be used to analyze tally results as well @@ -710,3 +715,56 @@ class StatePoint: tally_filter.paths = cell.paths self._summary = summary + + def get_kinetics_parameters(self) -> KineticsParameters: + """Get kinetics parameters from IFP tallies. + + This method searches the tallies in the statepoint for the tallies + required to compute kinetics parameters using the Iterated Fission + Probability (IFP) method. + + Returns + ------- + KineticsParameters + A named tuple containing the generation time and effective delayed + neutron fraction. If the necessary tallies for one or both + parameters are not found, that parameter is returned as None. + + """ + + denom_tally = None + gen_time_tally = None + beta_tally = None + for tally in self.tallies.values(): + if 'ifp-denominator' in tally.scores: + denom_tally = self.get_tally(scores=['ifp-denominator']) + if 'ifp-time-numerator' in tally.scores: + gen_time_tally = self.get_tally(scores=['ifp-time-numerator']) + if 'ifp-beta-numerator' in tally.scores: + beta_tally = self.get_tally(scores=['ifp-beta-numerator']) + + if denom_tally is None: + return KineticsParameters(None, None) + + def get_ufloat(tally, score): + return uarray(tally.get_values(scores=[score]), + tally.get_values(scores=[score], value='std_dev')) + + denom_values = get_ufloat(denom_tally, 'ifp-denominator') + if gen_time_tally is None: + generation_time = None + else: + gen_time_values = get_ufloat(gen_time_tally, 'ifp-time-numerator') + gen_time_values /= denom_values*self.keff + generation_time = gen_time_values.flatten()[0] + + if beta_tally is None: + beta_effective = None + else: + beta_values = get_ufloat(beta_tally, 'ifp-beta-numerator') + beta_values /= denom_values + beta_effective = beta_values.flatten() + if beta_effective.size == 1: + beta_effective = beta_effective[0] + + return KineticsParameters(generation_time, beta_effective) diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index ae0bffe6eb..b9c615ecb5 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -560,7 +560,8 @@ void Tally::set_scores(const vector& scores) // Make sure a delayed group filter wasn't used with an incompatible // score. if (delayedgroup_filter_ != C_NONE) { - if (score_str != "delayed-nu-fission" && score_str != "decay-rate") + if (score_str != "delayed-nu-fission" && score_str != "decay-rate" && + score_str != "ifp-beta-numerator") fatal_error("Cannot tally " + score_str + "with a delayedgroup filter"); } diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 0df80a2398..67e851644a 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -964,6 +964,15 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, if (delayed_groups.size() == settings::ifp_n_generation) { if (delayed_groups[0] > 0) { score = p.wgt_last(); + if (tally.delayedgroup_filter_ != C_NONE) { + auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_]; + const DelayedGroupFilter& filt { + *dynamic_cast( + model::tally_filters[i_dg_filt].get())}; + score_fission_delayed_dg(i_tally, delayed_groups[0] - 1, + score, score_index, p.filter_matches()); + continue; + } } } } diff --git a/tests/regression_tests/ifp/groupwise/__init__.py b/tests/regression_tests/ifp/groupwise/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/regression_tests/ifp/groupwise/inputs_true.dat b/tests/regression_tests/ifp/groupwise/inputs_true.dat new file mode 100644 index 0000000000..6d7e20717b --- /dev/null +++ b/tests/regression_tests/ifp/groupwise/inputs_true.dat @@ -0,0 +1,43 @@ + + + + + + + + + + + + + + eigenvalue + 1000 + 20 + 5 + + + -10.0 -10.0 -10.0 10.0 10.0 10.0 + + + true + + + 5 + + + + 1 2 3 4 5 6 + + + ifp-time-numerator + + + 1 + ifp-beta-numerator + + + ifp-denominator + + + diff --git a/tests/regression_tests/ifp/groupwise/results_true.dat b/tests/regression_tests/ifp/groupwise/results_true.dat new file mode 100644 index 0000000000..ea66a8de3c --- /dev/null +++ b/tests/regression_tests/ifp/groupwise/results_true.dat @@ -0,0 +1,21 @@ +k-combined: +1.006559E+00 5.389391E-03 +tally 1: +9.109384E-08 +5.667165E-16 +tally 2: +3.000000E-03 +9.000000E-06 +0.000000E+00 +0.000000E+00 +2.100000E-02 +1.370000E-04 +2.800000E-02 +2.220000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +tally 3: +1.489000E+01 +1.480036E+01 diff --git a/tests/regression_tests/ifp/groupwise/test.py b/tests/regression_tests/ifp/groupwise/test.py new file mode 100644 index 0000000000..a1a0ebefb8 --- /dev/null +++ b/tests/regression_tests/ifp/groupwise/test.py @@ -0,0 +1,40 @@ +"""Test the Iterated Fission Probability (IFP) method to compute adjoint-weighted +kinetics parameters using dedicated tallies.""" + +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + +@pytest.fixture() +def ifp_model(): + # Material + material = openmc.Material(name="core") + material.add_nuclide("U235", 1.0) + material.set_density('g/cm3', 16.0) + + # Geometry + radius = 10.0 + sphere = openmc.Sphere(r=radius, boundary_type="vacuum") + cell = openmc.Cell(region=-sphere, fill=material) + geometry = openmc.Geometry([cell]) + + # Settings + settings = openmc.Settings() + settings.particles = 1000 + settings.batches = 20 + settings.inactive = 5 + settings.ifp_n_generation = 5 + + model = openmc.Model(settings=settings, geometry=geometry) + + space = openmc.stats.Box(*cell.bounding_box) + model.settings.source = openmc.IndependentSource( + space=space, constraints={'fissionable': True}) + model.add_kinetics_parameters_tallies(num_groups=6) + return model + + +def test_iterated_fission_probability(ifp_model): + harness = PyAPITestHarness("statepoint.20.h5", model=ifp_model) + harness.main() diff --git a/tests/regression_tests/ifp/total/__init__.py b/tests/regression_tests/ifp/total/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/regression_tests/ifp/inputs_true.dat b/tests/regression_tests/ifp/total/inputs_true.dat similarity index 100% rename from tests/regression_tests/ifp/inputs_true.dat rename to tests/regression_tests/ifp/total/inputs_true.dat diff --git a/tests/regression_tests/ifp/results_true.dat b/tests/regression_tests/ifp/total/results_true.dat similarity index 100% rename from tests/regression_tests/ifp/results_true.dat rename to tests/regression_tests/ifp/total/results_true.dat diff --git a/tests/regression_tests/ifp/test.py b/tests/regression_tests/ifp/total/test.py similarity index 99% rename from tests/regression_tests/ifp/test.py rename to tests/regression_tests/ifp/total/test.py index 6969a54c49..18b89cfc0b 100644 --- a/tests/regression_tests/ifp/test.py +++ b/tests/regression_tests/ifp/total/test.py @@ -6,7 +6,6 @@ import pytest from tests.testing_harness import PyAPITestHarness - @pytest.fixture() def ifp_model(): model = openmc.Model() diff --git a/tests/unit_tests/test_ifp.py b/tests/unit_tests/test_ifp.py index 8d0fd98010..e527f16246 100644 --- a/tests/unit_tests/test_ifp.py +++ b/tests/unit_tests/test_ifp.py @@ -47,3 +47,42 @@ def test_exceptions(options, error, run_in_tmpdir, geometry): tallies = openmc.Tallies([tally]) model = openmc.Model(geometry=geometry, settings=settings, tallies=tallies) model.run() + + +@pytest.mark.parametrize( + "num_groups, use_auto_tallies", + [ + (None, True), + (None, False), + (6, True), + (6, False), + ], +) +def test_get_kinetics_parameters(run_in_tmpdir, geometry, num_groups, use_auto_tallies): + # Create basic model + model = openmc.Model(geometry=geometry) + model.settings.particles = 1000 + model.settings.batches = 20 + model.settings.inactive = 5 + model.settings.ifp_n_generation = 5 + + # Add IFP tallies either via the convenience method or manually + if use_auto_tallies: + model.add_kinetics_parameters_tallies(num_groups=num_groups) + else: + for score in ["ifp-time-numerator", "ifp-beta-numerator", "ifp-denominator"]: + tally = openmc.Tally() + tally.scores = [score] + if score == "ifp-beta-numerator" and num_groups is not None: + tally.filters = [openmc.DelayedGroupFilter(list(range(1, num_groups + 1)))] + model.tallies.append(tally) + + # Run and get kinetics parameters + sp_file = model.run() + with openmc.StatePoint(sp_file) as sp: + params = sp.get_kinetics_parameters() + assert isinstance(params, openmc.KineticsParameters) + assert params.generation_time is not None + assert params.beta_effective is not None + if num_groups is not None: + assert len(params.beta_effective) == num_groups