diff --git a/docs/source/usersguide/depletion.rst b/docs/source/usersguide/depletion.rst index 4e52fefde6a..e6736d81edf 100644 --- a/docs/source/usersguide/depletion.rst +++ b/docs/source/usersguide/depletion.rst @@ -270,6 +270,21 @@ transport-depletion calculation and follow the same steps from there. the depletion chain with at least one reaction, that reaction will not be simulated. +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) + +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/include/openmc/nuclide.h b/include/openmc/nuclide.h index ae39a53ddf5..aa64ca6048b 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -84,6 +84,20 @@ 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 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] + //! \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; + //! Return a ParticleType object representing this nuclide ParticleType particle_type() const { return {Z_, A_, metastable_}; } diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 9a35654a47e..fa3ba91557b 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 bdfc32763be..6231759f601 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 @@ -16,6 +17,7 @@ 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 @@ -54,7 +56,17 @@ 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. + 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 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 + given and the required cross sections are present. prev_results : Results, optional Results from a previous depletion calculation. normalization_mode : {"fission-q", "source-rate"} @@ -134,6 +146,15 @@ def __init__(self, self._keff = keff + # 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 len(micros) > 0 + 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 = {} helper_kwargs = {'normalization_mode': normalization_mode, @@ -154,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, @@ -411,14 +451,126 @@ 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 of the depletable + materials. + + The estimate is computed as the ratio of the neutron production 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_\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 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 + ------- + 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] + + # 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 + atoms = self.number[mat, nuc] + if atoms <= 0.0: + continue + for rxn in micro_xs.reactions: + rate = (atoms * (micro_xs[nuc, rxn] * flux).sum() + / volume_b_cm) + if rxn == 'nu-fission': + production += 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) + 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) + 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 0228a63bbaa..61e87e8d61f 100644 --- a/openmc/deplete/microxs.py +++ b/openmc/deplete/microxs.py @@ -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+ @@ -85,7 +86,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, no energy filter is applied to the flux tally. When @@ -427,7 +431,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.16.1 + Added support for 'nu-fission'. **init_kwargs : dict Keyword arguments passed to :func:`openmc.lib.init` @@ -460,10 +471,13 @@ 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. if reactions is None: reactions = chain.reactions - mts = [REACTION_MT[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)) @@ -478,15 +492,20 @@ def from_multigroup_flux( # 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): - microxs_arr[nuc_index, mt_index, 0] = lib_nuc.collapse_rate( - mt, temperature, energies, multigroup_flux - ) + for mt_index, (rxn_name, mt) in enumerate( + zip(reactions, mts)): + if rxn_name == 'nu-fission': + microxs_arr[nuc_index, mt_index, 0] = \ + 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) return cls(microxs_arr, nuclides, reactions) diff --git a/openmc/lib/nuclide.py b/openmc/lib/nuclide.py index ef1287cf34a..d91228aeaa2 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.16.1 + + 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..c3fd509fa48 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -1072,6 +1072,85 @@ 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); + + // 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. 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& grid = grid_[t_idx].energy; + // 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; + + 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); + + 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 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 = (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); + + rate_sum += flux_per_eV * nuxs_avg * (E_high - E_low); + } + + 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 +1294,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 aca83399a08..abd0e52f13a 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.chain import REACTIONS CHAIN_PATH = Path(__file__).parents[1] / "chain_simple.xml" ONE_GROUP_XS = Path(__file__).parents[1] / "micro_xs_simple.csv" @@ -53,3 +55,156 @@ 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_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(): + # 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) + + # 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') + 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') + 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_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) + op = IndependentOperator([_uranium_material()], [1.0], [micro_xs], + CHAIN_PATH) + assert not op._calculate_kinf + + +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']) + op = IndependentOperator([_uranium_material()], [1.0], [micro_xs], + CHAIN_PATH, keff=(1.0, 0.0)) + assert not op._calculate_kinf diff --git a/tests/unit_tests/test_deplete_microxs.py b/tests/unit_tests/test_deplete_microxs.py index c938f7fb06b..8a0815b2f5e 100644 --- a/tests/unit_tests/test_deplete_microxs.py +++ b/tests/unit_tests/test_deplete_microxs.py @@ -116,6 +116,63 @@ 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_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' + upper = [6.25e-1, 5.53e3, 8.21e5, 2.e7] + flux = [1.0, 0., 0., 0.] + + 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]) + + assert nu_bar(0.0) == pytest.approx(nu_bar(1.0e-5), rel=1e-6) + + def test_microxs_zero_flux(): chain_file = Path(__file__).parents[1] / 'chain_simple.xml'