From 72be5291cf6014f065c9ae3043a2084e1764832a Mon Sep 17 00:00:00 2001 From: Eden Rochman Date: Thu, 6 Aug 2026 16:27:18 +0300 Subject: [PATCH 1/6] Add k-infinity estimator for IndependentOperator Add calculate_kinf parameter to IndependentOperator that estimates the infinite multiplication factor from one-group reaction rates, following OpenMC's keff convention for (n,xn) reactions. Add nu-fission one-group cross section support to MicroXS. --- docs/source/usersguide/depletion.rst | 16 +++ openmc/deplete/independent_operator.py | 114 ++++++++++++++++- openmc/deplete/microxs.py | 115 ++++++++++++++++-- .../test_deplete_independent_operator.py | 97 +++++++++++++++ tests/unit_tests/test_deplete_microxs.py | 27 ++++ 5 files changed, 356 insertions(+), 13 deletions(-) diff --git a/docs/source/usersguide/depletion.rst b/docs/source/usersguide/depletion.rst index 261900ce61a..0459ade02d6 100644 --- a/docs/source/usersguide/depletion.rst +++ b/docs/source/usersguide/depletion.rst @@ -270,6 +270,22 @@ transport-depletion calculation and follow the same steps from there. the depletion chain with at least one reaction, that reaction will not be simulated. +If the microscopic cross section data includes 'fission' and 'nu-fission' +cross sections, :class:`~openmc.deplete.IndependentOperator` can also estimate +the infinite multiplication factor at each depletion step by passing +``calculate_kinf=True``:: + + op = openmc.deplete.IndependentOperator(materials, fluxes, micros, + calculate_kinf=True) + +The estimate is computed as the ratio of the neutron production rate to the +neutron loss rate based on the one-group reaction rates and is reported as the +eigenvalue in the depletion results, which can be retrieved with +:meth:`~openmc.deplete.Results.get_keff`. Consistent with the definition of +the multiplication factor used elsewhere in OpenMC, neutrons produced in +(n,xn) reactions are not counted as production; instead, each (n,xn) reaction +reduces the loss term by :math:`x - 1`. + .. _micros: Loading and Generating Microscopic Cross Sections diff --git a/openmc/deplete/independent_operator.py b/openmc/deplete/independent_operator.py index c12863956b9..42eb92c0b57 100644 --- a/openmc/deplete/independent_operator.py +++ b/openmc/deplete/independent_operator.py @@ -8,6 +8,7 @@ from __future__ import annotations from collections.abc import Iterable import copy +import re import numpy as np from uncertainties import ufloat @@ -22,6 +23,18 @@ from .results import Results from .helpers import ChainFissionHelper, ConstantFissionYieldHelper, SourceRateHelper +# Regular expression matching reactions that emit one or more neutrons, e.g., +# (n,2n) or (n,np), with the number of emitted neutrons captured +_XN_REACTION = re.compile(r'\(n,(\d*)n') + + +def _neutrons_emitted(reaction: str) -> int: + """Number of neutrons in the exit channel of a transmutation reaction.""" + match = _XN_REACTION.match(reaction) + if match is None: + return 0 + return int(match.group(1)) if match.group(1) else 1 + class IndependentOperator(OpenMCOperator): """Transport-independent transport operator based on multigroup data. @@ -55,6 +68,14 @@ class IndependentOperator(OpenMCOperator): Defaults to ``openmc.config['chain_file']``. keff : 2-tuple of float, optional keff eigenvalue and uncertainty from transport calculation. + calculate_kinf : bool, optional + If True, the infinite multiplication factor is estimated from the + material compositions and one-group cross sections at each depletion + step and reported as the eigenvalue in the depletion results. Requires + that each :class:`~openmc.deplete.MicroXS` instance contains 'fission' + and 'nu-fission' cross sections. Mutually exclusive with ``keff``. + + .. versionadded:: 0.15.4 prev_results : Results, optional Results from a previous depletion calculation. normalization_mode : {"fission-q", "source-rate"} @@ -116,7 +137,8 @@ def __init__(self, fission_q=None, prev_results=None, reduce_chain_level=None, - fission_yield_opts=None): + fission_yield_opts=None, + calculate_kinf=False): # Validate micro-xs parameters check_type('materials', materials, Iterable, openmc.Material) check_type('micros', micros, Iterable, MicroXS) @@ -134,6 +156,20 @@ def __init__(self, self._keff = keff + check_type('calculate_kinf', calculate_kinf, bool) + if calculate_kinf: + if keff is not None: + raise ValueError("The 'keff' and 'calculate_kinf' arguments " + "are mutually exclusive.") + for micro in micros: + missing = {'fission', 'nu-fission'} - set(micro.reactions) + if missing: + raise ValueError( + "Estimating k-infinity requires 'fission' and " + "'nu-fission' cross sections in each MicroXS " + f"instance (missing {sorted(missing)}).") + self._calculate_kinf = calculate_kinf + if fission_yield_opts is None: fission_yield_opts = {} helper_kwargs = {'normalization_mode': normalization_mode, @@ -165,7 +201,8 @@ def from_nuclides(cls, volume, nuclides, fission_q=None, prev_results=None, reduce_chain_level=None, - fission_yield_opts=None): + fission_yield_opts=None, + calculate_kinf=False): """ Alternate constructor from a dictionary of nuclide concentrations @@ -187,6 +224,13 @@ def from_nuclides(cls, volume, nuclides, keff : 2-tuple of float, optional keff eigenvalue and uncertainty from transport calculation. Default is None. + calculate_kinf : bool, optional + If True, the infinite multiplication factor is estimated from the + material compositions and one-group cross sections at each + depletion step. Requires that ``micro_xs`` contains 'fission' and + 'nu-fission' cross sections. Mutually exclusive with ``keff``. + + .. versionadded:: 0.15.4 normalization_mode : {"fission-q", "source-rate"} Indicate how reaction rates should be calculated. ``"fission-q"`` uses the fission Q values from the depletion @@ -222,7 +266,8 @@ def from_nuclides(cls, volume, nuclides, fission_q=fission_q, prev_results=prev_results, reduce_chain_level=reduce_chain_level, - fission_yield_opts=fission_yield_opts) + fission_yield_opts=fission_yield_opts, + calculate_kinf=calculate_kinf) @staticmethod def _consolidate_nuclides_to_material(nuclides, nuc_units, volume): @@ -407,14 +452,73 @@ def __call__(self, vec, source_rate) -> OperatorResult: if source_rate == 0.0: rates = self.reaction_rates.copy() rates.fill(0.0) - return OperatorResult(ufloat(0.0, 0.0), rates) + if self._calculate_kinf: + keff = self._estimate_k_inf() + else: + keff = ufloat(0.0, 0.0) + return OperatorResult(keff, rates) rates = self._calculate_reaction_rates(source_rate) - keff = self._keff + if self._calculate_kinf: + keff = self._estimate_k_inf() + else: + keff = self._keff op_result = OperatorResult(keff, rates) return copy.deepcopy(op_result) + def _estimate_k_inf(self): + r"""Estimate the infinite multiplication factor. + + The estimate is computed as the ratio of the neutron production rate + to the neutron loss rate: + + .. math:: + k_\infty = \frac{\sum_i N_i (\nu\sigma_f)_i} + {\sum_i N_i \sum_j (1 - x_j) \sigma_{i,j}} + + where :math:`N_i` is the number of atoms of nuclide :math:`i`, + :math:`(\nu\sigma_f)_i` is its one-group fission neutron production + cross section, :math:`\sigma_{i,j}` is the one-group cross section of + transmutation reaction :math:`j`, and :math:`x_j` is the number of + neutrons emitted by reaction :math:`j`. This is consistent with the + definition of the multiplication factor used elsewhere in OpenMC: + neutrons produced in (n,xn) reactions are not counted as production; + instead, each (n,xn) reaction reduces the loss term by :math:`x - 1`. + + Returns + ------- + uncertainties.UFloat + Estimated k-infinity with zero uncertainty + + """ + production = 0.0 + loss = 0.0 + for mat in self.local_mats: + i_mat = self._mat_index_map[mat] + flux = self.fluxes[i_mat] + micro_xs = self.cross_sections[i_mat] + for nuc in micro_xs.nuclides: + if nuc not in self.number.index_nuc: + continue + atoms = self.number[mat, nuc] + if atoms <= 0.0: + continue + for rxn in micro_xs.reactions: + rate = atoms * (micro_xs[nuc, rxn] * flux).sum() + if rxn == 'nu-fission': + production += rate + elif rxn != 'damage-energy': + loss += (1 - _neutrons_emitted(rxn)) * rate + + # Sum contributions over all MPI processes + production = comm.allreduce(production) + loss = comm.allreduce(loss) + + if loss <= 0.0: + return ufloat(0.0, 0.0) + return ufloat(production / loss, 0.0) + def _update_materials(self): """Updates material compositions in OpenMC on all processes.""" diff --git a/openmc/deplete/microxs.py b/openmc/deplete/microxs.py index 687cf646f29..b0279f8defa 100644 --- a/openmc/deplete/microxs.py +++ b/openmc/deplete/microxs.py @@ -17,7 +17,7 @@ from openmc.checkvalue import check_type, check_value, check_iterable_type, PathLike from openmc import StatePoint from openmc.mgxs import GROUP_STRUCTURES -from openmc.data import REACTION_MT +from openmc.data import DataLibrary, REACTION_MT, Reaction import openmc from .chain import Chain, REACTIONS, _get_chain from .coupled_operator import _find_cross_sections, _get_nuclides_with_data @@ -28,6 +28,7 @@ _valid_rxns = list(REACTIONS) _valid_rxns.append('fission') _valid_rxns.append('damage-energy') +_valid_rxns.append('nu-fission') # TODO: Replace with type statement when support is Python 3.12+ @@ -81,7 +82,10 @@ def get_microxs_and_flux( nuclides from the depletion chain file are used. reactions : list of str Reactions to get cross sections for. If not specified, all neutron - reactions listed in the depletion chain file are used. + reactions listed in the depletion chain file are used. In addition to + transmutation reactions, 'nu-fission' may be specified to obtain the + fission neutron production cross section, which is needed to estimate + k-infinity with :class:`~openmc.deplete.IndependentOperator`. energies : iterable of float or str Energy group boundaries in [eV] or the name of the group structure. If left as None energies will default to [0.0, 100e6] @@ -303,6 +307,82 @@ def get_microxs_and_flux( return fluxes, micros +def _collapse_nu_fission( + path: PathLike, + nuclide: str, + temperature: float, + energies: Sequence[float], + flux: np.ndarray +) -> float: + r"""Compute a one-group fission neutron production cross section. + + The fission neutron production cross section, + :math:`\nu(E)\sigma_f(E)`, is integrated against a flux that is assumed + to be constant in energy within each group, matching the treatment used + for other reactions in :meth:`openmc.lib.Nuclide.collapse_rate`. + + Parameters + ---------- + path : PathLike + Path to the HDF5 data file containing the nuclide. + nuclide : str + Name of the nuclide, e.g., 'U235'. + temperature : float + Temperature in [K]. The closest available temperature is used for the + fission cross section. + energies : iterable of float + Energy group boundaries in [eV] in ascending order. + flux : numpy.ndarray + Flux in each energy group, normalized to sum to unity. + + Returns + ------- + float + Flux-averaged fission neutron production cross section in [b]. Zero if + the nuclide has no fission data. + + """ + with h5py.File(path, 'r') as h5: + group = h5[nuclide] + if 'reactions/reaction_018' not in group: + return 0.0 + + # Select the available temperature closest to the requested one + temp_keys = list(group['energy']) + temps = np.array([float(t[:-1]) for t in temp_keys]) + temp_key = temp_keys[np.argmin(np.abs(temps - temperature))] + + energy_grid = {temp_key: group['energy'][temp_key][()]} + rx = Reaction.from_hdf5(group['reactions/reaction_018'], energy_grid) + + xs = rx.xs[temp_key] + + # Total nu(E) is the sum of the yields of all neutron products. If a + # product with emission mode 'total' is present, use it alone to avoid + # double counting prompt and delayed neutrons. + neutron_products = [p for p in rx.products if p.particle == 'neutron'] + total_products = [p for p in neutron_products if p.emission_mode == 'total'] + if total_products: + neutron_products = total_products + if not neutron_products: + return 0.0 + + def nu(e): + return sum(p.yield_(e) for p in neutron_products) + + # Integrate nu(E)*sigma_f(E) against a histogram flux + nu_fission = 0.0 + for g, flux_g in enumerate(flux): + if flux_g == 0.0: + continue + e_low, e_high = energies[g], energies[g + 1] + inside = xs.x[(xs.x > e_low) & (xs.x < e_high)] + e = np.concatenate([[e_low], inside, [e_high]]) + nu_fission += np.trapezoid(nu(e) * xs(e), e) * flux_g / (e_high - e_low) + + return nu_fission + + class MicroXS: """Microscopic cross section data for use in transport-independent depletion. @@ -385,7 +465,14 @@ def from_multigroup_flux( nuclides from the depletion chain file are used. reactions : list of str, optional Reactions to get cross sections for. If not specified, all neutron - reactions listed in the depletion chain file are used. + reactions listed in the depletion chain file are used. In addition + to transmutation reactions, 'nu-fission' may be specified to + obtain the fission neutron production cross section, which is + needed to estimate k-infinity with + :class:`~openmc.deplete.IndependentOperator`. + + .. versionchanged:: 0.15.4 + Added support for 'nu-fission'. **init_kwargs : dict Keyword arguments passed to :func:`openmc.lib.init` @@ -418,10 +505,12 @@ def from_multigroup_flux( nuclides = [nuc.name for nuc in nuclides] # Get reaction MT values. If no reactions specified, default to the - # reactions available in the chain file + # reactions available in the chain file. The 'nu-fission' reaction is + # handled separately since it does not correspond to a single MT value. if reactions is None: reactions = chain.reactions - mts = [REACTION_MT[name] for name in reactions] + mts = [REACTION_MT[name] if name != 'nu-fission' else None + for name in reactions] # Create 3D array for microscopic cross sections microxs_arr = np.zeros((len(nuclides), len(mts), 1)) @@ -434,6 +523,10 @@ def from_multigroup_flux( # Normalize multigroup flux multigroup_flux /= flux_sum + # If nu-fission was requested, get paths to pointwise data files + if 'nu-fission' in reactions: + data_library = DataLibrary.from_xml(cross_sections) + # Compute microscopic cross sections within a temporary session with openmc.lib.TemporarySession(**init_kwargs): # For each nuclide and reaction, compute the flux-averaged xs @@ -442,9 +535,15 @@ def from_multigroup_flux( continue lib_nuc = openmc.lib.load_nuclide(nuc) for mt_index, mt in enumerate(mts): - microxs_arr[nuc_index, mt_index, 0] = lib_nuc.collapse_rate( - mt, temperature, energies, multigroup_flux - ) + if mt is None: + path = data_library.get_by_material(nuc)['path'] + microxs_arr[nuc_index, mt_index, 0] = \ + _collapse_nu_fission(path, nuc, temperature, + energies, multigroup_flux) + else: + microxs_arr[nuc_index, mt_index, 0] = \ + lib_nuc.collapse_rate( + mt, temperature, energies, multigroup_flux) return cls(microxs_arr, nuclides, reactions) diff --git a/tests/unit_tests/test_deplete_independent_operator.py b/tests/unit_tests/test_deplete_independent_operator.py index aca83399a08..3f682b1d686 100644 --- a/tests/unit_tests/test_deplete_independent_operator.py +++ b/tests/unit_tests/test_deplete_independent_operator.py @@ -4,10 +4,12 @@ from pathlib import Path +import numpy as np import pytest from openmc import Material from openmc.deplete import IndependentOperator, MicroXS, Chain +from openmc.deplete.independent_operator import _neutrons_emitted CHAIN_PATH = Path(__file__).parents[1] / "chain_simple.xml" ONE_GROUP_XS = Path(__file__).parents[1] / "micro_xs_simple.csv" @@ -53,3 +55,98 @@ def test_error_handling(): micros = [micro_xs] with pytest.raises(ValueError, match=r"The length of fluxes \(2\)"): IndependentOperator(materials, fluxes, micros, CHAIN_PATH) + + +def _uranium_material(): + fuel = Material(name="u metal") + fuel.add_nuclide("U235", 0.05) + fuel.add_nuclide("U238", 0.95) + fuel.set_density("g/cc", 19.0) + fuel.depletable = True + fuel.volume = 1.0 + return fuel + + +def test_neutrons_emitted(): + assert _neutrons_emitted('fission') == 0 + assert _neutrons_emitted('(n,gamma)') == 0 + assert _neutrons_emitted('(n,p)') == 0 + assert _neutrons_emitted('(n,a)') == 0 + assert _neutrons_emitted('(n,3He)') == 0 + assert _neutrons_emitted('(n,2a)') == 0 + assert _neutrons_emitted('(n,np)') == 1 + assert _neutrons_emitted('(n,nd2a)') == 1 + assert _neutrons_emitted('(n,2n)') == 2 + assert _neutrons_emitted('(n,2nd)') == 2 + assert _neutrons_emitted('(n,3n)') == 3 + assert _neutrons_emitted('(n,3np)') == 3 + assert _neutrons_emitted('(n,4n)') == 4 + + +def test_calculate_kinf(): + # Only U235 has nonzero cross sections so that the k-infinity estimate + # does not depend on the material composition + nuclides = ['U235', 'U238'] + reactions = ['fission', 'nu-fission', '(n,gamma)', '(n,2n)'] + data = np.array([ + [[50.0], [120.0], [10.0], [2.0]], + [[0.0], [0.0], [0.0], [0.0]], + ]) + micro_xs = MicroXS(data, nuclides, reactions) + + op = IndependentOperator( + [_uranium_material()], [np.array([1.0])], [micro_xs], CHAIN_PATH, + normalization_mode='source-rate', calculate_kinf=True) + vec = op.initial_condition() + + # Production is nu-fission; loss is fission + (n,gamma) - (n,2n), since + # (n,2n) produces one net neutron and is not counted as production + expected = 120.0 / (50.0 + 10.0 - 2.0) + result = op(vec, 1.0) + assert result.k.n == pytest.approx(expected) + assert result.k.s == 0.0 + + # A decay step (zero source rate) reports the same estimate + result = op(vec, 0.0) + assert result.k.n == pytest.approx(expected) + + +def test_calculate_kinf_multigroup(): + # Two-group cross sections and flux, U235 only so that the estimate does + # not depend on the material composition + nuclides = ['U235'] + reactions = ['fission', 'nu-fission', '(n,gamma)'] + data = np.array([[[10.0, 50.0], [25.0, 120.0], [2.0, 10.0]]]) + micro_xs = MicroXS(data, nuclides, reactions) + flux = np.array([0.75, 0.25]) + + fuel = Material(name="u metal") + fuel.add_nuclide("U235", 1.0) + fuel.set_density("g/cc", 19.0) + fuel.depletable = True + fuel.volume = 1.0 + + op = IndependentOperator( + [fuel], [flux], [micro_xs], CHAIN_PATH, + normalization_mode='source-rate', calculate_kinf=True) + vec = op.initial_condition() + + production = 25.0*0.75 + 120.0*0.25 + loss = (10.0 + 2.0)*0.75 + (50.0 + 10.0)*0.25 + result = op(vec, 1.0) + assert result.k.n == pytest.approx(production / loss) + + +def test_calculate_kinf_errors(): + # MicroXS without nu-fission data cannot be used to estimate k-infinity + micro_xs = MicroXS.from_csv(ONE_GROUP_XS) + with pytest.raises(ValueError, match="nu-fission"): + IndependentOperator([_uranium_material()], [1.0], [micro_xs], + CHAIN_PATH, calculate_kinf=True) + + # keff and calculate_kinf are mutually exclusive + data = np.zeros((1, 2, 1)) + micro_xs = MicroXS(data, ['U235'], ['fission', 'nu-fission']) + with pytest.raises(ValueError, match="mutually exclusive"): + IndependentOperator([_uranium_material()], [1.0], [micro_xs], + CHAIN_PATH, keff=(1.0, 0.0), calculate_kinf=True) diff --git a/tests/unit_tests/test_deplete_microxs.py b/tests/unit_tests/test_deplete_microxs.py index 26529e6ce96..674f1488857 100644 --- a/tests/unit_tests/test_deplete_microxs.py +++ b/tests/unit_tests/test_deplete_microxs.py @@ -116,6 +116,33 @@ def test_multigroup_flux_same(): assert microxs_4g.data == pytest.approx(microxs_2g.data) +def test_from_multigroup_flux_nu_fission(): + chain_file = Path(__file__).parents[1] / 'chain_simple.xml' + energies = [0., 6.25e-1, 5.53e3, 8.21e5, 2.e7] + + # For a thermal flux, the average number of neutrons per U235 fission + # (including delayed neutrons) should be about 2.43 + flux = [1.0, 0., 0., 0.] + microxs = MicroXS.from_multigroup_flux( + energies=energies, multigroup_flux=flux, chain_file=chain_file, + nuclides=['U235', 'O16'], reactions=['fission', 'nu-fission']) + assert microxs.reactions == ['fission', 'nu-fission'] + nu_bar = microxs['U235', 'nu-fission'][0] / microxs['U235', 'fission'][0] + assert nu_bar == pytest.approx(2.43, abs=0.05) + + # Nuclides without fission data have zero nu-fission + assert microxs['O16', 'nu-fission'][0] == 0.0 + + # A fast flux should produce a larger nu-bar + flux = [0., 0., 0., 1.0] + microxs_fast = MicroXS.from_multigroup_flux( + energies=energies, multigroup_flux=flux, chain_file=chain_file, + nuclides=['U235'], reactions=['fission', 'nu-fission']) + nu_bar_fast = (microxs_fast['U235', 'nu-fission'][0] + / microxs_fast['U235', 'fission'][0]) + assert nu_bar_fast > nu_bar + + def test_microxs_zero_flux(): chain_file = Path(__file__).parents[1] / 'chain_simple.xml' From 3449f4e0f5f0273f5f7b849916f104e455449806 Mon Sep 17 00:00:00 2001 From: Eden Rochman Date: Wed, 16 Sep 2026 12:23:07 +0200 Subject: [PATCH 2/6] Refactor k-infinity estimator per review suggestions - Auto-detect kinf computation when keff is not provided and all MicroXS contain fission + nu-fission, removing the calculate_kinf parameter from IndependentOperator and from_nuclides - Add neutrons_out field to ReactionInfo in chain.py, replacing the _XN_REACTION regex and _neutrons_emitted() function with a lookup into the REACTIONS dict for neutron exit-channel counts - Move nu-fission cross section collapsing from Python to C++ via Nuclide::collapse_nu_fission_rate, with ctypes binding and C API --- docs/source/usersguide/depletion.rst | 11 +- include/openmc/nuclide.h | 13 ++ openmc/deplete/chain.py | 172 +++++++++--------- openmc/deplete/independent_operator.py | 77 +++----- openmc/deplete/microxs.py | 101 +--------- openmc/lib/nuclide.py | 32 ++++ src/nuclide.cpp | 108 +++++++++++ .../test_deplete_independent_operator.py | 54 +++--- 8 files changed, 308 insertions(+), 260 deletions(-) diff --git a/docs/source/usersguide/depletion.rst b/docs/source/usersguide/depletion.rst index 0459ade02d6..8c70a1a734d 100644 --- a/docs/source/usersguide/depletion.rst +++ b/docs/source/usersguide/depletion.rst @@ -270,13 +270,12 @@ transport-depletion calculation and follow the same steps from there. the depletion chain with at least one reaction, that reaction will not be simulated. -If the microscopic cross section data includes 'fission' and 'nu-fission' -cross sections, :class:`~openmc.deplete.IndependentOperator` can also estimate -the infinite multiplication factor at each depletion step by passing -``calculate_kinf=True``:: +When the microscopic cross section data includes both 'fission' and +'nu-fission' cross sections and no explicit ``keff`` value is provided, +:class:`~openmc.deplete.IndependentOperator` automatically estimates the +infinite multiplication factor at each depletion step:: - op = openmc.deplete.IndependentOperator(materials, fluxes, micros, - calculate_kinf=True) + op = openmc.deplete.IndependentOperator(materials, fluxes, micros) The estimate is computed as the ratio of the neutron production rate to the neutron loss rate based on the one-group reaction rates and is reported as the diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index 7a8b2acadd9..c83c44c19c3 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -84,6 +84,19 @@ class Nuclide { double collapse_rate(int MT, double temperature, span energy, span flux) const; + //! \brief Calculate flux-averaged nu-fission cross section + // + //! Computes the one-group nu(E)*sigma_f(E) collapsed against a multigroup + //! flux, using the same integration scheme as collapse_rate but weighting + //! the fission cross section by the total neutron yield at each energy. + //! + //! \param[in] temperature Temperature in [K] + //! \param[in] energy Energy group boundaries in [eV] + //! \param[in] flux Flux in each energy group (not normalized per eV) + //! \return Flux-averaged nu-fission cross section, or 0.0 if not fissionable + double collapse_nu_fission_rate(double temperature, + span energy, span flux) const; + //============================================================================ // Data members std::string name_; //!< Name of nuclide, e.g. "U235" diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index a2ff6dea0b6..895fd879ce3 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -27,94 +27,94 @@ import openmc.data -# tuple of (possible MT values, secondaries) -ReactionInfo = namedtuple('ReactionInfo', ('mts', 'secondaries')) +# tuple of (possible MT values, secondaries, neutrons emitted in exit channel) +ReactionInfo = namedtuple('ReactionInfo', ('mts', 'secondaries', 'neutrons_out')) REACTIONS = { - '(n,2nd)': ReactionInfo({11}, ('H2',)), - '(n,2n)': ReactionInfo(set(chain([16], range(875, 892))), ()), - '(n,3n)': ReactionInfo({17}, ()), - '(n,na)': ReactionInfo({22}, ('He4',)), - '(n,n3a)': ReactionInfo({23}, ('He4', 'He4', 'He4')), - '(n,2na)': ReactionInfo({24}, ('He4',)), - '(n,3na)': ReactionInfo({25}, ('He4',)), - '(n,np)': ReactionInfo({28}, ('H1',)), - '(n,n2a)': ReactionInfo({29}, ('He4', 'He4')), - '(n,2n2a)': ReactionInfo({30}, ('He4', 'He4')), - '(n,nd)': ReactionInfo({32}, ('H2',)), - '(n,nt)': ReactionInfo({33}, ('H3',)), - '(n,n3He)': ReactionInfo({34}, ('He3',)), - '(n,nd2a)': ReactionInfo({35}, ('H2', 'He4', 'He4')), - '(n,nt2a)': ReactionInfo({36}, ('H3', 'He4', 'He4')), - '(n,4n)': ReactionInfo({37}, ()), - '(n,2np)': ReactionInfo({41}, ('H1',)), - '(n,3np)': ReactionInfo({42}, ('H1',)), - '(n,n2p)': ReactionInfo({44}, ('H1', 'H1')), - '(n,npa)': ReactionInfo({45}, ('H1', 'He4')), - '(n,gamma)': ReactionInfo({102}, ()), - '(n,p)': ReactionInfo(set(chain([103], range(600, 650))), ('H1',)), - '(n,d)': ReactionInfo(set(chain([104], range(650, 700))), ('H2',)), - '(n,t)': ReactionInfo(set(chain([105], range(700, 750))), ('H3',)), - '(n,3He)': ReactionInfo(set(chain([106], range(750, 800))), ('He3',)), - '(n,a)': ReactionInfo(set(chain([107], range(800, 850))), ('He4',)), - '(n,2a)': ReactionInfo({108}, ('He4', 'He4')), - '(n,3a)': ReactionInfo({109}, ('He4', 'He4', 'He4')), - '(n,2p)': ReactionInfo({111}, ('H1', 'H1')), - '(n,pa)': ReactionInfo({112}, ('H1', 'He4')), - '(n,t2a)': ReactionInfo({113}, ('H3', 'He4', 'He4')), - '(n,d2a)': ReactionInfo({114}, ('H2', 'He4', 'He4')), - '(n,pd)': ReactionInfo({115}, ('H1', 'H2')), - '(n,pt)': ReactionInfo({116}, ('H1', 'H3')), - '(n,da)': ReactionInfo({117}, ('H2', 'He4')), - '(n,5n)': ReactionInfo({152}, ()), - '(n,6n)': ReactionInfo({153}, ()), - '(n,2nt)': ReactionInfo({154}, ('H3',)), - '(n,ta)': ReactionInfo({155}, ('H3', 'He4')), - '(n,4np)': ReactionInfo({156}, ('H1',)), - '(n,3nd)': ReactionInfo({157}, ('H2',)), - '(n,nda)': ReactionInfo({158}, ('H2', 'He4')), - '(n,2npa)': ReactionInfo({159}, ('H1', 'He4')), - '(n,7n)': ReactionInfo({160}, ()), - '(n,8n)': ReactionInfo({161}, ()), - '(n,5np)': ReactionInfo({162}, ('H1',)), - '(n,6np)': ReactionInfo({163}, ('H1',)), - '(n,7np)': ReactionInfo({164}, ('H1',)), - '(n,4na)': ReactionInfo({165}, ('He4',)), - '(n,5na)': ReactionInfo({166}, ('He4',)), - '(n,6na)': ReactionInfo({167}, ('He4',)), - '(n,7na)': ReactionInfo({168}, ('He4',)), - '(n,4nd)': ReactionInfo({169}, ('H2',)), - '(n,5nd)': ReactionInfo({170}, ('H2',)), - '(n,6nd)': ReactionInfo({171}, ('H2',)), - '(n,3nt)': ReactionInfo({172}, ('H3',)), - '(n,4nt)': ReactionInfo({173}, ('H3',)), - '(n,5nt)': ReactionInfo({174}, ('H3',)), - '(n,6nt)': ReactionInfo({175}, ('H3',)), - '(n,2n3He)': ReactionInfo({176}, ('He3',)), - '(n,3n3He)': ReactionInfo({177}, ('He3',)), - '(n,4n3He)': ReactionInfo({178}, ('He3',)), - '(n,3n2p)': ReactionInfo({179}, ('H1', 'H1')), - '(n,3n2a)': ReactionInfo({180}, ('He4', 'He4')), - '(n,3npa)': ReactionInfo({181}, ('H1', 'He4')), - '(n,dt)': ReactionInfo({182}, ('H2', 'H3')), - '(n,npd)': ReactionInfo({183}, ('H1', 'H2')), - '(n,npt)': ReactionInfo({184}, ('H1', 'H3')), - '(n,ndt)': ReactionInfo({185}, ('H2', 'H3')), - '(n,np3He)': ReactionInfo({186}, ('H1', 'He3')), - '(n,nd3He)': ReactionInfo({187}, ('H2', 'He3')), - '(n,nt3He)': ReactionInfo({188}, ('H3', 'He3')), - '(n,nta)': ReactionInfo({189}, ('H3', 'He4')), - '(n,2n2p)': ReactionInfo({190}, ('H1', 'H1')), - '(n,p3He)': ReactionInfo({191}, ('H1', 'He3')), - '(n,d3He)': ReactionInfo({192}, ('H2', 'He3')), - '(n,3Hea)': ReactionInfo({193}, ('He3', 'He4')), - '(n,4n2p)': ReactionInfo({194}, ('H1', 'H1')), - '(n,4n2a)': ReactionInfo({195}, ('He4', 'He4')), - '(n,4npa)': ReactionInfo({196}, ('H1', 'He4')), - '(n,3p)': ReactionInfo({197}, ('H1', 'H1', 'H1')), - '(n,n3p)': ReactionInfo({198}, ('H1', 'H1', 'H1')), - '(n,3n2pa)': ReactionInfo({199}, ('H1', 'H1', 'He4')), - '(n,5n2p)': ReactionInfo({200}, ('H1', 'H1')), + '(n,2nd)': ReactionInfo({11}, ('H2',), 2), + '(n,2n)': ReactionInfo(set(chain([16], range(875, 892))), (), 2), + '(n,3n)': ReactionInfo({17}, (), 3), + '(n,na)': ReactionInfo({22}, ('He4',), 1), + '(n,n3a)': ReactionInfo({23}, ('He4', 'He4', 'He4'), 1), + '(n,2na)': ReactionInfo({24}, ('He4',), 2), + '(n,3na)': ReactionInfo({25}, ('He4',), 3), + '(n,np)': ReactionInfo({28}, ('H1',), 1), + '(n,n2a)': ReactionInfo({29}, ('He4', 'He4'), 1), + '(n,2n2a)': ReactionInfo({30}, ('He4', 'He4'), 2), + '(n,nd)': ReactionInfo({32}, ('H2',), 1), + '(n,nt)': ReactionInfo({33}, ('H3',), 1), + '(n,n3He)': ReactionInfo({34}, ('He3',), 1), + '(n,nd2a)': ReactionInfo({35}, ('H2', 'He4', 'He4'), 1), + '(n,nt2a)': ReactionInfo({36}, ('H3', 'He4', 'He4'), 1), + '(n,4n)': ReactionInfo({37}, (), 4), + '(n,2np)': ReactionInfo({41}, ('H1',), 2), + '(n,3np)': ReactionInfo({42}, ('H1',), 3), + '(n,n2p)': ReactionInfo({44}, ('H1', 'H1'), 1), + '(n,npa)': ReactionInfo({45}, ('H1', 'He4'), 1), + '(n,gamma)': ReactionInfo({102}, (), 0), + '(n,p)': ReactionInfo(set(chain([103], range(600, 650))), ('H1',), 0), + '(n,d)': ReactionInfo(set(chain([104], range(650, 700))), ('H2',), 0), + '(n,t)': ReactionInfo(set(chain([105], range(700, 750))), ('H3',), 0), + '(n,3He)': ReactionInfo(set(chain([106], range(750, 800))), ('He3',), 0), + '(n,a)': ReactionInfo(set(chain([107], range(800, 850))), ('He4',), 0), + '(n,2a)': ReactionInfo({108}, ('He4', 'He4'), 0), + '(n,3a)': ReactionInfo({109}, ('He4', 'He4', 'He4'), 0), + '(n,2p)': ReactionInfo({111}, ('H1', 'H1'), 0), + '(n,pa)': ReactionInfo({112}, ('H1', 'He4'), 0), + '(n,t2a)': ReactionInfo({113}, ('H3', 'He4', 'He4'), 0), + '(n,d2a)': ReactionInfo({114}, ('H2', 'He4', 'He4'), 0), + '(n,pd)': ReactionInfo({115}, ('H1', 'H2'), 0), + '(n,pt)': ReactionInfo({116}, ('H1', 'H3'), 0), + '(n,da)': ReactionInfo({117}, ('H2', 'He4'), 0), + '(n,5n)': ReactionInfo({152}, (), 5), + '(n,6n)': ReactionInfo({153}, (), 6), + '(n,2nt)': ReactionInfo({154}, ('H3',), 2), + '(n,ta)': ReactionInfo({155}, ('H3', 'He4'), 0), + '(n,4np)': ReactionInfo({156}, ('H1',), 4), + '(n,3nd)': ReactionInfo({157}, ('H2',), 3), + '(n,nda)': ReactionInfo({158}, ('H2', 'He4'), 1), + '(n,2npa)': ReactionInfo({159}, ('H1', 'He4'), 2), + '(n,7n)': ReactionInfo({160}, (), 7), + '(n,8n)': ReactionInfo({161}, (), 8), + '(n,5np)': ReactionInfo({162}, ('H1',), 5), + '(n,6np)': ReactionInfo({163}, ('H1',), 6), + '(n,7np)': ReactionInfo({164}, ('H1',), 7), + '(n,4na)': ReactionInfo({165}, ('He4',), 4), + '(n,5na)': ReactionInfo({166}, ('He4',), 5), + '(n,6na)': ReactionInfo({167}, ('He4',), 6), + '(n,7na)': ReactionInfo({168}, ('He4',), 7), + '(n,4nd)': ReactionInfo({169}, ('H2',), 4), + '(n,5nd)': ReactionInfo({170}, ('H2',), 5), + '(n,6nd)': ReactionInfo({171}, ('H2',), 6), + '(n,3nt)': ReactionInfo({172}, ('H3',), 3), + '(n,4nt)': ReactionInfo({173}, ('H3',), 4), + '(n,5nt)': ReactionInfo({174}, ('H3',), 5), + '(n,6nt)': ReactionInfo({175}, ('H3',), 6), + '(n,2n3He)': ReactionInfo({176}, ('He3',), 2), + '(n,3n3He)': ReactionInfo({177}, ('He3',), 3), + '(n,4n3He)': ReactionInfo({178}, ('He3',), 4), + '(n,3n2p)': ReactionInfo({179}, ('H1', 'H1'), 3), + '(n,3n2a)': ReactionInfo({180}, ('He4', 'He4'), 3), + '(n,3npa)': ReactionInfo({181}, ('H1', 'He4'), 3), + '(n,dt)': ReactionInfo({182}, ('H2', 'H3'), 0), + '(n,npd)': ReactionInfo({183}, ('H1', 'H2'), 1), + '(n,npt)': ReactionInfo({184}, ('H1', 'H3'), 1), + '(n,ndt)': ReactionInfo({185}, ('H2', 'H3'), 1), + '(n,np3He)': ReactionInfo({186}, ('H1', 'He3'), 1), + '(n,nd3He)': ReactionInfo({187}, ('H2', 'He3'), 1), + '(n,nt3He)': ReactionInfo({188}, ('H3', 'He3'), 1), + '(n,nta)': ReactionInfo({189}, ('H3', 'He4'), 1), + '(n,2n2p)': ReactionInfo({190}, ('H1', 'H1'), 2), + '(n,p3He)': ReactionInfo({191}, ('H1', 'He3'), 0), + '(n,d3He)': ReactionInfo({192}, ('H2', 'He3'), 0), + '(n,3Hea)': ReactionInfo({193}, ('He3', 'He4'), 0), + '(n,4n2p)': ReactionInfo({194}, ('H1', 'H1'), 4), + '(n,4n2a)': ReactionInfo({195}, ('He4', 'He4'), 4), + '(n,4npa)': ReactionInfo({196}, ('H1', 'He4'), 4), + '(n,3p)': ReactionInfo({197}, ('H1', 'H1', 'H1'), 0), + '(n,n3p)': ReactionInfo({198}, ('H1', 'H1', 'H1'), 1), + '(n,3n2pa)': ReactionInfo({199}, ('H1', 'H1', 'He4'), 3), + '(n,5n2p)': ReactionInfo({200}, ('H1', 'H1'), 5), } __all__ = ["Chain", "REACTIONS"] diff --git a/openmc/deplete/independent_operator.py b/openmc/deplete/independent_operator.py index 42eb92c0b57..55630beca01 100644 --- a/openmc/deplete/independent_operator.py +++ b/openmc/deplete/independent_operator.py @@ -8,7 +8,6 @@ from __future__ import annotations from collections.abc import Iterable import copy -import re import numpy as np from uncertainties import ufloat @@ -17,24 +16,13 @@ from openmc.checkvalue import check_type from openmc.mpi import comm from .abc import ReactionRateHelper, OperatorResult +from .chain import REACTIONS from .openmc_operator import OpenMCOperator from .pool import _distribute from .microxs import MicroXS from .results import Results from .helpers import ChainFissionHelper, ConstantFissionYieldHelper, SourceRateHelper -# Regular expression matching reactions that emit one or more neutrons, e.g., -# (n,2n) or (n,np), with the number of emitted neutrons captured -_XN_REACTION = re.compile(r'\(n,(\d*)n') - - -def _neutrons_emitted(reaction: str) -> int: - """Number of neutrons in the exit channel of a transmutation reaction.""" - match = _XN_REACTION.match(reaction) - if match is None: - return 0 - return int(match.group(1)) if match.group(1) else 1 - class IndependentOperator(OpenMCOperator): """Transport-independent transport operator based on multigroup data. @@ -67,15 +55,15 @@ class IndependentOperator(OpenMCOperator): Path to the depletion chain XML file or instance of openmc.deplete.Chain. Defaults to ``openmc.config['chain_file']``. keff : 2-tuple of float, optional - keff eigenvalue and uncertainty from transport calculation. - calculate_kinf : bool, optional - If True, the infinite multiplication factor is estimated from the - material compositions and one-group cross sections at each depletion - step and reported as the eigenvalue in the depletion results. Requires - that each :class:`~openmc.deplete.MicroXS` instance contains 'fission' - and 'nu-fission' cross sections. Mutually exclusive with ``keff``. - - .. versionadded:: 0.15.4 + keff eigenvalue and uncertainty from transport calculation. When not + provided and every :class:`~openmc.deplete.MicroXS` instance contains + both 'fission' and 'nu-fission' cross sections, the infinite + multiplication factor is estimated automatically from the material + compositions and one-group cross sections at each depletion step. + + .. versionchanged:: 0.15.4 + k-infinity is now estimated automatically when ``keff`` is not + given and the required cross sections are present. prev_results : Results, optional Results from a previous depletion calculation. normalization_mode : {"fission-q", "source-rate"} @@ -137,8 +125,7 @@ def __init__(self, fission_q=None, prev_results=None, reduce_chain_level=None, - fission_yield_opts=None, - calculate_kinf=False): + fission_yield_opts=None): # Validate micro-xs parameters check_type('materials', materials, Iterable, openmc.Material) check_type('micros', micros, Iterable, MicroXS) @@ -156,19 +143,13 @@ def __init__(self, self._keff = keff - check_type('calculate_kinf', calculate_kinf, bool) - if calculate_kinf: - if keff is not None: - raise ValueError("The 'keff' and 'calculate_kinf' arguments " - "are mutually exclusive.") - for micro in micros: - missing = {'fission', 'nu-fission'} - set(micro.reactions) - if missing: - raise ValueError( - "Estimating k-infinity requires 'fission' and " - "'nu-fission' cross sections in each MicroXS " - f"instance (missing {sorted(missing)}).") - self._calculate_kinf = calculate_kinf + # Auto-detect k-infinity capability: estimate kinf when keff is not + # provided and every MicroXS contains fission + nu-fission data. + self._calculate_kinf = ( + keff is None + and all('fission' in m.reactions and 'nu-fission' in m.reactions + for m in micros) + ) if fission_yield_opts is None: fission_yield_opts = {} @@ -201,8 +182,7 @@ def from_nuclides(cls, volume, nuclides, fission_q=None, prev_results=None, reduce_chain_level=None, - fission_yield_opts=None, - calculate_kinf=False): + fission_yield_opts=None): """ Alternate constructor from a dictionary of nuclide concentrations @@ -224,13 +204,6 @@ def from_nuclides(cls, volume, nuclides, keff : 2-tuple of float, optional keff eigenvalue and uncertainty from transport calculation. Default is None. - calculate_kinf : bool, optional - If True, the infinite multiplication factor is estimated from the - material compositions and one-group cross sections at each - depletion step. Requires that ``micro_xs`` contains 'fission' and - 'nu-fission' cross sections. Mutually exclusive with ``keff``. - - .. versionadded:: 0.15.4 normalization_mode : {"fission-q", "source-rate"} Indicate how reaction rates should be calculated. ``"fission-q"`` uses the fission Q values from the depletion @@ -266,8 +239,7 @@ def from_nuclides(cls, volume, nuclides, fission_q=fission_q, prev_results=prev_results, reduce_chain_level=reduce_chain_level, - fission_yield_opts=fission_yield_opts, - calculate_kinf=calculate_kinf) + fission_yield_opts=fission_yield_opts) @staticmethod def _consolidate_nuclides_to_material(nuclides, nuc_units, volume): @@ -508,8 +480,13 @@ def _estimate_k_inf(self): rate = atoms * (micro_xs[nuc, rxn] * flux).sum() if rxn == 'nu-fission': production += rate - elif rxn != 'damage-energy': - loss += (1 - _neutrons_emitted(rxn)) * rate + elif rxn == 'damage-energy': + pass + elif rxn in REACTIONS: + n_out = REACTIONS[rxn].neutrons_out + loss += (1 - n_out) * rate + else: + loss += rate # Sum contributions over all MPI processes production = comm.allreduce(production) diff --git a/openmc/deplete/microxs.py b/openmc/deplete/microxs.py index b0279f8defa..4d34a603235 100644 --- a/openmc/deplete/microxs.py +++ b/openmc/deplete/microxs.py @@ -17,7 +17,7 @@ from openmc.checkvalue import check_type, check_value, check_iterable_type, PathLike from openmc import StatePoint from openmc.mgxs import GROUP_STRUCTURES -from openmc.data import DataLibrary, REACTION_MT, Reaction +from openmc.data import REACTION_MT import openmc from .chain import Chain, REACTIONS, _get_chain from .coupled_operator import _find_cross_sections, _get_nuclides_with_data @@ -307,82 +307,6 @@ def get_microxs_and_flux( return fluxes, micros -def _collapse_nu_fission( - path: PathLike, - nuclide: str, - temperature: float, - energies: Sequence[float], - flux: np.ndarray -) -> float: - r"""Compute a one-group fission neutron production cross section. - - The fission neutron production cross section, - :math:`\nu(E)\sigma_f(E)`, is integrated against a flux that is assumed - to be constant in energy within each group, matching the treatment used - for other reactions in :meth:`openmc.lib.Nuclide.collapse_rate`. - - Parameters - ---------- - path : PathLike - Path to the HDF5 data file containing the nuclide. - nuclide : str - Name of the nuclide, e.g., 'U235'. - temperature : float - Temperature in [K]. The closest available temperature is used for the - fission cross section. - energies : iterable of float - Energy group boundaries in [eV] in ascending order. - flux : numpy.ndarray - Flux in each energy group, normalized to sum to unity. - - Returns - ------- - float - Flux-averaged fission neutron production cross section in [b]. Zero if - the nuclide has no fission data. - - """ - with h5py.File(path, 'r') as h5: - group = h5[nuclide] - if 'reactions/reaction_018' not in group: - return 0.0 - - # Select the available temperature closest to the requested one - temp_keys = list(group['energy']) - temps = np.array([float(t[:-1]) for t in temp_keys]) - temp_key = temp_keys[np.argmin(np.abs(temps - temperature))] - - energy_grid = {temp_key: group['energy'][temp_key][()]} - rx = Reaction.from_hdf5(group['reactions/reaction_018'], energy_grid) - - xs = rx.xs[temp_key] - - # Total nu(E) is the sum of the yields of all neutron products. If a - # product with emission mode 'total' is present, use it alone to avoid - # double counting prompt and delayed neutrons. - neutron_products = [p for p in rx.products if p.particle == 'neutron'] - total_products = [p for p in neutron_products if p.emission_mode == 'total'] - if total_products: - neutron_products = total_products - if not neutron_products: - return 0.0 - - def nu(e): - return sum(p.yield_(e) for p in neutron_products) - - # Integrate nu(E)*sigma_f(E) against a histogram flux - nu_fission = 0.0 - for g, flux_g in enumerate(flux): - if flux_g == 0.0: - continue - e_low, e_high = energies[g], energies[g + 1] - inside = xs.x[(xs.x > e_low) & (xs.x < e_high)] - e = np.concatenate([[e_low], inside, [e_high]]) - nu_fission += np.trapezoid(nu(e) * xs(e), e) * flux_g / (e_high - e_low) - - return nu_fission - - class MicroXS: """Microscopic cross section data for use in transport-independent depletion. @@ -505,12 +429,10 @@ def from_multigroup_flux( nuclides = [nuc.name for nuc in nuclides] # Get reaction MT values. If no reactions specified, default to the - # reactions available in the chain file. The 'nu-fission' reaction is - # handled separately since it does not correspond to a single MT value. + # reactions available in the chain file. if reactions is None: reactions = chain.reactions - mts = [REACTION_MT[name] if name != 'nu-fission' else None - for name in reactions] + mts = [REACTION_MT.get(name) for name in reactions] # Create 3D array for microscopic cross sections microxs_arr = np.zeros((len(nuclides), len(mts), 1)) @@ -523,24 +445,19 @@ def from_multigroup_flux( # Normalize multigroup flux multigroup_flux /= flux_sum - # If nu-fission was requested, get paths to pointwise data files - if 'nu-fission' in reactions: - data_library = DataLibrary.from_xml(cross_sections) - # Compute microscopic cross sections within a temporary session with openmc.lib.TemporarySession(**init_kwargs): - # For each nuclide and reaction, compute the flux-averaged xs for nuc_index, nuc in enumerate(nuclides): if nuc not in nuclides_with_data: continue lib_nuc = openmc.lib.load_nuclide(nuc) - for mt_index, mt in enumerate(mts): - if mt is None: - path = data_library.get_by_material(nuc)['path'] + for mt_index, (rxn_name, mt) in enumerate( + zip(reactions, mts)): + if rxn_name == 'nu-fission': microxs_arr[nuc_index, mt_index, 0] = \ - _collapse_nu_fission(path, nuc, temperature, - energies, multigroup_flux) - else: + lib_nuc.collapse_nu_fission_rate( + temperature, energies, multigroup_flux) + elif mt is not None: microxs_arr[nuc_index, mt_index, 0] = \ lib_nuc.collapse_rate( mt, temperature, energies, multigroup_flux) diff --git a/openmc/lib/nuclide.py b/openmc/lib/nuclide.py index ef1287cf34a..3404f6e4bb8 100644 --- a/openmc/lib/nuclide.py +++ b/openmc/lib/nuclide.py @@ -29,6 +29,10 @@ _array_1d_dble, _array_1d_dble, c_int, POINTER(c_double)] _dll.openmc_nuclide_collapse_rate.restype = c_int _dll.openmc_nuclide_collapse_rate.errcheck = _error_handler +_dll.openmc_nuclide_collapse_nu_fission_rate.argtypes = [c_int, c_double, + _array_1d_dble, _array_1d_dble, c_int, POINTER(c_double)] +_dll.openmc_nuclide_collapse_nu_fission_rate.restype = c_int +_dll.openmc_nuclide_collapse_nu_fission_rate.errcheck = _error_handler _dll.nuclides_size.restype = c_size_t @@ -112,6 +116,34 @@ def collapse_rate(self, MT, temperature, energy, flux): flux, len(flux), xs) return xs.value + def collapse_nu_fission_rate(self, temperature, energy, flux): + r"""Calculate flux-averaged :math:`\nu\sigma_f` cross section + + .. versionadded:: 0.15.4 + + Parameters + ---------- + temperature : float + Temperature in [K] at which to evaluate cross sections + energy : iterable of float + Energy group boundaries in [eV] + flux : iterable of float + Flux in each energy group (not normalized per eV) + + Returns + ------- + float + Flux-averaged nu-fission cross section, or 0.0 if nuclide is + not fissionable + + """ + energy = np.asarray(energy, dtype=float) + flux = np.asarray(flux, dtype=float) + xs = c_double() + _dll.openmc_nuclide_collapse_nu_fission_rate( + self._index, temperature, energy, flux, len(flux), xs) + return xs.value + class _NuclideMapping(Mapping): """Provide mapping from nuclide name to index in nuclides array.""" diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 17d6e952c39..00de05554e1 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -1072,6 +1072,95 @@ double Nuclide::collapse_rate(int MT, double temperature, } } +double Nuclide::collapse_nu_fission_rate(double temperature, + span energy, span flux) const +{ + if (!fissionable_) + return 0.0; + + assert(energy.size() > 0); + assert(energy.size() == flux.size() + 1); + + int i_rx = reaction_index_[18]; + if (i_rx < 0) + return 0.0; + const auto& rx = reactions_[i_rx]; + + // Determine temperature index + int64_t i_temp; + double f; + std::tie(i_temp, f) = this->find_temperature(temperature); + + // Helper to integrate nu(E)*sigma_f(E)*phi(E) at a given temperature + auto compute = [&](int64_t t_idx) -> double { + const auto& xs = rx->xs_[t_idx].value; + const auto& grid = grid_[t_idx].energy; + int i_low = lower_bound_index(grid.cbegin(), grid.cend(), energy.front()); + + int j_start = 0; + int i_threshold = rx->xs_[t_idx].threshold; + if (i_low < i_threshold) { + i_low = i_threshold; + while (energy[j_start + 1] < grid[i_low]) { + ++j_start; + if (static_cast(j_start + 1) == energy.size()) + return 0.0; + } + } + + double rate_sum = 0.0; + + for (size_t j = j_start; j < flux.size(); ++j) { + double E_group_low = energy[j]; + double E_group_high = energy[j + 1]; + double flux_per_eV = flux[j] / (E_group_high - E_group_low); + + int i_high = i_low; + while (grid[i_high + 1] < E_group_high && + static_cast(i_high + 1) < grid.size() - 1) + ++i_high; + + for (; i_low <= i_high; ++i_low) { + double E_l = grid[i_low]; + double E_r = grid[i_low + 1]; + if (E_l == E_r) + continue; + + double xs_l = xs[i_low - i_threshold]; + double xs_r = xs[i_low + 1 - i_threshold]; + + double E_low = std::max(E_group_low, E_l); + double E_high = std::min(E_group_high, E_r); + + double m = (xs_r - xs_l) / (E_r - E_l); + double sig_low = xs_l + m * (E_low - E_l); + double sig_high = xs_l + m * (E_high - E_l); + + double nu_low = this->nu(E_low, EmissionMode::total); + double nu_high = this->nu(E_high, EmissionMode::total); + + double nuxs_avg = 0.5 * (nu_low * sig_low + nu_high * sig_high); + + double dE = (E_high - E_low); + rate_sum += flux_per_eV * nuxs_avg * dE; + } + + i_low = i_high; + if (static_cast(i_low + 1) == grid.size()) + break; + } + + return rate_sum; + }; + + double rr_low = compute(i_temp); + if (f > 0.0) { + double rr_high = compute(i_temp + 1); + return rr_low + f * (rr_high - rr_low); + } + return rr_low; +} + //============================================================================== // Non-member functions //============================================================================== @@ -1215,6 +1304,25 @@ extern "C" int openmc_nuclide_collapse_rate(int index, int MT, return 0; } +extern "C" int openmc_nuclide_collapse_nu_fission_rate(int index, + double temperature, const double* energy, const double* flux, int n, + double* xs) +{ + if (index < 0 || index >= data::nuclides.size()) { + set_errmsg("Index in nuclides vector is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + + try { + *xs = data::nuclides[index]->collapse_nu_fission_rate( + temperature, {energy, energy + n + 1}, {flux, flux + n}); + } catch (const std::out_of_range& e) { + set_errmsg(e.what()); + return OPENMC_E_OUT_OF_BOUNDS; + } + return 0; +} + void nuclides_clear() { data::nuclides.clear(); diff --git a/tests/unit_tests/test_deplete_independent_operator.py b/tests/unit_tests/test_deplete_independent_operator.py index 3f682b1d686..d1c0a867c86 100644 --- a/tests/unit_tests/test_deplete_independent_operator.py +++ b/tests/unit_tests/test_deplete_independent_operator.py @@ -9,7 +9,7 @@ from openmc import Material from openmc.deplete import IndependentOperator, MicroXS, Chain -from openmc.deplete.independent_operator import _neutrons_emitted +from openmc.deplete.chain import REACTIONS CHAIN_PATH = Path(__file__).parents[1] / "chain_simple.xml" ONE_GROUP_XS = Path(__file__).parents[1] / "micro_xs_simple.csv" @@ -67,20 +67,19 @@ def _uranium_material(): return fuel -def test_neutrons_emitted(): - assert _neutrons_emitted('fission') == 0 - assert _neutrons_emitted('(n,gamma)') == 0 - assert _neutrons_emitted('(n,p)') == 0 - assert _neutrons_emitted('(n,a)') == 0 - assert _neutrons_emitted('(n,3He)') == 0 - assert _neutrons_emitted('(n,2a)') == 0 - assert _neutrons_emitted('(n,np)') == 1 - assert _neutrons_emitted('(n,nd2a)') == 1 - assert _neutrons_emitted('(n,2n)') == 2 - assert _neutrons_emitted('(n,2nd)') == 2 - assert _neutrons_emitted('(n,3n)') == 3 - assert _neutrons_emitted('(n,3np)') == 3 - assert _neutrons_emitted('(n,4n)') == 4 +def test_neutrons_out_in_reactions(): + assert REACTIONS['(n,gamma)'].neutrons_out == 0 + assert REACTIONS['(n,p)'].neutrons_out == 0 + assert REACTIONS['(n,a)'].neutrons_out == 0 + assert REACTIONS['(n,3He)'].neutrons_out == 0 + assert REACTIONS['(n,2a)'].neutrons_out == 0 + assert REACTIONS['(n,np)'].neutrons_out == 1 + assert REACTIONS['(n,nd2a)'].neutrons_out == 1 + assert REACTIONS['(n,2n)'].neutrons_out == 2 + assert REACTIONS['(n,2nd)'].neutrons_out == 2 + assert REACTIONS['(n,3n)'].neutrons_out == 3 + assert REACTIONS['(n,3np)'].neutrons_out == 3 + assert REACTIONS['(n,4n)'].neutrons_out == 4 def test_calculate_kinf(): @@ -94,9 +93,10 @@ def test_calculate_kinf(): ]) micro_xs = MicroXS(data, nuclides, reactions) + # With fission + nu-fission present and no keff, kinf is auto-detected op = IndependentOperator( [_uranium_material()], [np.array([1.0])], [micro_xs], CHAIN_PATH, - normalization_mode='source-rate', calculate_kinf=True) + normalization_mode='source-rate') vec = op.initial_condition() # Production is nu-fission; loss is fission + (n,gamma) - (n,2n), since @@ -128,7 +128,7 @@ def test_calculate_kinf_multigroup(): op = IndependentOperator( [fuel], [flux], [micro_xs], CHAIN_PATH, - normalization_mode='source-rate', calculate_kinf=True) + normalization_mode='source-rate') vec = op.initial_condition() production = 25.0*0.75 + 120.0*0.25 @@ -137,16 +137,18 @@ def test_calculate_kinf_multigroup(): assert result.k.n == pytest.approx(production / loss) -def test_calculate_kinf_errors(): - # MicroXS without nu-fission data cannot be used to estimate k-infinity +def test_kinf_not_computed_without_nu_fission(): + # When MicroXS lacks nu-fission, kinf is not computed and keff stays None micro_xs = MicroXS.from_csv(ONE_GROUP_XS) - with pytest.raises(ValueError, match="nu-fission"): - IndependentOperator([_uranium_material()], [1.0], [micro_xs], - CHAIN_PATH, calculate_kinf=True) + op = IndependentOperator([_uranium_material()], [1.0], [micro_xs], + CHAIN_PATH) + assert not op._calculate_kinf - # keff and calculate_kinf are mutually exclusive + +def test_kinf_not_computed_when_keff_given(): + # When keff is explicitly given, kinf auto-detection is disabled data = np.zeros((1, 2, 1)) micro_xs = MicroXS(data, ['U235'], ['fission', 'nu-fission']) - with pytest.raises(ValueError, match="mutually exclusive"): - IndependentOperator([_uranium_material()], [1.0], [micro_xs], - CHAIN_PATH, keff=(1.0, 0.0), calculate_kinf=True) + op = IndependentOperator([_uranium_material()], [1.0], [micro_xs], + CHAIN_PATH, keff=(1.0, 0.0)) + assert not op._calculate_kinf From 3649e66fe69965081ad317de18bba9f272fb8aa7 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 16 Sep 2026 19:10:49 +0300 Subject: [PATCH 3/6] some mechanical fixes --- include/openmc/nuclide.h | 4 ++-- openmc/deplete/independent_operator.py | 7 ++++++- openmc/deplete/microxs.py | 7 +++++-- openmc/lib/nuclide.py | 2 +- src/nuclide.cpp | 6 +++--- 5 files changed, 17 insertions(+), 9 deletions(-) diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index 64b3f66d571..141b78f5f8f 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -94,8 +94,8 @@ class Nuclide { //! \param[in] energy Energy group boundaries in [eV] //! \param[in] flux Flux in each energy group (not normalized per eV) //! \return Flux-averaged nu-fission cross section, or 0.0 if not fissionable - double collapse_nu_fission_rate(double temperature, - span energy, span flux) const; + double collapse_nu_fission_rate(double temperature, span energy, + span flux) const; //! Return a ParticleType object representing this nuclide ParticleType particle_type() const { return {Z_, A_, metastable_}; } diff --git a/openmc/deplete/independent_operator.py b/openmc/deplete/independent_operator.py index 7587a3d7f04..33c8ee29964 100644 --- a/openmc/deplete/independent_operator.py +++ b/openmc/deplete/independent_operator.py @@ -8,6 +8,7 @@ from __future__ import annotations from collections.abc import Iterable import copy +from warnings import warn import numpy as np from uncertainties import ufloat @@ -61,7 +62,7 @@ class IndependentOperator(OpenMCOperator): multiplication factor is estimated automatically from the material compositions and one-group cross sections at each depletion step. - .. versionchanged:: 0.15.4 + .. versionchanged:: 0.16.1 k-infinity is now estimated automatically when ``keff`` is not given and the required cross sections are present. prev_results : Results, optional @@ -147,6 +148,7 @@ def __init__(self, # provided and every MicroXS contains fission + nu-fission data. self._calculate_kinf = ( keff is None + and len(micros) > 0 and all('fission' in m.reactions and 'nu-fission' in m.reactions for m in micros) ) @@ -497,6 +499,9 @@ def _estimate_k_inf(self): loss = comm.allreduce(loss) if loss <= 0.0: + warn('Unable to estimate k-infinity because the total neutron ' + 'loss rate is zero. Check that the supplied MicroXS data ' + 'contains absorption reactions for the nuclides present.') return ufloat(0.0, 0.0) return ufloat(production / loss, 0.0) diff --git a/openmc/deplete/microxs.py b/openmc/deplete/microxs.py index b6cc3c1ecd6..61e87e8d61f 100644 --- a/openmc/deplete/microxs.py +++ b/openmc/deplete/microxs.py @@ -437,7 +437,7 @@ def from_multigroup_flux( needed to estimate k-infinity with :class:`~openmc.deplete.IndependentOperator`. - .. versionchanged:: 0.15.4 + .. versionchanged:: 0.16.1 Added support for 'nu-fission'. **init_kwargs : dict Keyword arguments passed to :func:`openmc.lib.init` @@ -474,7 +474,10 @@ def from_multigroup_flux( # reactions available in the chain file. if reactions is None: reactions = chain.reactions - mts = [REACTION_MT.get(name) for name in reactions] + # 'nu-fission' has no MT of its own and is collapsed separately + # below; every other reaction name must map to an MT + mts = [None if name == 'nu-fission' else REACTION_MT[name] + for name in reactions] # Create 3D array for microscopic cross sections microxs_arr = np.zeros((len(nuclides), len(mts), 1)) diff --git a/openmc/lib/nuclide.py b/openmc/lib/nuclide.py index 3404f6e4bb8..d91228aeaa2 100644 --- a/openmc/lib/nuclide.py +++ b/openmc/lib/nuclide.py @@ -119,7 +119,7 @@ def collapse_rate(self, MT, temperature, energy, flux): def collapse_nu_fission_rate(self, temperature, energy, flux): r"""Calculate flux-averaged :math:`\nu\sigma_f` cross section - .. versionadded:: 0.15.4 + .. versionadded:: 0.16.1 Parameters ---------- diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 00de05554e1..17dbd45d113 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -1072,8 +1072,8 @@ double Nuclide::collapse_rate(int MT, double temperature, } } -double Nuclide::collapse_nu_fission_rate(double temperature, - span energy, span flux) const +double Nuclide::collapse_nu_fission_rate( + double temperature, span energy, span flux) const { if (!fissionable_) return 0.0; @@ -1081,7 +1081,7 @@ double Nuclide::collapse_nu_fission_rate(double temperature, assert(energy.size() > 0); assert(energy.size() == flux.size() + 1); - int i_rx = reaction_index_[18]; + int i_rx = reaction_index_[N_FISSION]; if (i_rx < 0) return 0.0; const auto& rx = reactions_[i_rx]; From 3006080a767b8108bb316eb246266002b41377a4 Mon Sep 17 00:00:00 2001 From: Eden Rochman Date: Thu, 17 Sep 2026 11:59:49 +0200 Subject: [PATCH 4/6] Address review: volume fix, precondition, C++ refactor, docstrings - Fix volume normalization in _estimate_k_inf: divide rate by volume_b_cm so multi-material estimates are not biased by material volume - Add precondition after super().__init__() checking that MicroXS covers all chain reactions, disabling kinf estimate with a warning when absorption channels are missing (prevents silent nu-bar result) - Refactor collapse_nu_fission_rate to collapse the pre-tabulated XS_NU_FISSION column on the nuclide union grid instead of duplicating the integration loop with per-segment nu() calls - Rewrite _estimate_k_inf docstring with full formula, index ranges, k-inf vs k-eff distinction, (n,xn) derivation, and assumptions list - Add tests: two-material volume independence and disabled kinf when reactions are missing --- include/openmc/nuclide.h | 5 +- openmc/deplete/independent_operator.py | 98 +++++++++++++++---- src/nuclide.cpp | 46 +++------ .../test_deplete_independent_operator.py | 56 +++++++++++ 4 files changed, 155 insertions(+), 50 deletions(-) diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index 141b78f5f8f..aa64ca6048b 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -87,8 +87,9 @@ class Nuclide { //! \brief Calculate flux-averaged nu-fission cross section // //! Computes the one-group nu(E)*sigma_f(E) collapsed against a multigroup - //! flux, using the same integration scheme as collapse_rate but weighting - //! the fission cross section by the total neutron yield at each energy. + //! flux, using the same integration scheme as collapse_rate but reading + //! the pre-tabulated nu-fission cross section (XS_NU_FISSION), which sums + //! over all partial fission reactions. //! //! \param[in] temperature Temperature in [K] //! \param[in] energy Energy group boundaries in [eV] diff --git a/openmc/deplete/independent_operator.py b/openmc/deplete/independent_operator.py index 33c8ee29964..82c6a4bc8e3 100644 --- a/openmc/deplete/independent_operator.py +++ b/openmc/deplete/independent_operator.py @@ -58,9 +58,11 @@ class IndependentOperator(OpenMCOperator): keff : 2-tuple of float, optional keff eigenvalue and uncertainty from transport calculation. When not provided and every :class:`~openmc.deplete.MicroXS` instance contains - both 'fission' and 'nu-fission' cross sections, the infinite - multiplication factor is estimated automatically from the material - compositions and one-group cross sections at each depletion step. + both 'fission' and 'nu-fission' cross sections as well as all + transmutation reactions defined by the depletion chain, the infinite + multiplication factor of the depletable materials is estimated + automatically from the material compositions and multigroup cross + sections at each depletion step. .. versionchanged:: 0.16.1 k-infinity is now estimated automatically when ``keff`` is not @@ -173,6 +175,25 @@ def __init__(self, helper_kwargs=helper_kwargs, reduce_chain_level=reduce_chain_level) + # The k-infinity estimate divides the neutron production rate by the + # neutron loss rate, so the loss term must include every absorption + # channel that the depletion chain will use. If a MicroXS is missing + # some of the chain's transmutation reactions (e.g., only 'fission' + # and 'nu-fission' were tallied), the ratio would silently degenerate + # toward nu-bar rather than k-infinity. Note that self.chain only + # exists after the super().__init__() call above (which also applies + # any chain reduction), so this check must come here. + if self._calculate_kinf: + chain_rxns = set(self.chain.reactions) + for m in micros: + if not chain_rxns <= set(m.reactions): + missing = chain_rxns - set(m.reactions) + warn(f'Disabling k-infinity estimate: MicroXS is missing ' + f'chain reactions {missing}. The estimate requires ' + f'all absorption channels to be present.') + self._calculate_kinf = False + break + @classmethod def from_nuclides(cls, volume, nuclides, flux, @@ -446,23 +467,60 @@ def __call__(self, vec, source_rate) -> OperatorResult: return copy.deepcopy(op_result) def _estimate_k_inf(self): - r"""Estimate the infinite multiplication factor. + r"""Estimate the infinite multiplication factor of the depletable + materials. The estimate is computed as the ratio of the neutron production rate - to the neutron loss rate: + to the neutron loss rate summed over the *depletable materials only*: + + .. math:: + k_\infty = \frac{\displaystyle\sum_m \frac{1}{V_m} \sum_i N_{m,i} + \sum_g (\nu\sigma_f)_{m,i,g}\, \phi_{m,g}} + {\displaystyle\sum_m \frac{1}{V_m} \sum_i N_{m,i} + \sum_j (1 - x_j) \sum_g \sigma_{m,i,j,g}\, + \phi_{m,g}} + + where the index :math:`m` runs over the depletable materials, + :math:`i` over the nuclides with cross-section data, :math:`j` over + the transmutation reactions, and :math:`g` over the energy groups. + :math:`N_{m,i}` is the number of atoms of nuclide :math:`i` in + material :math:`m`, :math:`V_m` is the material volume, + :math:`\phi_{m,g}` is the volume-integrated multigroup flux from the + transport run, :math:`(\nu\sigma_f)_{m,i,g}` is the fission neutron + production cross section, :math:`\sigma_{m,i,j,g}` is the cross + section of transmutation reaction :math:`j`, and :math:`x_j` is the + number of neutrons emitted by reaction :math:`j`. + + **This is not k-eff.** The balance above contains no leakage term, so + it relates to the effective multiplication factor as + :math:`k_\infty = k_\text{eff} / (1 - L)` where :math:`L` is the + leakage fraction. Moreover, only depletable materials contribute to + the loss term: for models that also contain non-depletable materials + (moderator, cladding, reflector, ...), absorption in those materials + is not accounted for, and the estimate will be *higher* than the true + k-infinity of the full system. In other words, the estimate assumes + that all relevant absorption happens in the depletable materials. + + The treatment of (n,xn) reactions follows from writing the + multiplication factor as .. math:: - k_\infty = \frac{\sum_i N_i (\nu\sigma_f)_i} - {\sum_i N_i \sum_j (1 - x_j) \sigma_{i,j}} - - where :math:`N_i` is the number of atoms of nuclide :math:`i`, - :math:`(\nu\sigma_f)_i` is its one-group fission neutron production - cross section, :math:`\sigma_{i,j}` is the one-group cross section of - transmutation reaction :math:`j`, and :math:`x_j` is the number of - neutrons emitted by reaction :math:`j`. This is consistent with the - definition of the multiplication factor used elsewhere in OpenMC: - neutrons produced in (n,xn) reactions are not counted as production; - instead, each (n,xn) reaction reduces the loss term by :math:`x - 1`. + k_\text{eff} = \frac{P}{A + L - X} + + where :math:`P` is the fission neutron production rate, :math:`A` the + absorption rate, :math:`L` the leakage rate, and :math:`X` the net + neutron production rate from (n,xn) reactions. The denominator uses + "reduced absorption" :math:`A - X`, which is exactly the convention + used by OpenMC's k-eff estimators: neutrons produced in (n,xn) + reactions are not counted as production; instead each (n,xn) reaction + with :math:`x` neutrons out contributes :math:`(1 - x)` times its + rate to the loss term, giving :math:`A - X` in a single pass over the + reactions. + + Assumptions: the multigroup fluxes are those obtained from the + transport run (and are not recomputed as the compositions change), + non-depletable materials do not deplete, and any background + absorption outside the depletable materials is constant and ignored. Returns ------- @@ -476,6 +534,11 @@ def _estimate_k_inf(self): i_mat = self._mat_index_map[mat] flux = self.fluxes[i_mat] micro_xs = self.cross_sections[i_mat] + + # Convert total atoms and volume-integrated flux to rates per + # unit volume, consistent with _calculate_reaction_rates + volume_b_cm = 1e24 * self.number.get_mat_volume(mat) + for nuc in micro_xs.nuclides: if nuc not in self.number.index_nuc: continue @@ -483,7 +546,8 @@ def _estimate_k_inf(self): if atoms <= 0.0: continue for rxn in micro_xs.reactions: - rate = atoms * (micro_xs[nuc, rxn] * flux).sum() + rate = (atoms * (micro_xs[nuc, rxn] * flux).sum() + / volume_b_cm) if rxn == 'nu-fission': production += rate elif rxn == 'damage-energy': diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 17dbd45d113..5b8996d0f42 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -1081,36 +1081,25 @@ double Nuclide::collapse_nu_fission_rate( assert(energy.size() > 0); assert(energy.size() == flux.size() + 1); - int i_rx = reaction_index_[N_FISSION]; - if (i_rx < 0) - return 0.0; - const auto& rx = reactions_[i_rx]; - // Determine temperature index int64_t i_temp; double f; std::tie(i_temp, f) = this->find_temperature(temperature); - // Helper to integrate nu(E)*sigma_f(E)*phi(E) at a given temperature + // Helper to integrate nu(E)*sigma_f(E)*phi(E) at a given temperature. The + // nu-fission cross section is pre-tabulated on the nuclide's union energy + // grid (see Nuclide::init_grid), so we collapse the XS_NU_FISSION column + // directly with the same trapezoidal scheme as Reaction::collapse_rate. + // Unlike a Reaction's cross section array, the nuclide grid has no + // threshold offset, and XS_NU_FISSION already sums over all partial + // fission reactions. auto compute = [&](int64_t t_idx) -> double { - const auto& xs = rx->xs_[t_idx].value; const auto& grid = grid_[t_idx].energy; int i_low = lower_bound_index(grid.cbegin(), grid.cend(), energy.front()); - int j_start = 0; - int i_threshold = rx->xs_[t_idx].threshold; - if (i_low < i_threshold) { - i_low = i_threshold; - while (energy[j_start + 1] < grid[i_low]) { - ++j_start; - if (static_cast(j_start + 1) == energy.size()) - return 0.0; - } - } - double rate_sum = 0.0; - for (size_t j = j_start; j < flux.size(); ++j) { + for (size_t j = 0; j < flux.size(); ++j) { double E_group_low = energy[j]; double E_group_high = energy[j + 1]; double flux_per_eV = flux[j] / (E_group_high - E_group_low); @@ -1126,23 +1115,18 @@ double Nuclide::collapse_nu_fission_rate( if (E_l == E_r) continue; - double xs_l = xs[i_low - i_threshold]; - double xs_r = xs[i_low + 1 - i_threshold]; + double nuxs_l = xs_[t_idx](i_low, XS_NU_FISSION); + double nuxs_r = xs_[t_idx](i_low + 1, XS_NU_FISSION); double E_low = std::max(E_group_low, E_l); double E_high = std::min(E_group_high, E_r); - double m = (xs_r - xs_l) / (E_r - E_l); - double sig_low = xs_l + m * (E_low - E_l); - double sig_high = xs_l + m * (E_high - E_l); - - double nu_low = this->nu(E_low, EmissionMode::total); - double nu_high = this->nu(E_high, EmissionMode::total); - - double nuxs_avg = 0.5 * (nu_low * sig_low + nu_high * sig_high); + double m = (nuxs_r - nuxs_l) / (E_r - E_l); + double sig_low = nuxs_l + m * (E_low - E_l); + double sig_high = nuxs_l + m * (E_high - E_l); + double nuxs_avg = 0.5 * (sig_low + sig_high); - double dE = (E_high - E_low); - rate_sum += flux_per_eV * nuxs_avg * dE; + rate_sum += flux_per_eV * nuxs_avg * (E_high - E_low); } i_low = i_high; diff --git a/tests/unit_tests/test_deplete_independent_operator.py b/tests/unit_tests/test_deplete_independent_operator.py index d1c0a867c86..abd0e52f13a 100644 --- a/tests/unit_tests/test_deplete_independent_operator.py +++ b/tests/unit_tests/test_deplete_independent_operator.py @@ -137,6 +137,62 @@ def test_calculate_kinf_multigroup(): assert result.k.n == pytest.approx(production / loss) +def test_calculate_kinf_two_materials(): + # Two depletable materials with different volumes. Material 1 contains + # only U235 (the fissile inventory) and material 2 contains only U238 + # (a pure absorber here). With proper per-volume normalization the + # estimate depends on the atom densities, not on the material volumes. + nuclides = ['U235', 'U238'] + reactions = ['fission', 'nu-fission', '(n,gamma)'] + data = np.array([ + [[50.0], [120.0], [10.0]], + [[0.0], [0.0], [20.0]], + ]) + micro_xs = MicroXS(data, nuclides, reactions) + + def build_op(vol1, vol2): + mat1 = Material(name="fissile") + mat1.add_nuclide("U235", 0.05) + mat1.set_density("sum") + mat1.depletable = True + mat1.volume = vol1 + + mat2 = Material(name="absorber") + mat2.add_nuclide("U238", 0.05) + mat2.set_density("sum") + mat2.depletable = True + mat2.volume = vol2 + + return IndependentOperator( + [mat1, mat2], [np.array([1.0]), np.array([1.0])], + [micro_xs, micro_xs], CHAIN_PATH, + normalization_mode='source-rate') + + # Both materials have the same atom density (0.05 atom/b-cm), so the + # volume-normalized rates weight both materials equally: + # production = 120, loss = (50 + 10) + 20 = 80 + expected = 120.0 / 80.0 + + for vol1, vol2 in [(1.0, 4.0), (2.0, 2.0), (3.0, 0.5)]: + op = build_op(vol1, vol2) + vec = op.initial_condition() + result = op(vec, 1.0) + assert result.k.n == pytest.approx(expected), (vol1, vol2) + + +def test_kinf_disabled_when_reactions_missing(): + # A MicroXS narrowed to only fission channels lacks the chain's + # absorption reactions (e.g. (n,gamma)); the resulting ratio would just + # be nu-bar, so the estimate must be disabled with a warning. + data = np.array([[[50.0], [120.0]]]) + micro_xs = MicroXS(data, ['U235'], ['fission', 'nu-fission']) + with pytest.warns(UserWarning, match='Disabling k-infinity'): + op = IndependentOperator( + [_uranium_material()], [np.array([1.0])], [micro_xs], CHAIN_PATH, + normalization_mode='source-rate') + assert not op._calculate_kinf + + def test_kinf_not_computed_without_nu_fission(): # When MicroXS lacks nu-fission, kinf is not computed and keff stays None micro_xs = MicroXS.from_csv(ONE_GROUP_XS) From bac36b66779dea8d5dd7c7b7b9208dfbba32944b Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 17 Sep 2026 19:24:54 +0300 Subject: [PATCH 5/6] some more fixes --- openmc/deplete/independent_operator.py | 6 ++- src/nuclide.cpp | 8 +++- tests/unit_tests/test_deplete_microxs.py | 55 ++++++++++++++++++++++++ 3 files changed, 66 insertions(+), 3 deletions(-) diff --git a/openmc/deplete/independent_operator.py b/openmc/deplete/independent_operator.py index 82c6a4bc8e3..6231759f601 100644 --- a/openmc/deplete/independent_operator.py +++ b/openmc/deplete/independent_operator.py @@ -519,8 +519,10 @@ def _estimate_k_inf(self): Assumptions: the multigroup fluxes are those obtained from the transport run (and are not recomputed as the compositions change), - non-depletable materials do not deplete, and any background - absorption outside the depletable materials is constant and ignored. + non-depletable materials do not deplete, and absorption outside the + depletable materials is ignored entirely. The bias introduced by that + last assumption stays roughly constant over the depletion only to the + extent that the background absorption itself does. Returns ------- diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 5b8996d0f42..c3fd509fa48 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -1095,7 +1095,13 @@ double Nuclide::collapse_nu_fission_rate( // fission reactions. auto compute = [&](int64_t t_idx) -> double { const auto& grid = grid_[t_idx].energy; - int i_low = lower_bound_index(grid.cbegin(), grid.cend(), energy.front()); + // lower_bound_index returns -1 when the first group boundary lies below + // the nuclide's energy grid, which is the usual case for a group + // structure starting at 0 eV. Reaction::collapse_rate happens to mask + // this with its threshold adjustment; the nuclide grid has no threshold, + // so clamp explicitly rather than indexing grid[-1] below. + int i_low = std::max(0, static_cast(lower_bound_index( + grid.cbegin(), grid.cend(), energy.front()))); double rate_sum = 0.0; diff --git a/tests/unit_tests/test_deplete_microxs.py b/tests/unit_tests/test_deplete_microxs.py index 2b78a3b1420..c7a4a8f5afa 100644 --- a/tests/unit_tests/test_deplete_microxs.py +++ b/tests/unit_tests/test_deplete_microxs.py @@ -143,6 +143,61 @@ def test_from_multigroup_flux_nu_fission(): assert nu_bar_fast > nu_bar +def test_nu_fission_collapse_below_energy_grid(): + """Group structure starting below the nuclide's energy grid. + + Regression test: when the first group boundary lies below the nuclide's + first grid energy, which is the usual case for a structure starting at + 0 eV, the index lookup used to collapse nu-fission returns -1 and the + integration must not read before the start of the grid. + + U238 is the interesting case because it is fissionable, so the collapse + runs the full integration loop, but its fission threshold is around + 1 MeV. Against a purely thermal flux the answer is therefore exactly + zero, and any spurious contribution from below the grid shows up as a + nonzero result. O16 does not exercise this path: it is not fissionable + and returns early. + """ + chain_file = Path(__file__).parents[1] / 'chain_simple.xml' + energies = [0., 6.25e-1, 5.53e3, 8.21e5, 2.e7] + flux = [1.0, 0., 0., 0.] + + microxs = MicroXS.from_multigroup_flux( + energies=energies, multigroup_flux=flux, chain_file=chain_file, + nuclides=['U238'], reactions=['fission', 'nu-fission']) + + assert microxs['U238', 'fission'][0] == 0.0 + assert microxs['U238', 'nu-fission'][0] == 0.0 + + +def test_nu_fission_collapse_lower_bound_invariance(): + """Extending the lowest group below the grid must not change the result. + + There is no cross section data below the nuclide's first grid energy, so + collapsing over [0, 0.625] and over [1e-5, 0.625] with the same flux per + eV must give the same answer. Assumes the neutron grid starts at or below + 1e-5 eV, which holds for the standard libraries. + """ + chain_file = Path(__file__).parents[1] / 'chain_simple.xml' + upper = [6.25e-1, 5.53e3, 8.21e5, 2.e7] + + # Both cases give a flux per eV of exactly 1.0 in the lowest group, so a + # correct implementation integrates an identical set of segments. + from_zero = MicroXS.from_multigroup_flux( + energies=[0.] + upper, multigroup_flux=[6.25e-1, 0., 0., 0.], + chain_file=chain_file, nuclides=['U235'], + reactions=['fission', 'nu-fission']) + from_grid = MicroXS.from_multigroup_flux( + energies=[1.0e-5] + upper, + multigroup_flux=[6.25e-1 - 1.0e-5, 0., 0., 0.], + chain_file=chain_file, nuclides=['U235'], + reactions=['fission', 'nu-fission']) + + for rxn in ('fission', 'nu-fission'): + assert from_zero['U235', rxn][0] == pytest.approx( + from_grid['U235', rxn][0], rel=1e-6) + + def test_microxs_zero_flux(): chain_file = Path(__file__).parents[1] / 'chain_simple.xml' From b83ba271b4b3124536e3fa1e21e4fbd2d2df51e3 Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 17 Sep 2026 19:56:03 +0300 Subject: [PATCH 6/6] fix the test --- tests/unit_tests/test_deplete_microxs.py | 71 ++++++++---------------- 1 file changed, 23 insertions(+), 48 deletions(-) diff --git a/tests/unit_tests/test_deplete_microxs.py b/tests/unit_tests/test_deplete_microxs.py index c7a4a8f5afa..8a0815b2f5e 100644 --- a/tests/unit_tests/test_deplete_microxs.py +++ b/tests/unit_tests/test_deplete_microxs.py @@ -143,59 +143,34 @@ def test_from_multigroup_flux_nu_fission(): assert nu_bar_fast > nu_bar -def test_nu_fission_collapse_below_energy_grid(): - """Group structure starting below the nuclide's energy grid. - - Regression test: when the first group boundary lies below the nuclide's - first grid energy, which is the usual case for a structure starting at - 0 eV, the index lookup used to collapse nu-fission returns -1 and the - integration must not read before the start of the grid. - - U238 is the interesting case because it is fissionable, so the collapse - runs the full integration loop, but its fission threshold is around - 1 MeV. Against a purely thermal flux the answer is therefore exactly - zero, and any spurious contribution from below the grid shows up as a - nonzero result. O16 does not exercise this path: it is not fissionable - and returns early. +def test_nu_fission_collapse_lower_bound(): + """Collapsed nu-bar must not depend on the lowest group boundary. + + ``from_multigroup_flux`` normalizes the flux and returns a flux-averaged + cross section, so moving the lowest boundary changes the group width and + rescales both channels. Taking the ratio cancels that, leaving a quantity + that is genuinely invariant: there is no cross section data below the + nuclide's first grid energy, so both structures integrate the same range. + + This guards the index handling in ``Nuclide::collapse_nu_fission_rate``. + The energy lookup returns -1 when the first boundary lies below the grid, + and unlike ``Reaction::collapse_rate`` that function has no threshold + adjustment to mask it. The fission channel goes through the guarded path, + so a mismatch between the two shows up in the ratio. """ chain_file = Path(__file__).parents[1] / 'chain_simple.xml' - energies = [0., 6.25e-1, 5.53e3, 8.21e5, 2.e7] + upper = [6.25e-1, 5.53e3, 8.21e5, 2.e7] flux = [1.0, 0., 0., 0.] - microxs = MicroXS.from_multigroup_flux( - energies=energies, multigroup_flux=flux, chain_file=chain_file, - nuclides=['U238'], reactions=['fission', 'nu-fission']) - - assert microxs['U238', 'fission'][0] == 0.0 - assert microxs['U238', 'nu-fission'][0] == 0.0 - - -def test_nu_fission_collapse_lower_bound_invariance(): - """Extending the lowest group below the grid must not change the result. - - There is no cross section data below the nuclide's first grid energy, so - collapsing over [0, 0.625] and over [1e-5, 0.625] with the same flux per - eV must give the same answer. Assumes the neutron grid starts at or below - 1e-5 eV, which holds for the standard libraries. - """ - chain_file = Path(__file__).parents[1] / 'chain_simple.xml' - upper = [6.25e-1, 5.53e3, 8.21e5, 2.e7] + def nu_bar(first_boundary): + microxs = MicroXS.from_multigroup_flux( + energies=[first_boundary] + upper, multigroup_flux=flux, + chain_file=chain_file, nuclides=['U235'], + reactions=['fission', 'nu-fission']) + return (microxs['U235', 'nu-fission'][0] + / microxs['U235', 'fission'][0]) - # Both cases give a flux per eV of exactly 1.0 in the lowest group, so a - # correct implementation integrates an identical set of segments. - from_zero = MicroXS.from_multigroup_flux( - energies=[0.] + upper, multigroup_flux=[6.25e-1, 0., 0., 0.], - chain_file=chain_file, nuclides=['U235'], - reactions=['fission', 'nu-fission']) - from_grid = MicroXS.from_multigroup_flux( - energies=[1.0e-5] + upper, - multigroup_flux=[6.25e-1 - 1.0e-5, 0., 0., 0.], - chain_file=chain_file, nuclides=['U235'], - reactions=['fission', 'nu-fission']) - - for rxn in ('fission', 'nu-fission'): - assert from_zero['U235', rxn][0] == pytest.approx( - from_grid['U235', rxn][0], rel=1e-6) + assert nu_bar(0.0) == pytest.approx(nu_bar(1.0e-5), rel=1e-6) def test_microxs_zero_flux():