Add (n,n') inelastic isomeric activation to depletion chains

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

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)
 

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,

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;

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 = """<depletion_chain>
+  <nuclide name="A" reactions="2">
+    <reaction type="(n,n')" Q="0.0" target="A" branching_ratio="0.9"/>
+    <reaction type="(n,n')" Q="0.0" target="A_m1" branching_ratio="0.1"/>
+  </nuclide>
+  <nuclide name="A_m1" reactions="1">
+    <reaction type="(n,n')" Q="0.0" target="A" branching_ratio="1.0"/>
+  </nuclide>
+</depletion_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)