Multi-group capability for kinetics parameter calculations with Iterated Fission Probability (#3425)

Co-authored-by: GuySten <guyste@post.bgu.ac.il>
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

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.

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,

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)

src/tallies/tally.cpp

@@ -560,7 +560,8 @@ void Tally::set_scores(const vector<std::string>& 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");
     }
 

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<DelayedGroupFilter*>(
+                      model::tally_filters[i_dg_filt].get())};
+                  score_fission_delayed_dg(i_tally, delayed_groups[0] - 1,
+                    score, score_index, p.filter_matches());
+                  continue;
+                }
               }
             }
           }

tests/regression_tests/ifp/groupwise/inputs_true.dat

@@ -0,0 +1,43 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1" name="core" depletable="true">
+      <density value="16.0" units="g/cm3"/>
+      <nuclide name="U235" ao="1.0"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1" universe="1"/>
+    <surface id="1" type="sphere" boundary="vacuum" coeffs="0.0 0.0 0.0 10.0"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>1000</particles>
+    <batches>20</batches>
+    <inactive>5</inactive>
+    <source type="independent" strength="1.0" particle="neutron">
+      <space type="box">
+        <parameters>-10.0 -10.0 -10.0 10.0 10.0 10.0</parameters>
+      </space>
+      <constraints>
+        <fissionable>true</fissionable>
+      </constraints>
+    </source>
+    <ifp_n_generation>5</ifp_n_generation>
+  </settings>
+  <tallies>
+    <filter id="1" type="delayedgroup">
+      <bins>1 2 3 4 5 6</bins>
+    </filter>
+    <tally id="1" name="IFP time numerator">
+      <scores>ifp-time-numerator</scores>
+    </tally>
+    <tally id="2" name="IFP beta numerator">
+      <filters>1</filters>
+      <scores>ifp-beta-numerator</scores>
+    </tally>
+    <tally id="3" name="IFP denominator">
+      <scores>ifp-denominator</scores>
+    </tally>
+  </tallies>
+</model>

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

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()

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()

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