From 178a898b59bb80611bbd6570dd406f9c10d1cc88 Mon Sep 17 00:00:00 2001 From: shimwell Date: Mon, 6 Jul 2026 12:51:31 +0200 Subject: [PATCH] Add (n,n') inelastic isomeric activation to depletion chains Registers (n,n') as a chain transmutation reaction (MT=4) with entries in REACTIONS, DADZ, REACTION_MT, and the C++ REACTION_TYPE_MAP so that existing depletion operators can tally it by name. Chain.from_endf emits (n,n') entries only when isomeric production data exists for the parent, always including the ground-state self-loop so the Bateman loss and gain terms stay exact; parents without data get no self-loop at all. Covers cases like Nb93(n,n')Nb93_m1 and In115(n,n')In115_m1 that are pure isomeric-state changes with no nuclide change. Tests cover the registry entries, atom conservation of self-loop matrices, and an integration test on the real ENDF/B-VIII.1 Nb93 evaluation. See #121. --- openmc/data/data.py | 1 + openmc/data/reaction.py | 1 + openmc/deplete/chain.py | 9 +++ src/reaction.cpp | 1 + tests/unit_tests/test_deplete_nnprime.py | 93 ++++++++++++++++++++++++ 5 files changed, 105 insertions(+) create mode 100644 tests/unit_tests/test_deplete_nnprime.py diff --git a/openmc/data/data.py b/openmc/data/data.py index c22e54e7dc8..bf5bc58d585 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -119,6 +119,7 @@ } DADZ = { + "(n,n')": (0, 0), '(n,2nd)': (-3, -1), '(n,2n)': (-1, 0), '(n,3n)': (-2, 0), diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index 44ced5511f3..24f9fd2b11c 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -78,6 +78,7 @@ REACTION_MT['fission'] = 18 REACTION_MT['absorption'] = 27 REACTION_MT['capture'] = 102 +REACTION_MT["(n,n')"] = 4 FISSION_MTS = (18, 19, 20, 21, 38) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index dd60a2b7e77..1e80b7066e6 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -56,6 +56,10 @@ '(n,3np)': ReactionInfo({42}, ('H1',)), '(n,n2p)': ReactionInfo({44}, ('H1', 'H1')), '(n,npa)': ReactionInfo({45}, ('H1', 'He4')), + # Inelastic scattering only changes the isomeric state of the target, + # so (n,n') entries are added by Chain.from_endf only when isomeric + # production data exists for the parent + "(n,n')": ReactionInfo({4}, ()), '(n,gamma)': ReactionInfo({102}, ()), '(n,p)': ReactionInfo(set(chain([103], range(600, 650))), ('H1',)), '(n,d)': ReactionInfo(set(chain([104], range(650, 700))), ('H2',)), @@ -586,6 +590,11 @@ def from_endf(cls, decay_files, fpy_files, neutron_files, + ' not used (only the first MT with ' 'data is read)') + # Self-loop inelastic entries are only meaningful + # when they split the parent into isomeric states + if name == "(n,n')" and not records: + continue + if records: cls._add_isomeric_reactions( nuclide, name, daughter, q_value, records, diff --git a/src/reaction.cpp b/src/reaction.cpp index c02e0cc407e..3f6164a6437 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -357,6 +357,7 @@ void initialize_maps() // Alternate names REACTION_TYPE_MAP["elastic"] = ELASTIC; + REACTION_TYPE_MAP["(n,n')"] = N_LEVEL; REACTION_TYPE_MAP["n2n"] = N_2N; REACTION_TYPE_MAP["n3n"] = N_3N; REACTION_TYPE_MAP["n4n"] = N_4N; diff --git a/tests/unit_tests/test_deplete_nnprime.py b/tests/unit_tests/test_deplete_nnprime.py new file mode 100644 index 00000000000..f83adfee8f9 --- /dev/null +++ b/tests/unit_tests/test_deplete_nnprime.py @@ -0,0 +1,93 @@ +"""Tests for (n,n') inelastic isomeric activation in depletion chains.""" + +import os +from pathlib import Path + +import numpy as np +import pytest + +import openmc.data +from openmc.deplete import Chain, reaction_rates +from openmc.deplete.chain import REACTIONS + +NNPRIME_CHAIN = """ + + + + + + + + +""" + + +def test_reaction_registry(): + assert REACTIONS["(n,n')"].mts == {4} + assert REACTIONS["(n,n')"].secondaries == () + assert openmc.data.DADZ["(n,n')"] == (0, 0) + assert openmc.data.REACTION_MT["(n,n')"] == 4 + + +def test_self_loop_matrix_conserves_atoms(tmp_path): + chain_file = tmp_path / "chain.xml" + chain_file.write_text(NNPRIME_CHAIN) + chain = Chain.from_xml(chain_file) + assert chain.reactions == ["(n,n')"] + + rates = reaction_rates.ReactionRates(["1"], ["A", "A_m1"], ["(n,n')"]) + rates.set("1", "A", "(n,n')", 2.0) + rates.set("1", "A_m1", "(n,n')", 0.5) + matrix = chain.form_matrix(rates[0]).toarray() + + # A loses at the full rate but 90% returns as a self-loop; A_m1 fully + # de-excites back to ground + assert matrix[0, 0] == pytest.approx(-2.0 + 2.0 * 0.9) + assert matrix[1, 0] == pytest.approx(2.0 * 0.1) + assert matrix[0, 1] == pytest.approx(0.5) + assert matrix[1, 1] == pytest.approx(-0.5) + + # Atom conservation: every column sums to zero + assert np.allclose(matrix.sum(axis=0), 0.0) + + +def _find_endf(directory, patterns): + for pattern in patterns: + matches = sorted(Path(directory).glob(pattern)) + if matches: + return matches[0] + pytest.skip(f'No file matching {patterns} under {directory}') + + +@pytest.mark.skipif( + 'OPENMC_ENDF_DATA' not in os.environ, + reason='OPENMC_ENDF_DATA environment variable must be set') +def test_from_endf_nnprime(): + endf_data = Path(os.environ['OPENMC_ENDF_DATA']) + decay_files = sorted((endf_data / 'decay').glob('*.endf')) + fpy = _find_endf(endf_data / 'nfy', ['*U*235*']) + nb93 = _find_endf(endf_data / 'neutrons', ['*Nb*93*']) + am241 = _find_endf(endf_data / 'neutrons', ['*Am*241*']) + + chain = Chain.from_endf( + decay_files, [fpy], [nb93, am241], + reactions=["(n,n')"], + progress=False, + isomeric_branching=True) + + # Nb93 has MF=10 data for MT=4: a ground self-loop plus the m1 target + nb = chain['Nb93'] + targets = {rx.target: rx.branching_ratio for rx in nb.reactions + if rx.type == "(n,n')"} + assert targets == {'Nb93': 1.0, 'Nb93_m1': 0.0} + production = nb.isomeric_production[("(n,n')", 'Nb93_m1')] + assert production[0].level == 1 + table = production[0].tables[0] + assert table.mf == 10 + assert len(table.data.x) == 34 + assert table.data(14.356e6) == pytest.approx(3.756e-2) + + # Am241 has no MF=8/9/10 for MT=4, so no self-loop entries are added + assert all(rx.type != "(n,n')" for rx in chain['Am241'].reactions) + + assert chain.validate(strict=True)