From abd377a2982ca5ac69382667006bff51c5e58797 Mon Sep 17 00:00:00 2001 From: shimwell Date: Mon, 6 Jul 2026 11:33:12 +0200 Subject: [PATCH 1/4] Store energy-dependent isomeric production data in depletion chain files Adds ProductionTable and IsomericProduction dataclasses to openmc.deplete.nuclide, a Nuclide.isomeric_production mapping keyed by (reaction type, target), XML round-trip via a new child element of , validation of the stored data, format documentation, and unit tests. Also fixes a latent NameError in Nuclide.validate where the reaction branch ratio message referenced the decay branch ratio sum. See #121. --- docs/source/io_formats/depletion_chain.rst | 78 +++++ openmc/deplete/nuclide.py | 289 +++++++++++++++++- .../test_deplete_isomeric_production.py | 258 ++++++++++++++++ 3 files changed, 623 insertions(+), 2 deletions(-) create mode 100644 tests/unit_tests/test_deplete_isomeric_production.py diff --git a/docs/source/io_formats/depletion_chain.rst b/docs/source/io_formats/depletion_chain.rst index 74413e7b610..eb2b02c179c 100644 --- a/docs/source/io_formats/depletion_chain.rst +++ b/docs/source/io_formats/depletion_chain.rst @@ -96,6 +96,84 @@ element has the following attributes: :branching_ratio: The branching ratio for the reaction +In addition to these attributes, a ```` element may contain zero or +more :ref:`io_chain_isomeric` sub-elements carrying verbatim energy-dependent +isomeric production data from the source evaluation. + +.. _io_chain_isomeric: + +--------------------------------- +```` Element +--------------------------------- + +.. versionadded:: 0.15.4 + +The ```` element stores energy-dependent isomeric +production data (ENDF MF=9 yields and/or MF=10 partial cross sections) for the +final state produced by the parent ```` element, exactly as given in +the source evaluation. The scalar ``branching_ratio`` attributes on +```` elements remain authoritative for consumers that do not use +this data, and the scalar branching ratios of a reaction type always sum to +about one. Note that the stored functions are verbatim evaluation data: an +MF=10 partial cross section requires an external total cross section to form a +branching ratio, and the per-state MF=9 yields of a reaction do not +necessarily sum to one when the ground-state share is implicit. + +A ```` element may hold several ```` elements, +e.g. when a final state that could not be matched to a known metastable state +is folded onto the ground-state target. Consumers should sum the tables of all +instances attached to one target; OpenMC never performs arithmetic on the +stored data. + +This element has the following attributes: + + :level: + ENDF LFS level number of the final state, verbatim from the source + evaluation. This is a level index, not a metastable-state index. + + :excitation_energy: + Excitation energy of the final state in [eV] (the MF=8 ELFS value when + present, otherwise QM minus QI). Zero for the ground state. + +Each ```` element contains one or more ```` +sub-elements, each holding one tabulated function with the following +attributes and sub-elements: + + :mf: + ENDF file number the data comes from: 9 for energy-dependent yields + (dimensionless multiplicities of the reaction cross section) or 10 for + partial production cross sections in [b] + + :mt: + ENDF reaction number of the section that supplied the data + + :source: + Source library identifier, e.g. 'ENDF/B-VIII.1' or 'TENDL-2025' + + :QM: + Mass-difference Q value in [eV], verbatim from the TAB1 header + + :QI: + Reaction Q value for this particular state in [eV], verbatim from the + TAB1 header + + :breakpoints: + Whitespace-separated breakpoints of the interpolation regions + + :interpolation: + Whitespace-separated ENDF interpolation scheme codes for each region, + following the same convention as :class:`openmc.data.Tabulated1D` + + :energies: + Sub-element listing the incident neutron energies in [eV] + + :values: + Sub-element listing the tabulated yields or cross sections + +Only ENDF "scheme 1" isomeric data (MF=8/9/10) is represented in the chain +file. Evaluations that encode isomer production via MF=6 product distributions +or via discrete-level reaction sections are not captured. + .. _io_chain_nfy: ------------------------------------ diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index 95881483474..1d386ccebab 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -6,6 +6,7 @@ import bisect from collections.abc import Mapping from collections import namedtuple, defaultdict +from dataclasses import dataclass, field from warnings import warn from numbers import Real @@ -13,12 +14,13 @@ import numpy as np from openmc.checkvalue import check_type +from openmc.data import Tabulated1D from openmc.stats import Univariate from .._xml import get_elem_list, get_text __all__ = [ "DecayTuple", "ReactionTuple", "Nuclide", "FissionYield", - "FissionYieldDistribution"] + "FissionYieldDistribution", "ProductionTable", "IsomericProduction"] DecayTuple = namedtuple('DecayTuple', 'type target branching_ratio') @@ -74,6 +76,194 @@ pass +@dataclass +class ProductionTable: + """Energy-dependent isomeric production data from a single ENDF section. + + Stores one MF=9 yield function or MF=10 partial cross section for the + production of a single final state, exactly as given in the source + evaluation (original energy grid, values, and interpolation regions). + The stored data is never collapsed, renormalized, or re-interpolated. + + .. versionadded:: 0.15.4 + + Parameters + ---------- + mf : int + ENDF file number the data comes from: 9 for energy-dependent yields + (dimensionless multiplicities of the reaction cross section) or 10 + for partial production cross sections in [b]. + mt : int + ENDF reaction number of the section that supplied the data. + source : str + Source library identifier, e.g. 'ENDF/B-VIII.1' or 'TENDL-2025'. + QM : float + Mass-difference Q value in [eV], verbatim from the TAB1 header. + QI : float + Reaction Q value for this particular state in [eV], verbatim from + the TAB1 header. + data : openmc.data.Tabulated1D + Tabulated function of incident neutron energy in [eV]. + + """ + + mf: int + mt: int + source: str + QM: float + QI: float + data: Tabulated1D + + def __repr__(self): + return (f"") + + def __eq__(self, other): + if not isinstance(other, ProductionTable): + return NotImplemented + return ( + self.mf == other.mf + and self.mt == other.mt + and self.source == other.source + and self.QM == other.QM + and self.QI == other.QI + and np.array_equal(self.data.x, other.data.x) + and np.array_equal(self.data.y, other.data.y) + and np.array_equal(self.data.breakpoints, other.data.breakpoints) + and np.array_equal(self.data.interpolation, + other.data.interpolation) + ) + + def to_xml_element(self): + """Write production table to an XML element. + + Returns + ------- + elem : lxml.etree._Element + XML element containing the table data + + """ + elem = ET.Element('table') + elem.set('mf', str(self.mf)) + elem.set('mt', str(self.mt)) + elem.set('source', self.source) + elem.set('QM', str(self.QM)) + elem.set('QI', str(self.QI)) + elem.set('breakpoints', + ' '.join(str(b) for b in self.data.breakpoints)) + elem.set('interpolation', + ' '.join(str(i) for i in self.data.interpolation)) + energy_elem = ET.SubElement(elem, 'energies') + energy_elem.text = ' '.join(str(x) for x in self.data.x) + values_elem = ET.SubElement(elem, 'values') + values_elem.text = ' '.join(str(y) for y in self.data.y) + return elem + + @classmethod + def from_xml_element(cls, element): + """Read production table from an XML element. + + Parameters + ---------- + element : lxml.etree._Element + XML element to read table data from + + Returns + ------- + ProductionTable + + """ + x = get_elem_list(element, 'energies', float) + y = get_elem_list(element, 'values', float) + breakpoints = [int(b) for b in get_text(element, 'breakpoints').split()] + interpolation = [ + int(i) for i in get_text(element, 'interpolation').split()] + return cls( + mf=int(get_text(element, 'mf')), + mt=int(get_text(element, 'mt')), + source=get_text(element, 'source'), + QM=float(get_text(element, 'QM')), + QI=float(get_text(element, 'QI')), + data=Tabulated1D(x, y, breakpoints, interpolation), + ) + + +@dataclass +class IsomericProduction: + """Production data for one evaluation final state routed to one target. + + A single transmutation reaction target may carry more than one instance, + e.g. when a level that could not be matched to a known metastable state + is folded onto the ground-state target. Consumers should sum the tables + of all instances attached to a target; the stored data itself is always + verbatim evaluation data. + + .. versionadded:: 0.15.4 + + Parameters + ---------- + level : int + ENDF LFS level number of the final state, verbatim from the source + evaluation. Note that this is a level index, not a metastable-state + index; the mapping to a metastable target is made by matching + excitation energies against decay data. + excitation_energy : float + Excitation energy of the final state in [eV], taken from the MF=8 + ELFS value when present and otherwise from QM minus QI. Zero for + the ground state. + tables : list of ProductionTable + Production data for this level. May hold both an MF=9 and an MF=10 + table when the evaluation provides both. + + """ + + level: int + excitation_energy: float + tables: list = field(default_factory=list) + + def __repr__(self): + return (f"") + + def to_xml_element(self): + """Write isomeric production data to an XML element. + + Returns + ------- + elem : lxml.etree._Element + XML element containing the production data + + """ + elem = ET.Element('isomeric_production') + elem.set('level', str(self.level)) + elem.set('excitation_energy', str(self.excitation_energy)) + for table in self.tables: + elem.append(table.to_xml_element()) + return elem + + @classmethod + def from_xml_element(cls, element): + """Read isomeric production data from an XML element. + + Parameters + ---------- + element : lxml.etree._Element + XML element to read production data from + + Returns + ------- + IsomericProduction + + """ + return cls( + level=int(get_text(element, 'level')), + excitation_energy=float(get_text(element, 'excitation_energy')), + tables=[ProductionTable.from_xml_element(e) + for e in element.findall('table')], + ) + + class Nuclide: """Decay modes, reactions, and fission yields for a single nuclide. @@ -108,6 +298,13 @@ class Nuclide: treated as a nested dictionary ``{energy: {product: yield}}`` yield_energies : tuple of float or None Energies at which fission product yields exist + isomeric_production : dict + Dictionary mapping ``(reaction type, target)`` tuples to lists of + :class:`IsomericProduction` instances carrying verbatim + energy-dependent isomeric production data from the source + evaluations. + + .. versionadded:: 0.15.4 """ def __init__(self, name=None): @@ -122,6 +319,10 @@ def __init__(self, name=None): # Reaction paths self.reactions = [] + # Energy-dependent isomeric production data, keyed by + # (reaction type, target) + self.isomeric_production = {} + # Decay sources self.sources = {} @@ -270,6 +471,11 @@ def from_xml(cls, element, root=None, fission_q=None): nuc.reactions.append(ReactionTuple( r_type, target, Q, branching_ratio)) + # Check for energy-dependent isomeric production data + for iso_elem in reaction_elem.findall('isomeric_production'): + nuc.isomeric_production.setdefault((r_type, target), []).append( + IsomericProduction.from_xml_element(iso_elem)) + fpy_elem = element.find('neutron_fission_yields') if fpy_elem is not None: # Check for use of FPY from other nuclide @@ -330,6 +536,8 @@ def to_xml_element(self): rx_elem.set('target', daughter) if br != 1.0: rx_elem.set('branching_ratio', str(br)) + for iso in self.isomeric_production.get((rx, daughter), []): + rx_elem.append(iso.to_xml_element()) if self.yield_data: fpy_elem = ET.SubElement(elem, 'neutron_fission_yields') @@ -353,6 +561,9 @@ def validate(self, strict=True, quiet=False, tolerance=1e-4): does the sum of branching ratios equal about one? 2) for fission reactions, does the sum of fission yield fractions equal about two? + 3) does every isomeric production entry correspond to a + reaction present on this nuclide, with yields and cross + sections in physical ranges? Parameters ---------- @@ -414,7 +625,7 @@ def validate(self, strict=True, quiet=False, tolerance=1e-4): if stat: continue msg = msg_func( - name=self.name, actual=sum_br, expected=1.0, tol=tolerance, + name=self.name, actual=sum_rxn, expected=1.0, tol=tolerance, prop=f"{rxn_type} reaction branch ratios") if strict: raise ValueError(msg) @@ -423,6 +634,14 @@ def validate(self, strict=True, quiet=False, tolerance=1e-4): warn(msg) valid = False + for msg in self._check_isomeric_production(tolerance): + if strict: + raise ValueError(msg) + elif quiet: + return False + warn(msg) + valid = False + if self.yield_data: for energy, fission_yield in self.yield_data.items(): sum_yield = fission_yield.yields.sum() @@ -442,6 +661,72 @@ def validate(self, strict=True, quiet=False, tolerance=1e-4): return valid + def _check_isomeric_production(self, tolerance): + """Generate messages describing inconsistent isomeric production data + + Parameters + ---------- + tolerance : float + Absolute tolerance for comparisons + + Yields + ------ + str + Description of each inconsistency found + + """ + if not self.isomeric_production: + return + + rx_pairs = {(rx.type, rx.target) for rx in self.reactions} + for (rxn_type, target), productions in self.isomeric_production.items(): + if (rxn_type, target) not in rx_pairs: + yield (f"Nuclide {self.name} has isomeric production data " + f"for {rxn_type} to {target}, which is not a reaction " + "present on this nuclide") + for prod in productions: + if prod.level < 0 or prod.excitation_energy < 0.0: + yield (f"Nuclide {self.name} has an isomeric production " + f"entry for {rxn_type} to {target} with a negative " + "level or excitation energy") + for table in prod.tables: + y = table.data.y + if table.mf == 9 and ( + y.min() < 0.0 or y.max() > 1.0 + tolerance): + yield (f"Nuclide {self.name} has MF=9 isomeric " + f"yields for {rxn_type} to {target} outside " + "of [0, 1]") + elif table.mf == 10 and y.min() < 0.0: + yield (f"Nuclide {self.name} has negative MF=10 " + f"isomeric production cross sections for " + f"{rxn_type} to {target}") + + # Pointwise sum of MF=9 yields across all states of a reaction type. + # Only checked when every state has MF=9 data and a ground-state + # entry exists; otherwise the ground share is implicit and no sum + # rule holds. Evaluated on the intersection of the tabulated ranges + # so that differing thresholds do not produce clamp artifacts. + by_type = defaultdict(list) + for (rxn_type, _), productions in self.isomeric_production.items(): + by_type[rxn_type].extend(productions) + for rxn_type, productions in by_type.items(): + tables = [next((t for t in prod.tables if t.mf == 9), None) + for prod in productions] + if None in tables or all(p.level != 0 for p in productions): + continue + low = max(t.data.x[0] for t in tables) + high = min(t.data.x[-1] for t in tables) + union = np.unique(np.concatenate([t.data.x for t in tables])) + union = union[(union >= low) & (union <= high)] + if union.size == 0: + continue + total = sum(t.data(union) for t in tables) + if total.max() > 1.0 + tolerance: + yield ("Nuclide {} has MF=9 isomeric yields for {} that " + "sum to {:7.4e} at some energies instead of at most " + "1 +/- {:7.4e}").format( + self.name, rxn_type, total.max(), tolerance) + class FissionYieldDistribution(Mapping): """Energy-dependent fission product yields for a single nuclide diff --git a/tests/unit_tests/test_deplete_isomeric_production.py b/tests/unit_tests/test_deplete_isomeric_production.py new file mode 100644 index 00000000000..1d08075556a --- /dev/null +++ b/tests/unit_tests/test_deplete_isomeric_production.py @@ -0,0 +1,258 @@ +"""Tests for energy-dependent isomeric production data in depletion chains.""" + +import lxml.etree as ET +import numpy as np +import pytest + +from openmc.data import Tabulated1D +from openmc.deplete import Chain, IsomericProduction, Nuclide, ProductionTable +from openmc.deplete.nuclide import ReactionTuple + +# Values from the ENDF/B-VIII.1 Am241 evaluation, MF=9, MT=102 +AM241_ENERGIES = ("1e-05 0.369 1000.0 100000.0 600000.1 1000000.0 " + "2000000.0 4000000.1 30000000.0") +AM241_GROUND = "0.9 0.9 0.8667 0.842 0.81533 0.74382 0.5703 0.52 0.52" +AM241_M1 = "0.1 0.1 0.1333 0.158 0.18467 0.25618 0.4297 0.48 0.48" + +CHAIN_XML = f""" + + + + +
+ {AM241_ENERGIES} + {AM241_GROUND} +
+ + + + + + {AM241_ENERGIES} + {AM241_M1} +
+ + + + + + + + 1e-05 1000.0 100000.0 20000000.0 + 0.5 0.55 0.6 0.7 +
+ + 1e-05 1000.0 100000.0 20000000.0 + 0.5 0.02 0.001 1e-05 +
+ + + + 1500000.0 5000000.0 20000000.0 + 0.0 0.001 1e-05 +
+ + + + + + 1e-05 1000.0 100000.0 20000000.0 + 0.5 0.45 0.4 0.3 +
+ + + + + + + + +""" + + +@pytest.fixture +def chain(tmp_path): + chain_file = tmp_path / "chain.xml" + chain_file.write_text(CHAIN_XML) + return Chain.from_xml(chain_file) + + +def test_from_xml(chain): + am241 = chain["Am241"] + + # Reaction tuples are unchanged by the new child elements + assert am241.reactions == [ + ReactionTuple("(n,gamma)", "Am242", 5537755.0, 1.0), + ReactionTuple("(n,gamma)", "Am242_m1", 5537755.0, 0.0), + ] + + # Ground and metastable entries each carry their own verbatim table + ground = am241.isomeric_production[("(n,gamma)", "Am242")] + meta = am241.isomeric_production[("(n,gamma)", "Am242_m1")] + assert len(ground) == 1 and len(meta) == 1 + assert ground[0].level == 0 + assert ground[0].excitation_energy == 0.0 + assert meta[0].level == 2 + assert meta[0].excitation_energy == 48630.0 + + table = meta[0].tables[0] + assert table.mf == 9 + assert table.mt == 102 + assert table.source == "ENDF/B-VIII.1" + assert table.QM == 5537755.0 + assert table.QI == 5489125.0 + np.testing.assert_array_equal(table.data.breakpoints, [9]) + np.testing.assert_array_equal(table.data.interpolation, [3]) + np.testing.assert_array_equal( + table.data.x, [float(x) for x in AM241_ENERGIES.split()]) + np.testing.assert_array_equal( + table.data.y, [float(y) for y in AM241_M1.split()]) + + # The tables are callable functions of energy + assert table.data(0.0253) == pytest.approx(0.1) + assert ground[0].tables[0].data(0.0253) == pytest.approx(0.9) + + +def test_multiple_tables_and_folded_levels(chain): + nb93 = chain["Nb93"] + ground = nb93.isomeric_production[("(n,gamma)", "Nb94")] + + # Ground target holds its own level plus a folded unmappable level + assert [p.level for p in ground] == [0, 7] + + # MF=9 and MF=10 tables coexist on one level + assert [t.mf for t in ground[0].tables] == [9, 10] + + # Multi-region interpolation survives parsing + table = ground[0].tables[0] + np.testing.assert_array_equal(table.data.breakpoints, [2, 4]) + np.testing.assert_array_equal(table.data.interpolation, [2, 5]) + + # The folded level keeps its own excitation energy and threshold grid + folded = ground[1] + assert folded.excitation_energy == 1300000.0 + assert folded.tables[0].data.x[0] == 1500000.0 + + +def test_xml_roundtrip(chain, tmp_path): + out = tmp_path / "chain_out.xml" + chain.export_to_xml(out) + reread = Chain.from_xml(out) + + for nuclide in chain.nuclides: + other = reread[nuclide.name] + assert other.reactions == nuclide.reactions + assert other.isomeric_production == nuclide.isomeric_production + + # A second write produces identical bytes (bit-exact round-trip) + out2 = tmp_path / "chain_out2.xml" + reread.export_to_xml(out2) + assert out.read_text() == out2.read_text() + + +def test_legacy_reader_compatibility(chain, tmp_path): + # A reader that only looks at attributes sees a normal chain + out = tmp_path / "chain_out.xml" + chain.export_to_xml(out) + root = ET.parse(str(out)) + am241 = root.find('nuclide[@name="Am241"]') + reactions = am241.findall("reaction") + assert len(reactions) == 2 + assert reactions[0].get("target") == "Am242" + assert reactions[0].get("branching_ratio") is None + assert reactions[1].get("target") == "Am242_m1" + assert reactions[1].get("branching_ratio") == "0.0" + + # Scalar branching ratios per reaction type still sum to one + assert chain.validate(strict=True) + + +def test_validate_orphan_entry(): + nuc = Nuclide("A") + nuc.add_reaction("(n,gamma)", "B", 0.0, 1.0) + table = ProductionTable( + mf=9, mt=102, source="test", QM=0.0, QI=0.0, + data=Tabulated1D([1.0, 2.0], [0.5, 0.5])) + nuc.isomeric_production[("(n,gamma)", "B_m1")] = [ + IsomericProduction(1, 100.0, [table])] + with pytest.raises(ValueError, match="not a reaction present"): + nuc.validate(strict=True) + with pytest.warns(UserWarning, match="not a reaction present"): + assert not nuc.validate(strict=False) + + +def test_validate_mf9_range(): + nuc = Nuclide("A") + nuc.add_reaction("(n,gamma)", "B", 0.0, 1.0) + table = ProductionTable( + mf=9, mt=102, source="test", QM=0.0, QI=0.0, + data=Tabulated1D([1.0, 2.0], [0.5, 1.2])) + nuc.isomeric_production[("(n,gamma)", "B")] = [ + IsomericProduction(0, 0.0, [table])] + with pytest.raises(ValueError, match="outside of"): + nuc.validate(strict=True) + + +def test_validate_mf10_negative(): + nuc = Nuclide("A") + nuc.add_reaction("(n,gamma)", "B", 0.0, 1.0) + table = ProductionTable( + mf=10, mt=102, source="test", QM=0.0, QI=0.0, + data=Tabulated1D([1.0, 2.0], [0.1, -0.1])) + nuc.isomeric_production[("(n,gamma)", "B")] = [ + IsomericProduction(0, 0.0, [table])] + with pytest.raises(ValueError, match="negative MF=10"): + nuc.validate(strict=True) + + +def test_validate_mf9_sum(): + # Ground and metastable MF=9 yields that sum above one must fail + nuc = Nuclide("A") + nuc.add_reaction("(n,gamma)", "B", 0.0, 1.0) + nuc.add_reaction("(n,gamma)", "B_m1", 0.0, 0.0) + ground = ProductionTable( + mf=9, mt=102, source="test", QM=0.0, QI=0.0, + data=Tabulated1D([1.0, 2.0], [0.9, 0.9])) + meta = ProductionTable( + mf=9, mt=102, source="test", QM=0.0, QI=-100.0, + data=Tabulated1D([1.0, 2.0], [0.15, 0.1])) + nuc.isomeric_production[("(n,gamma)", "B")] = [ + IsomericProduction(0, 0.0, [ground])] + nuc.isomeric_production[("(n,gamma)", "B_m1")] = [ + IsomericProduction(1, 100.0, [meta])] + with pytest.raises(ValueError, match="sum to"): + nuc.validate(strict=True) + + +def test_validate_implicit_ground_share(): + # A metastable-only MF=9 yield (like In115 capture, Y=0.79 to In116_m1) + # has an implicit ground share, so no sum rule applies + nuc = Nuclide("In115") + nuc.add_reaction("(n,gamma)", "In116", 0.0, 1.0) + nuc.add_reaction("(n,gamma)", "In116_m1", 0.0, 0.0) + table = ProductionTable( + mf=9, mt=102, source="test", QM=6784719.0, QI=6657319.0, + data=Tabulated1D([1e-5, 3e7], [0.79, 0.79])) + nuc.isomeric_production[("(n,gamma)", "In116_m1")] = [ + IsomericProduction(2, 127400.0, [table])] + assert nuc.validate(strict=True) + + +def test_validate_quiet(): + nuc = Nuclide("A") + nuc.add_reaction("(n,gamma)", "B", 0.0, 1.0) + table = ProductionTable( + mf=9, mt=102, source="test", QM=0.0, QI=0.0, + data=Tabulated1D([1.0, 2.0], [0.5, 1.2])) + nuc.isomeric_production[("(n,gamma)", "B")] = [ + IsomericProduction(0, 0.0, [table])] + assert not nuc.validate(strict=False, quiet=True) + + +def test_no_isomeric_data(): + nuc = Nuclide("A") + nuc.add_reaction("(n,gamma)", "B", 0.0, 1.0) + assert nuc.isomeric_production == {} + assert nuc.validate(strict=True) + elem = nuc.to_xml_element() + assert elem.find("reaction").find("isomeric_production") is None From 9ade545a060994e5a7fddf71b1b64c4cb7e67ab4 Mon Sep 17 00:00:00 2001 From: shimwell Date: Mon, 6 Jul 2026 11:36:34 +0200 Subject: [PATCH 2/4] Carry isomeric production data through chain reduce and branch ratio edits Chain.reduce deep-copies isomeric production data for retained targets and warns when data is dropped with an excluded target. Chain.set_branch_ratios keeps attached data for targets that survive a rewrite and, controlled by a new preserve_isomeric_data argument, raises or warns when a rewrite would remove a target carrying data. Adds a Chain.get_isomeric_production accessor and tests. See #121. --- openmc/deplete/chain.py | 107 +++++++++++++++++- .../test_deplete_isomeric_production.py | 66 +++++++++++ 2 files changed, 169 insertions(+), 4 deletions(-) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index a2ff6dea0b6..24e91e53123 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -4,6 +4,7 @@ loaded from an .xml file and all the nuclides are linked together. """ +from copy import deepcopy from io import StringIO from itertools import chain import math @@ -963,8 +964,37 @@ def get_branch_ratios(self, reaction="(n,gamma)"): capt[nuclide.name] = nuc_capt return capt + def get_isomeric_production(self, nuclide, reaction): + """Return energy-dependent isomeric production data for one reaction + + .. versionadded:: 0.15.4 + + Parameters + ---------- + nuclide : str + Name of the parent nuclide, e.g. ``"Am241"`` + reaction : str + Reaction name, e.g. ``"(n,gamma)"`` + + Returns + ------- + dict + Mapping of target nuclide names to lists of + :class:`~openmc.deplete.IsomericProduction` instances. Empty if + the nuclide carries no data for this reaction. + + """ + parent = self[nuclide] + return { + target: productions + for (rx_type, target), productions + in parent.isomeric_production.items() + if rx_type == reaction + } + def set_branch_ratios(self, branch_ratios, reaction="(n,gamma)", - strict=True, tolerance=1e-5): + strict=True, tolerance=1e-5, + preserve_isomeric_data=True): """Set the branching ratios for a given reactions Parameters @@ -988,6 +1018,15 @@ def set_branch_ratios(self, branch_ratios, reaction="(n,gamma)", 1 - tol < sum_br < 1 + tol + preserve_isomeric_data : bool, optional + Energy-dependent isomeric production data attached to targets + that remain present is always kept. If a rewrite would remove a + target that carries such data, an error is raised when this + evaluates to ``True`` [default]; otherwise a warning is issued + and the data is dropped. + + .. versionadded:: 0.15.4 + Raises ------ IndexError @@ -1001,7 +1040,10 @@ def set_branch_ratios(self, branch_ratios, reaction="(n,gamma)", ``branch_ratios`` does not have the requested reaction ValueError If ``strict`` evalutes to ``False`` and the sum of one parents - branch ratios is outside 1 +/- ``tolerance`` + branch ratios is outside 1 +/- ``tolerance``, or if + ``preserve_isomeric_data`` evaluates to ``True`` and a target + carrying energy-dependent isomeric production data would be + removed See Also -------- @@ -1106,6 +1148,38 @@ def set_branch_ratios(self, branch_ratios, reaction="(n,gamma)", "with a sum outside tolerance of 1 +/- {:5.3e}:\n{}".format( reaction, tolerance, "\n".join(tail))) + # Check up front whether the rewrite would remove targets that carry + # energy-dependent isomeric production data, before any mutation + + discarded = [] + for parent_name, rxn_index in rxn_ix_map.items(): + parent = self[parent_name] + new_ratios = branch_ratios[parent_name] + kept = set(new_ratios) + # Account for the ground target that will be added automatically + if all("_m" in t for t in new_ratios) and sums[parent_name] != 1.0: + ground_target = grounds.get(parent_name) + if ground_target is None: + pz, pa, pm = zam(parent_name) + ground_target = gnds_name(pz, pa + 1, 0) + kept.add(ground_target) + for ix in rxn_index: + target = parent.reactions[ix].target + if (target not in kept + and (reaction, target) in parent.isomeric_production): + discarded.append((parent_name, target)) + + if discarded: + tail = ", ".join(f"{p} -> {t}" for p, t in discarded) + msg = (f"Setting {reaction} branch ratios would remove targets " + "that carry energy-dependent isomeric production data: " + f"{tail}") + if preserve_isomeric_data: + raise ValueError( + msg + ". Pass preserve_isomeric_data=False to drop the " + "data.") + warn(msg) + # Insert new ReactionTuples with updated branch ratios for parent_name, rxn_index in rxn_ix_map.items(): @@ -1117,9 +1191,15 @@ def set_branch_ratios(self, branch_ratios, reaction="(n,gamma)", # Assume Q value is independent of target state rxn_Q = parent.reactions[rxn_index[0]].Q - # Remove existing reactions + # Remove existing reactions, saving attached isomeric + # production data + saved_production = {} for ix in reversed(rxn_index): - parent.reactions.pop(ix) + popped = parent.reactions.pop(ix) + data = parent.isomeric_production.pop( + (reaction, popped.target), None) + if data is not None: + saved_production[popped.target] = data # Add new reactions all_meta = True @@ -1138,6 +1218,11 @@ def set_branch_ratios(self, branch_ratios, reaction="(n,gamma)", new_ratios[ground_target] = ground_br parent.add_reaction(reaction, ground_target, rxn_Q, ground_br) + # Reattach isomeric production data for surviving targets + for target, data in saved_production.items(): + if target in new_ratios: + parent.isomeric_production[(reaction, target)] = data + @property def fission_yields(self): if self._fission_yields is None: @@ -1262,6 +1347,7 @@ def reduce(self, initial_isotopes, level=None): name_sort = sorted(all_isotopes) new_chain = type(self)() + dropped_production = [] for idx, iso in enumerate(sorted(all_isotopes, key=openmc.data.zam)): previous = self[iso] @@ -1281,6 +1367,10 @@ def reduce(self, initial_isotopes, level=None): for rx in previous.reactions: if rx.target in all_isotopes: new_nuclide.add_reaction(*rx) + key = (rx.type, rx.target) + if key in previous.isomeric_production: + new_nuclide.isomeric_production[key] = deepcopy( + previous.isomeric_production[key]) elif rx.type == "fission": new_yields = new_nuclide.yield_data = ( previous.yield_data.restrict_products(name_sort)) @@ -1289,9 +1379,18 @@ def reduce(self, initial_isotopes, level=None): # Maintain total destruction rates but set no target else: new_nuclide.add_reaction(rx.type, None, rx.Q, rx.branching_ratio) + if (rx.type, rx.target) in previous.isomeric_production: + dropped_production.append( + (previous.name, rx.type, rx.target)) new_chain.add_nuclide(new_nuclide) + if dropped_production: + tail = ", ".join(f"{parent} {rx_type} -> {target}" + for parent, rx_type, target in dropped_production) + warn("Energy-dependent isomeric production data was dropped " + "because the target is not in the reduced chain: " + tail) + # Doesn't appear that the ordering matters for the reactions, # just the contents new_chain.reactions = sorted(new_chain.reactions) diff --git a/tests/unit_tests/test_deplete_isomeric_production.py b/tests/unit_tests/test_deplete_isomeric_production.py index 1d08075556a..5c5e17b2517 100644 --- a/tests/unit_tests/test_deplete_isomeric_production.py +++ b/tests/unit_tests/test_deplete_isomeric_production.py @@ -256,3 +256,69 @@ def test_no_isomeric_data(): assert nuc.validate(strict=True) elem = nuc.to_xml_element() assert elem.find("reaction").find("isomeric_production") is None + + +def test_get_isomeric_production(chain): + data = chain.get_isomeric_production("Am241", "(n,gamma)") + assert set(data) == {"Am242", "Am242_m1"} + assert data["Am242_m1"][0].level == 2 + assert chain.get_isomeric_production("Am242", "(n,gamma)") == {} + + +def test_reduce_carries_data(chain): + reduced = chain.reduce(["Am241"]) + original = chain["Am241"].isomeric_production + carried = reduced["Am241"].isomeric_production + assert carried == original + # Deep copy, not shared references + key = ("(n,gamma)", "Am242_m1") + assert carried[key][0] is not original[key][0] + + +def test_reduce_drops_data_with_warning(chain): + with pytest.warns(UserWarning, match="dropped"): + reduced = chain.reduce(["Nb93"], level=0) + nb93 = reduced["Nb93"] + assert nb93.isomeric_production == {} + # Total destruction rate entries are still present, with no target + assert all(rx.target is None for rx in nb93.reactions) + + +def test_set_branch_ratios_preserves_data(chain): + before = dict(chain["Am241"].isomeric_production) + chain.set_branch_ratios({"Am241": {"Am242": 0.89, "Am242_m1": 0.11}}) + am241 = chain["Am241"] + ratios = {rx.target: rx.branching_ratio for rx in am241.reactions + if rx.type == "(n,gamma)"} + assert ratios == {"Am242": 0.89, "Am242_m1": 0.11} + assert am241.isomeric_production == before + + +def test_set_branch_ratios_infers_ground_and_preserves(chain): + before = dict(chain["Am241"].isomeric_production) + chain.set_branch_ratios({"Am241": {"Am242_m1": 0.1}}) + am241 = chain["Am241"] + ratios = {rx.target: rx.branching_ratio for rx in am241.reactions + if rx.type == "(n,gamma)"} + assert ratios == {"Am242": pytest.approx(0.9), "Am242_m1": 0.1} + assert am241.isomeric_production == before + + +def test_set_branch_ratios_discard_protection(chain): + # Removing a target that carries data raises by default + with pytest.raises(ValueError, match="Am241 -> Am242_m1"): + chain.set_branch_ratios({"Am241": {"Am242": 1.0}}) + + # Nothing was mutated by the failed call + assert ("(n,gamma)", "Am242_m1") in chain["Am241"].isomeric_production + + # Opting out drops the data with a warning + with pytest.warns(UserWarning, match="Am241 -> Am242_m1"): + chain.set_branch_ratios({"Am241": {"Am242": 1.0}}, + preserve_isomeric_data=False) + am241 = chain["Am241"] + assert ("(n,gamma)", "Am242_m1") not in am241.isomeric_production + targets = [rx.target for rx in am241.reactions if rx.type == "(n,gamma)"] + assert targets == ["Am242"] + # The surviving target keeps its data + assert ("(n,gamma)", "Am242") in am241.isomeric_production From 670515292b977826e0d103781d4abc7321929f5b Mon Sep 17 00:00:00 2001 From: shimwell Date: Mon, 6 Jul 2026 11:53:37 +0200 Subject: [PATCH 3/4] Read MF=8/9/10 isomeric production data in Chain.from_endf Adds an openmc.deplete._isomeric module that parses ENDF MF=8/9/10 radioactive production sections from raw section text, keeping the MF=9 yields and MF=10 partial cross sections verbatim, and maps final state level numbers to metastable targets by matching excitation energies against decay data (with positional pairing as fallback and fold-to-ground for unmappable levels). Chain.from_endf gains isomeric_branching, branching_files (extra evaluations consulted where the primary neutron file has no data, e.g. TENDL), scalar_branching ('none' default keeps historical behavior; 'thermal' and multigroup flux collapse available), elis_rtol, and isomer_mapping_log arguments. Metastable targets are emitted as additional reaction entries and every mapping decision can be written to a report file. Tested against verbatim ENDF/B-VIII.1 Am241 and Nb93 section fixtures and an integration test on real decay plus neutron files reproducing the Am241 thermal capture split of 0.90 to ground and 0.10 to Am242_m1 from level LFS=2. See #121. --- openmc/deplete/_isomeric.py | 430 ++++++++++++++++++ openmc/deplete/chain.py | 180 +++++++- .../test_deplete_isomeric_extraction.py | 305 +++++++++++++ 3 files changed, 908 insertions(+), 7 deletions(-) create mode 100644 openmc/deplete/_isomeric.py create mode 100644 tests/unit_tests/test_deplete_isomeric_extraction.py diff --git a/openmc/deplete/_isomeric.py b/openmc/deplete/_isomeric.py new file mode 100644 index 00000000000..daf4814f9af --- /dev/null +++ b/openmc/deplete/_isomeric.py @@ -0,0 +1,430 @@ +"""Support for extracting energy-dependent isomeric production data. + +This module reads ENDF MF=8/9/10 radioactive nuclide production sections +from incident neutron evaluations and maps the final-state level numbers +(LFS) found there to metastable nuclide names by matching excitation +energies against decay data. The tabulated MF=9 yields and MF=10 partial +cross sections are kept verbatim; no arithmetic is performed on them. + +Note that an LFS value is a level index of the product nuclide, not a +metastable-state index, so targets are never named directly from it. +""" + +from collections import Counter, defaultdict +from dataclasses import dataclass +from io import StringIO + +import numpy as np + +from openmc.data import gnds_name, zam +import openmc.data.endf as endf6 + +from .nuclide import IsomericProduction, ProductionTable + +# Assignment status values used in the mapping report +MATCHED = 'matched' +POSITIONAL = 'positional' +GROUND = 'ground' +FOLDED = 'folded' + +# Final states with an excitation energy below this value [eV] are treated +# as the ground state regardless of their level number +GROUND_ENERGY_CUTOFF = 1.0 + +THERMAL_ENERGY = 0.0253 + + +@dataclass +class LevelRecord: + """Raw MF=9 or MF=10 data for the production of one final state.""" + + mf: int + mt: int + zap: int # ZA of the product; may be 0 when not given + lfs: int # ENDF level number of the final state + elfs: float # MF=8 excitation energy [eV]; None when absent + qm: float # [eV] + qi: float # [eV] + data: object # openmc.data.Tabulated1D + source: str # library label, e.g. 'ENDF/B-8.1' + + @property + def excitation_energy(self): + """Excitation energy [eV], preferring the MF=8 ELFS value.""" + if self.elfs is not None: + return self.elfs + return self.qm - self.qi + + +@dataclass +class LevelAssignment: + """The chain target chosen for one level record.""" + + record: LevelRecord + target: str + status: str + liso: int = None + decay_elis: float = None + + +def library_label(evaluation): + """Return a short provenance label for an evaluation. + + Parameters + ---------- + evaluation : openmc.data.endf.Evaluation + Evaluation to label + + Returns + ------- + str + Label such as 'ENDF/B-8.1' or 'TENDL-2025.0' + + """ + library, version, release = evaluation.info['library'] + return f'{library}-{version}.{release}' + + +def extract_isomeric_production(evaluation, source=None): + """Read MF=8/9/10 production sections from an evaluation. + + MF=8 subsections that point at other files for the production data + (LMF not equal to 9 or 10, e.g. the LMF=6 product distributions found + on TENDL MT=5) are skipped. + + Parameters + ---------- + evaluation : openmc.data.endf.Evaluation + Incident neutron evaluation to read + source : str, optional + Provenance label stored on each record. Defaults to + :func:`library_label` of the evaluation. + + Returns + ------- + dict + Mapping of MT numbers to lists of :class:`LevelRecord` + + """ + if source is None: + source = library_label(evaluation) + + by_mt = defaultdict(set) + for mf, mt in evaluation.section: + if mf in (9, 10): + by_mt[mt].add(mf) + + result = {} + for mt in sorted(by_mt): + # MF=8 links each (ZAP, LFS) pair to an excitation energy and says + # in which file the production data lives (LMF) + mf8_info = {} + if (8, mt) in evaluation.section: + file_obj = StringIO(evaluation.section[8, mt]) + items = endf6.get_head_record(file_obj) + n_states, complete_flag = items[4], items[5] + for _ in range(n_states): + if complete_flag == 0: + sub, _values = endf6.get_list_record(file_obj) + else: + sub = endf6.get_cont_record(file_obj) + zap, elfs = int(sub[0]), float(sub[1]) + lmf, lfs = int(sub[2]), int(sub[3]) + mf8_info[zap, lfs] = (elfs, lmf) + + levels = [] + for mf in sorted(by_mt[mt]): + file_obj = StringIO(evaluation.section[mf, mt]) + items = endf6.get_head_record(file_obj) + n_states = items[4] + for _ in range(n_states): + params, func = endf6.get_tab1_record(file_obj) + qm, qi = params[0], params[1] + zap, lfs = int(params[2]), int(params[3]) + elfs, lmf = mf8_info.get((zap, lfs), (None, None)) + if lmf is not None and lmf not in (9, 10): + continue + levels.append(LevelRecord( + mf=mf, mt=mt, zap=zap, lfs=lfs, elfs=elfs, + qm=qm, qi=qi, data=func, source=source)) + if levels: + result[mt] = levels + return result + + +def assign_levels(records, ground_daughter, isomer_energies, elis_rtol): + """Assign the level records of one reaction to chain target names. + + Levels are first matched to metastable states by excitation energy + (nearest decay-library state within ``elis_rtol`` relative tolerance). + Levels that fail the energy match are paired positionally with the + remaining metastables in order of increasing excitation energy, and + anything left over is folded onto the ground-state target. + + Parameters + ---------- + records : list of LevelRecord + Level records for a single transmutation reaction + ground_daughter : str + Ground-state product from mass and charge arithmetic, used when a + record does not carry a product ZA + isomer_energies : dict + Mapping ``{(Z, A): [(liso, elis), ...]}`` of metastable states with + positive excitation energies from the decay data + elis_rtol : float + Relative tolerance on the excitation energy match + + Returns + ------- + list of LevelAssignment + + """ + assignments = [None] * len(records) + unmatched = [] + matched_liso = defaultdict(set) + + for i, record in enumerate(records): + if record.zap > 0: + z, a = divmod(record.zap, 1000) + else: + z, a, _ = zam(ground_daughter) + target_elis = record.excitation_energy + if record.lfs == 0 or target_elis < GROUND_ENERGY_CUTOFF: + assignments[i] = LevelAssignment( + record, gnds_name(z, a, 0), GROUND) + continue + + best = None + for liso, elis in isomer_energies.get((z, a), []): + diff = abs(target_elis - elis) + if diff <= elis_rtol * elis and (best is None or diff < best[2]): + best = (liso, elis, diff) + if best is not None: + liso, elis, _diff = best + assignments[i] = LevelAssignment( + record, gnds_name(z, a, liso), MATCHED, liso, elis) + matched_liso[z, a].add(liso) + else: + unmatched.append((i, record, z, a, target_elis)) + + # Positional fallback per product nuclide + by_za = defaultdict(list) + for entry in unmatched: + by_za[entry[2], entry[3]].append(entry) + for (z, a), entries in by_za.items(): + available = [ + (liso, elis) + for liso, elis in sorted(isomer_energies.get((z, a), [])) + if liso not in matched_liso[z, a]] + entries.sort(key=lambda entry: entry[4]) + for (i, record, *_), (liso, elis) in zip(entries, available): + assignments[i] = LevelAssignment( + record, gnds_name(z, a, liso), POSITIONAL, liso, elis) + for i, record, *_ in entries[len(available):]: + assignments[i] = LevelAssignment( + record, gnds_name(z, a, 0), FOLDED) + + return assignments + + +def group_records(records_and_targets): + """Group level records into IsomericProduction objects per target. + + Records with the same level number attached to the same target (an + MF=9 and an MF=10 table for one state) share one + :class:`~openmc.deplete.IsomericProduction` instance. + + Parameters + ---------- + records_and_targets : list of (LevelRecord, str) + Level records with the final chain target each is attached to + + Returns + ------- + dict + Mapping of target names to lists of IsomericProduction + + """ + result = defaultdict(list) + index = {} + for record, target in records_and_targets: + key = (target, record.lfs) + production = index.get(key) + if production is None: + elfs = record.excitation_energy + if record.lfs == 0: + elfs = max(elfs, 0.0) + production = IsomericProduction(record.lfs, elfs, []) + index[key] = production + result[target].append(production) + production.tables.append(ProductionTable( + mf=record.mf, mt=record.mt, source=record.source, + QM=record.qm, QI=record.qi, data=record.data)) + return dict(result) + + +def _collapse(table, mode): + """Collapse an MF=9 yield table to a scalar for the requested mode.""" + x = table.data.x + if mode == 'thermal': + if x[0] > THERMAL_ENERGY or x[-1] < THERMAL_ENERGY: + return None + return float(table.data(THERMAL_ENERGY)) + + energies, flux = mode + energies = np.asarray(energies, dtype=float) + flux = np.asarray(flux, dtype=float) + total = flux.sum() + if total <= 0.0: + return None + midpoints = np.sqrt(energies[:-1] * energies[1:]) + inside = (midpoints >= x[0]) & (midpoints <= x[-1]) + values = np.zeros_like(midpoints) + values[inside] = table.data(midpoints[inside]) + return float(values @ flux / total) + + +def compute_scalar_ratios(productions_by_target, ground_target, mode): + """Compute scalar branching ratios per target for one reaction. + + The default mode ``'none'`` reproduces the historical behavior: the + full reaction rate goes to the ground-state target. The ``'thermal'`` + mode evaluates the MF=9 yields at 0.0253 eV, and an ``(energies, + flux)`` tuple collapses them with a multigroup flux (evaluated at the + geometric midpoints of the group boundaries). The ground-state share + is always taken as one minus the metastable sum, which also covers + evaluations where the ground yield is implicit. When the requested + mode cannot be supported by the stored data (an MF=10-only state, a + grid not covering the spectrum, or metastable yields summing above + one), the ratios fall back to the ``'none'`` values and a note is + returned. + + Parameters + ---------- + productions_by_target : dict + Mapping of target names to lists of IsomericProduction + ground_target : str + Name of the ground-state target + mode : {'none', 'thermal'} or tuple + Scalar branching mode + + Returns + ------- + ratios : dict + Mapping of target names to scalar branching ratios summing to one + note : str or None + Explanation when the ratios fell back to the default + + """ + default = {target: 0.0 for target in productions_by_target} + default[ground_target] = 1.0 + if mode == 'none' or set(default) == {ground_target}: + return default, None + + ratios = dict(default) + meta_sum = 0.0 + for target, productions in productions_by_target.items(): + if target == ground_target: + continue + value = 0.0 + for production in productions: + table = next( + (t for t in production.tables if t.mf == 9), None) + if table is None: + return default, ( + f'no MF=9 yield for {target}; scalar branching ratios ' + 'left at the default') + collapsed = _collapse(table, mode) + if collapsed is None: + return default, ( + f'MF=9 grid for {target} does not cover the requested ' + 'spectrum; scalar branching ratios left at the default') + value += collapsed + ratios[target] = value + meta_sum += value + + if meta_sum > 1.0: + return default, ( + 'metastable MF=9 yields sum above one; scalar branching ' + 'ratios left at the default') + ratios[ground_target] = 1.0 - meta_sum + return ratios, None + + +class IsomerMappingReport: + """Collects level-to-isomer mapping decisions made during a chain build. + + Parameters + ---------- + elis_rtol : float + Relative tolerance that was used for excitation energy matching + + """ + + def __init__(self, elis_rtol): + self.elis_rtol = elis_rtol + self.rows = [] + self.notes = [] + + def add(self, parent, reaction, assignment, target, status, note=''): + record = assignment.record + self.rows.append({ + 'parent': parent, + 'reaction': reaction, + 'mt': record.mt, + 'mf': record.mf, + 'lfs': record.lfs, + 'elfs': record.excitation_energy, + 'target': target, + 'status': status, + 'liso': assignment.liso, + 'decay_elis': assignment.decay_elis, + 'source': record.source, + 'note': note, + }) + + def add_note(self, parent, reaction, text): + self.notes.append((parent, reaction, text)) + + @property + def counts(self): + return Counter(row['status'] for row in self.rows) + + def write(self, path): + """Write a human readable mapping report. + + Parameters + ---------- + path : str or os.PathLike + File to write the report to + + """ + columns = ('parent', 'reaction', 'mt', 'mf', 'lfs', 'elfs', + 'target', 'status', 'liso', 'decay_elis', 'source', + 'note') + header = ('ISOMER MAPPING REPORT\n' + f'ELIS matching relative tolerance: {self.elis_rtol}\n') + counts = self.counts + summary = ['Level records: {}'.format(sum(counts.values()))] + for status in (MATCHED, POSITIONAL, GROUND, FOLDED): + summary.append(f' {status}: {counts.get(status, 0)}') + + lines = [header, '\n'.join(summary), ''] + lines.append(' '.join(f'{c:>12}' for c in columns)) + for row in self.rows: + formatted = [] + for column in columns: + value = row[column] + if value is None: + value = '-' + elif isinstance(value, float): + value = f'{value:.6g}' + formatted.append(f'{value!s:>12}') + lines.append(' '.join(formatted)) + if self.notes: + lines.append('') + lines.append('NOTES') + for parent, reaction, text in self.notes: + lines.append(f' {parent} {reaction}: {text}') + with open(path, 'w') as fh: + fh.write('\n'.join(lines) + '\n') diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 24e91e53123..8601b87a7ca 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -23,6 +23,9 @@ from openmc.data import gnds_name, zam from openmc.exceptions import DataError from .nuclide import FissionYieldDistribution, Nuclide +from ._isomeric import ( + MATCHED, POSITIONAL, GROUND, FOLDED, IsomerMappingReport, assign_levels, + compute_scalar_ratios, extract_isomeric_production, group_records) from .._xml import get_text from .._sparse_compat import csc_array, dok_array import openmc.data @@ -314,14 +317,16 @@ def add_nuclide(self, nuclide: Nuclide): @classmethod def from_endf(cls, decay_files, fpy_files, neutron_files, reactions=('(n,2n)', '(n,3n)', '(n,4n)', '(n,gamma)', '(n,p)', '(n,a)'), - progress=True + progress=True, isomeric_branching=False, branching_files=None, + scalar_branching='none', elis_rtol=0.5, isomer_mapping_log=None ): """Create a depletion chain from ENDF files. - String arguments in ``decay_files``, ``fpy_files``, and - ``neutron_files`` will be treated as file names to be read. - Alternatively, :class:`openmc.data.endf.Evaluation` or - ``endf.Material`` instances can be included in these arguments. + String arguments in ``decay_files``, ``fpy_files``, + ``neutron_files``, and ``branching_files`` will be treated as file + names to be read. Alternatively, + :class:`openmc.data.endf.Evaluation` or ``endf.Material`` instances + can be included in these arguments. Parameters ---------- @@ -341,6 +346,43 @@ def from_endf(cls, decay_files, fpy_files, neutron_files, progress : bool, optional Flag to print status messages during processing. Does not effect warning messages + isomeric_branching : bool, optional + Read energy-dependent isomeric production data (ENDF MF=8/9/10) + from the neutron files, add metastable reaction targets, and + store the MF=9 yields and MF=10 partial cross sections verbatim + on the chain. + + .. versionadded:: 0.15.4 + branching_files : list of str, openmc.data.endf.Evaluation, or endf.Material, optional + Additional ENDF evaluations consulted for MF=8/9/10 data for + parent nuclides whose primary neutron file has none, e.g. TENDL + files. Consulted in list order; only their MF=8/9/10 sections + are read, so they cannot change which reactions exist. Each + stored table is tagged with the library it came from. + + .. versionadded:: 0.15.4 + scalar_branching : {'none', 'thermal'} or tuple, optional + How to compute the scalar ``branching_ratio`` attributes for + reactions with isomeric targets. ``'none'`` [default] keeps the + historical behavior (ground state gets 1.0, metastables 0.0). + ``'thermal'`` evaluates the MF=9 yields at 0.0253 eV. A tuple + of ``(energies, flux)`` collapses the MF=9 yields with a + multigroup flux given on group boundaries ``energies`` in [eV]. + Reactions whose data cannot support the requested mode fall + back to ``'none'`` with a note in the mapping report. + + .. versionadded:: 0.15.4 + elis_rtol : float, optional + Relative tolerance used when matching the excitation energy of + a produced level against the excitation energies of metastable + states in the decay data. + + .. versionadded:: 0.15.4 + isomer_mapping_log : str, optional + Path to write a report of every level-to-isomer mapping + decision. + + .. versionadded:: 0.15.4 Returns ------- @@ -359,10 +401,18 @@ def from_endf(cls, decay_files, fpy_files, neutron_files, """ transmutation_reactions = reactions + if (scalar_branching not in ('none', 'thermal') + and not isinstance(scalar_branching, tuple)): + raise ValueError( + "scalar_branching must be 'none', 'thermal', or a tuple " + "of (energies, flux)") + report = IsomerMappingReport(elis_rtol) if isomeric_branching else None + # Create dictionary mapping target to filename if progress: print('Processing neutron sub-library files...') reactions = {} + iso_records = {} for f in neutron_files: evaluation = openmc.data.endf.as_evaluation(f) name = evaluation.gnds_name @@ -373,17 +423,48 @@ def from_endf(cls, decay_files, fpy_files, neutron_files, openmc.data.endf.get_head_record(file_obj) q_value = openmc.data.endf.get_cont_record(file_obj)[1] reactions[name][mt] = q_value + if isomeric_branching: + levels = extract_isomeric_production(evaluation) + if levels: + iso_records[name] = levels + + # Read supplementary evaluations, consulted only for parents whose + # primary neutron file has no MF=8/9/10 data. The first file in the + # list providing a section wins. + iso_supplement = {} + if isomeric_branching and branching_files: + if progress: + print('Processing supplementary isomeric branching files...') + for f in branching_files: + evaluation = openmc.data.endf.as_evaluation(f) + levels = extract_isomeric_production(evaluation) + if not levels: + continue + store = iso_supplement.setdefault(evaluation.gnds_name, {}) + for mt, records in levels.items(): + store.setdefault(mt, records) # Determine what decay and FPY nuclides are available if progress: print('Processing decay sub-library files...') decay_data = {} + isomer_energies = defaultdict(list) for f in decay_files: - data = openmc.data.Decay(f) + ev = openmc.data.endf.as_evaluation(f) + data = openmc.data.Decay(ev) # Skip decay data for neutron itself if data.nuclide['atomic_number'] == 0: continue decay_data[data.nuclide['name']] = data + # Record metastable excitation energies for level matching. + # States reported with a zero excitation energy cannot be + # matched and are skipped. + liso = data.nuclide['isomeric_state'] + elis = ev.target['excitation_energy'] + if liso > 0 and elis > 0.0: + key = (data.nuclide['atomic_number'], + data.nuclide['mass_number']) + isomer_energies[key].append((liso, elis)) if progress: print('Processing fission product yield sub-library files...') @@ -463,7 +544,24 @@ def from_endf(cls, decay_files, fpy_files, neutron_files, else: q_value = 0.0 - nuclide.add_reaction(name, daughter, q_value, 1.0) + # Gather energy-dependent isomeric production data, + # with the primary neutron file taking precedence + # over supplementary files + records = [] + if isomeric_branching and daughter is not None: + available = iso_records.get(parent, {}) + if not (mts & set(available)): + available = iso_supplement.get(parent, {}) + for mt in sorted(mts & set(available)): + records.extend(available[mt]) + + if records: + cls._add_isomeric_reactions( + nuclide, name, daughter, q_value, records, + decay_data, isomer_energies, elis_rtol, + scalar_branching, report) + else: + nuclide.add_reaction(name, daughter, q_value, 1.0) if any(mt in reactions_available for mt in openmc.data.FISSION_MTS): q_value = reactions[parent][18] @@ -533,8 +631,76 @@ def from_endf(cls, decay_files, fpy_files, neutron_files, for vals in missing_fp: print(' {}, E={} eV (total yield={})'.format(*vals)) + if report is not None: + counts = report.counts + if progress and counts: + print('Isomeric level mapping: ' + ', '.join( + f'{status}={counts.get(status, 0)}' + for status in (MATCHED, POSITIONAL, GROUND, FOLDED))) + if isomer_mapping_log is not None: + report.write(isomer_mapping_log) + return chain + @staticmethod + def _add_isomeric_reactions(nuclide, name, daughter, q_value, records, + decay_data, isomer_energies, elis_rtol, + scalar_branching, report): + """Add ground and metastable entries for one transmutation reaction. + + Level records are assigned to targets by excitation energy + matching, grouped into verbatim production data per target, and + emitted as one ReactionTuple per target with scalar branching + ratios according to ``scalar_branching``. + """ + parent = nuclide.name + assignments = assign_levels( + records, daughter, isomer_energies, elis_rtol) + + # Resolve assignments to nuclides that exist in the decay data, + # folding anything that cannot be resolved onto the ground target + resolved = [] + target_levels = defaultdict(set) + for assignment in assignments: + target, status = assignment.target, assignment.status + if status in (MATCHED, POSITIONAL): + if target not in decay_data: + replacement = replace_missing(target, decay_data) + if replacement is None or '_m' not in replacement: + target, status = daughter, FOLDED + else: + target = replacement + target_levels[target].add(assignment.record.lfs) + else: + target = daughter + resolved.append((assignment, target, status)) + + if report is not None: + for assignment, target, status in resolved: + note = '' + if (status in (MATCHED, POSITIONAL) + and len(target_levels[target]) > 1): + note = 'shares target with another level' + report.add(parent, name, assignment, target, status, note) + + grouped = group_records( + [(assignment.record, target) + for assignment, target, _status in resolved]) + ratios, note = compute_scalar_ratios( + grouped, daughter, scalar_branching) + if note and report is not None: + report.add_note(parent, name, note) + + nuclide.add_reaction(name, daughter, q_value, + ratios.get(daughter, 1.0)) + if daughter in grouped: + nuclide.isomeric_production[(name, daughter)] = grouped[daughter] + for target in sorted(grouped): + if target == daughter: + continue + nuclide.add_reaction(name, target, q_value, ratios[target]) + nuclide.isomeric_production[(name, target)] = grouped[target] + @classmethod def from_xml(cls, filename, fission_q=None): """Reads a depletion chain XML file. diff --git a/tests/unit_tests/test_deplete_isomeric_extraction.py b/tests/unit_tests/test_deplete_isomeric_extraction.py new file mode 100644 index 00000000000..c1cb08cacc9 --- /dev/null +++ b/tests/unit_tests/test_deplete_isomeric_extraction.py @@ -0,0 +1,305 @@ +"""Tests for ENDF MF=8/9/10 extraction and level-to-isomer mapping. + +The fixture strings below are verbatim section text from the +ENDF/B-VIII.1 Am241 and Nb93 incident neutron evaluations. +""" + +import os +from pathlib import Path + +import numpy as np +import pytest + +import openmc.deplete +from openmc.deplete import Chain +from openmc.deplete._isomeric import ( + FOLDED, GROUND, MATCHED, POSITIONAL, assign_levels, + compute_scalar_ratios, extract_isomeric_production, group_records) + +AM241_MF8_MT102 = """\ + 9.524100+4 2.389860+2 0 0 2 09543 8102 1 + 9.524201+4 0.000000+0 9 0 24 09543 8102 2 + 5.767200+4 1.400000+0 9.624201+4 3.686630-1 6.648200+5 2.010000+09543 8102 3 + 5.767200+4 1.400000+0 9.624201+4 4.583370-1 6.226900+5 2.020000+09543 8102 4 + 5.767200+4 2.400000+0 9.424201+4 6.700000-2 7.509630+5 2.010000+09543 8102 5 + 5.767200+4 2.400000+0 9.424201+4 1.060000-1 7.064230+5 2.020000+09543 8102 6 + 9.524201+4 4.863000+4 9 2 120 09543 8102 7 + 4.446580+9 3.100000+0 9.524201+4 8.232040-1 4.863000+4 2.000000+09543 8102 8 + 4.446580+9 3.200000+0 9.524201+4 1.722060-1 4.863000+4 2.010000+09543 8102 9 + 4.446580+9 4.100000+0 9.323801+4 2.747740-7 5.561940+6 2.020000+09543 8102 10 + 4.446580+9 4.100000+0 9.323801+4 6.411380-6 5.501660+6 2.040000+09543 8102 11 + 4.446580+9 4.100000+0 9.323801+4 4.762740-5 5.452290+6 2.070000+09543 8102 12 + 4.446580+9 4.100000+0 9.323801+4 5.358080-5 5.409190+6 2.090000+09543 8102 13 + 4.446580+9 4.100000+0 9.323801+4 3.159900-5 5.355540+6 2.110000+09543 8102 14 + 4.446580+9 4.100000+0 9.323801+4 3.938420-5 5.312840+6 2.130000+09543 8102 15 + 4.446580+9 4.100000+0 9.323801+4 1.831820-6 5.290040+6 2.140000+09543 8102 16 + 4.446580+9 4.100000+0 9.323801+4 5.037510-6 5.287640+6 2.150000+09543 8102 17 + 4.446580+9 4.100000+0 9.323801+4 1.373870-6 5.254040+6 2.160000+09543 8102 18 + 4.446580+9 4.100000+0 9.323801+4 4.114280-3 5.245940+6 2.170000+09543 8102 19 + 4.446580+9 4.100000+0 9.323801+4 9.159120-7 5.191740+6 2.200000+09543 8102 20 + 4.446580+9 4.100000+0 9.323801+4 2.665300-4 5.179840+6 2.210000+09543 8102 21 + 4.446580+9 4.100000+0 9.323801+4 8.701160-6 5.125740+6 2.220000+09543 8102 22 + 4.446580+9 4.100000+0 9.323801+4 1.373870-6 5.119740+6 2.230000+09543 8102 23 + 4.446580+9 4.100000+0 9.323801+4 1.007500-5 5.101140+6 2.240000+09543 8102 24 + 4.446580+9 4.100000+0 9.323801+4 9.159120-7 5.063440+6 2.250000+09543 8102 25 + 4.446580+9 4.100000+0 9.323801+4 9.159120-8 5.010340+6 2.260000+09543 8102 26 + 4.446580+9 6.000000+0 0.000000+0 4.70000-11 0.000000+0 1.000000+09543 8102 27 +""" + +AM241_MF9_MT102 = """\ + 9.524100+4 2.389860+2 0 0 2 09543 9102 1 + 5.537755+6 5.537755+6 95242 0 1 99543 9102 2 + 9 3 9543 9102 3 + 1.000000-5 9.000000-1 3.690000-1 9.000000-1 1.000000+3 8.667000-19543 9102 4 + 1.000000+5 8.420000-1 6.000001+5 8.153300-1 1.000000+6 7.438200-19543 9102 5 + 2.000000+6 5.703000-1 4.000001+6 5.200000-1 3.000000+7 5.200000-19543 9102 6 + 5.537755+6 5.489125+6 95242 2 1 99543 9102 7 + 9 3 9543 9102 8 + 1.000000-5 1.000000-1 3.690000-1 1.000000-1 1.000000+3 1.333000-19543 9102 9 + 1.000000+5 1.580000-1 6.000001+5 1.846700-1 1.000000+6 2.561800-19543 9102 10 + 2.000000+6 4.297000-1 4.000001+6 4.800000-1 3.000000+7 4.800000-19543 9102 11 +""" + +NB93_MF8_MT4 = """\ + 4.109300+4 9.210510+1 0 0 1 14125 8 4 1 + 4.109300+4 3.073000+4 10 1 0 04125 8 4 2 +""" + +NB93_MF10_MT4 = """\ + 4.109300+4 9.210510+1 0 0 1 0412510 4 1 + 0.000000+0-3.073000+4 41093 1 1 34412510 4 2 + 34 2 0 0 0 0412510 4 3 + 3.115000+4 0.000000+0 4.000000+4 6.053000-5 5.000000+4 1.513000-4412510 4 4 + 6.000000+4 2.522000-4 8.000000+4 5.246000-4 1.000000+5 8.676000-4412510 4 5 + 1.500000+5 1.987000-3 2.000000+5 3.460000-3 3.000000+5 7.253000-3412510 4 6 + 4.000000+5 1.198000-2 5.000000+5 1.742000-2 6.000000+5 2.337000-2412510 4 7 + 7.000000+5 2.921000-2 7.625000+5 4.035000-2 1.000000+6 7.484000-2412510 4 8 + 1.288000+6 1.122000-1 1.750000+6 2.110000-1 2.063000+6 2.060000-1412510 4 9 + 2.563000+6 2.486000-1 4.050000+6 2.782000-1 4.588000+6 2.629000-1412510 4 10 + 6.156000+6 2.455000-1 8.344000+6 2.446000-1 9.375000+6 2.089000-1412510 4 11 + 1.037500+7 1.349000-1 1.139400+7 8.441000-2 1.239400+7 5.764000-2412510 4 12 + 1.340600+7 4.385000-2 1.435600+7 3.756000-2 1.552000+7 3.088000-2412510 4 13 + 1.700000+7 2.718000-2 2.000000+7 2.265000-2 2.000000+7 0.000000+0412510 4 14 + 1.500000+8 0.000000+0 412510 4 15 +""" + + +class StubEvaluation: + """Minimal object exposing the raw section text of an evaluation.""" + + def __init__(self, sections): + self.section = sections + + +@pytest.fixture +def am241_records(): + stub = StubEvaluation({ + (8, 102): AM241_MF8_MT102, + (9, 102): AM241_MF9_MT102, + }) + return extract_isomeric_production(stub, source='ENDF/B-8.1') + + +def test_extract_am241(am241_records): + assert set(am241_records) == {102} + ground, meta = am241_records[102] + + assert ground.mf == 9 + assert ground.zap == 95242 + assert ground.lfs == 0 + assert ground.elfs == 0.0 + assert ground.qm == 5537755.0 + assert ground.qi == 5537755.0 + assert ground.source == 'ENDF/B-8.1' + + assert meta.lfs == 2 + assert meta.elfs == 48630.0 + assert meta.qi == 5489125.0 + assert meta.excitation_energy == 48630.0 + + # Verbatim 9-point single-region lin-log tables + np.testing.assert_array_equal(ground.data.breakpoints, [9]) + np.testing.assert_array_equal(ground.data.interpolation, [3]) + assert len(meta.data.x) == 9 + assert ground.data(0.0253) == pytest.approx(0.9) + assert meta.data(0.0253) == pytest.approx(0.1) + + +def test_extract_nb93(): + stub = StubEvaluation({ + (8, 4): NB93_MF8_MT4, + (10, 4): NB93_MF10_MT4, + }) + records = extract_isomeric_production(stub, source='ENDF/B-8.1') + (record,) = records[4] + assert record.mf == 10 + assert record.zap == 41093 + assert record.lfs == 1 + assert record.elfs == 30730.0 + assert record.qm == 0.0 + assert record.qi == -30730.0 + assert len(record.data.x) == 34 + assert record.data.x[0] == 31150.0 + # Partial cross section is kept in barns, never divided by a total + assert record.data(14.356e6) == pytest.approx(3.756e-2) + + +def test_extract_skips_other_schemes(): + # An MF=8 subsection pointing at LMF=6 (scheme 2, e.g. TENDL MT=5) + # must not produce a record even if an MF=10 section exists + mf8 = NB93_MF8_MT4.replace( + ' 4.109300+4 3.073000+4 10 1', + ' 4.109300+4 3.073000+4 6 1') + stub = StubEvaluation({(8, 4): mf8, (10, 4): NB93_MF10_MT4}) + assert extract_isomeric_production(stub, source='test') == {} + + +def test_assign_matched(am241_records): + # ELIS matching: level 2 at 48630 eV matches the decay metastable at + # 48600 eV, giving Am242_m1 (not the level-numbered Am242_m2) + assignments = assign_levels( + am241_records[102], 'Am242', {(95, 242): [(1, 48600.0)]}, 0.5) + ground, meta = assignments + assert ground.target == 'Am242' + assert ground.status == GROUND + assert meta.target == 'Am242_m1' + assert meta.status == MATCHED + assert meta.liso == 1 + assert meta.decay_elis == 48600.0 + + +def test_assign_positional_fallback(am241_records): + # An excitation energy failing the relative tolerance falls back to + # positional pairing with the unused metastables + assignments = assign_levels( + am241_records[102], 'Am242', {(95, 242): [(1, 10000.0)]}, 0.5) + meta = assignments[1] + assert meta.target == 'Am242_m1' + assert meta.status == POSITIONAL + + +def test_assign_folded(am241_records): + # With no metastables in the decay data the level folds onto ground + assignments = assign_levels(am241_records[102], 'Am242', {}, 0.5) + meta = assignments[1] + assert meta.target == 'Am242' + assert meta.status == FOLDED + + +def test_assign_two_isomers_no_cross_match(): + # Ir192-like case: two levels must map to m1 and m2 without crossing + records = extract_isomeric_production(StubEvaluation({ + (8, 102): AM241_MF8_MT102, + (9, 102): AM241_MF9_MT102, + }), source='test')[102] + # Reuse the Am241 records but pretend there are two isomers, with the + # level at 48630 eV clearly closer to the first + isomers = {(95, 242): [(1, 48600.0), (2, 2200000.0)]} + assignments = assign_levels(records, 'Am242', isomers, 0.5) + assert assignments[1].target == 'Am242_m1' + assert assignments[1].liso == 1 + + +def test_group_records_merges_mf9_mf10(): + stub = StubEvaluation({ + (8, 4): NB93_MF8_MT4, + (10, 4): NB93_MF10_MT4, + }) + (record,) = extract_isomeric_production(stub, source='test')[4] + duplicate = extract_isomeric_production(StubEvaluation({ + (8, 4): NB93_MF8_MT4, + (10, 4): NB93_MF10_MT4, + }), source='test')[4][0] + duplicate.mf = 9 + grouped = group_records([ + (record, 'Nb93_m1'), (duplicate, 'Nb93_m1')]) + (production,) = grouped['Nb93_m1'] + assert production.level == 1 + assert [t.mf for t in production.tables] == [10, 9] + + +def test_scalar_ratios_thermal(am241_records): + assignments = assign_levels( + am241_records[102], 'Am242', {(95, 242): [(1, 48600.0)]}, 0.5) + grouped = group_records( + [(a.record, a.target) for a in assignments]) + ratios, note = compute_scalar_ratios(grouped, 'Am242', 'thermal') + assert note is None + assert ratios == {'Am242': pytest.approx(0.9), + 'Am242_m1': pytest.approx(0.1)} + + +def test_scalar_ratios_flux_collapse(am241_records): + assignments = assign_levels( + am241_records[102], 'Am242', {(95, 242): [(1, 48600.0)]}, 0.5) + grouped = group_records( + [(a.record, a.target) for a in assignments]) + # All flux in a single low-energy group where the yield is constant + mode = ([1.0e-5, 1.0, 3.0e7], [1.0, 0.0]) + ratios, note = compute_scalar_ratios(grouped, 'Am242', mode) + assert note is None + assert ratios['Am242_m1'] == pytest.approx(0.1) + assert ratios['Am242'] == pytest.approx(0.9) + + +def test_scalar_ratios_mf10_fallback(): + stub = StubEvaluation({ + (8, 4): NB93_MF8_MT4, + (10, 4): NB93_MF10_MT4, + }) + (record,) = extract_isomeric_production(stub, source='test')[4] + grouped = group_records([(record, 'Nb93_m1')]) + ratios, note = compute_scalar_ratios(grouped, 'Nb93', 'thermal') + assert 'no MF=9 yield' in note + assert ratios == {'Nb93': 1.0, 'Nb93_m1': 0.0} + + +def test_scalar_ratios_none_mode(am241_records): + assignments = assign_levels( + am241_records[102], 'Am242', {(95, 242): [(1, 48600.0)]}, 0.5) + grouped = group_records( + [(a.record, a.target) for a in assignments]) + ratios, note = compute_scalar_ratios(grouped, 'Am242', 'none') + assert note is None + assert ratios == {'Am242': 1.0, 'Am242_m1': 0.0} + + +def test_from_endf_invalid_scalar_mode(): + with pytest.raises(ValueError, match='scalar_branching'): + Chain.from_endf([], [], [], scalar_branching='banana') + + +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_isomeric(): + endf_data = Path(os.environ['OPENMC_ENDF_DATA']) + neutron = _find_endf(endf_data / 'neutrons', ['*Am*241*']) + decay_files = sorted((endf_data / 'decay').glob('*.endf')) + fpy = _find_endf(endf_data / 'nfy', ['*U*235*']) + + chain = Chain.from_endf( + decay_files, [fpy], [neutron], + reactions=['(n,gamma)'], + progress=False, + isomeric_branching=True, + scalar_branching='thermal') + + am241 = chain['Am241'] + targets = {rx.target: rx.branching_ratio for rx in am241.reactions + if rx.type == '(n,gamma)'} + assert targets == {'Am242': pytest.approx(0.9), + 'Am242_m1': pytest.approx(0.1)} + production = am241.isomeric_production[('(n,gamma)', 'Am242_m1')] + assert production[0].level == 2 + assert production[0].tables[0].mf == 9 + assert production[0].tables[0].data(0.0253) == pytest.approx(0.1) + assert chain.validate(strict=True) From 42238b162f62c41214ac8e7a2ef85427c425e6f8 Mon Sep 17 00:00:00 2001 From: shimwell Date: Mon, 6 Jul 2026 12:45:49 +0200 Subject: [PATCH 4/4] Apply review fixes to isomeric production extraction Fixes from an adversarial review of the branch: - Multigroup scalar collapse now clamps group midpoints into the MF=9 grid instead of silently zeroing groups outside it, which discarded the entire thermal group for group structures starting at 0 eV. - The positional level-to-isomer fallback now pairs whole levels, not individual records, so an MF=9 and MF=10 record of one state can no longer be split across two different metastable targets. - Production data attached to a (type, target) pair appearing on more than one reaction element no longer multiplies on XML round trip. - from_endf validates the scalar_branching tuple shape and flux sum, validates elis_rtol, and raises when branching_files, scalar_branching, or isomer_mapping_log are passed without isomeric_branching=True. Reactions with MF=9/10 data on more than one MT use the first (summary) MT and note the skipped ones in the mapping report. - Deduplicated the ground-target inference in set_branch_ratios, added ProductionTable and IsomericProduction to the Python API docs, fixed source label examples to match evaluation headers, corrected type annotations, and added tests for all of the above plus the mapping report writer, the replace_missing folding path, and the validate message fix. See #121. --- docs/source/io_formats/depletion_chain.rst | 3 +- docs/source/pythonapi/deplete.rst | 2 + openmc/deplete/_isomeric.py | 123 +++++++++----- openmc/deplete/chain.py | 86 +++++++--- openmc/deplete/nuclide.py | 17 +- .../test_deplete_isomeric_extraction.py | 157 +++++++++++++++++- .../test_deplete_isomeric_production.py | 40 +++++ 7 files changed, 354 insertions(+), 74 deletions(-) diff --git a/docs/source/io_formats/depletion_chain.rst b/docs/source/io_formats/depletion_chain.rst index eb2b02c179c..4f354106aad 100644 --- a/docs/source/io_formats/depletion_chain.rst +++ b/docs/source/io_formats/depletion_chain.rst @@ -148,7 +148,8 @@ attributes and sub-elements: ENDF reaction number of the section that supplied the data :source: - Source library identifier, e.g. 'ENDF/B-VIII.1' or 'TENDL-2025' + Source library identifier as recorded in the evaluation header, + e.g. 'ENDF/B-8.1' or 'TENDL-2023.1' :QM: Mass-difference Q value in [eV], verbatim from the TAB1 header diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index 25fcd898f49..1c98559b633 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -124,6 +124,8 @@ for a depletion chain: ReactionTuple FissionYieldDistribution FissionYield + ProductionTable + IsomericProduction The :class:`Chain` class uses information from the following module variable: diff --git a/openmc/deplete/_isomeric.py b/openmc/deplete/_isomeric.py index daf4814f9af..bbc8eddac50 100644 --- a/openmc/deplete/_isomeric.py +++ b/openmc/deplete/_isomeric.py @@ -16,7 +16,7 @@ import numpy as np -from openmc.data import gnds_name, zam +from openmc.data import Tabulated1D, gnds_name, zam import openmc.data.endf as endf6 from .nuclide import IsomericProduction, ProductionTable @@ -40,13 +40,13 @@ class LevelRecord: mf: int mt: int - zap: int # ZA of the product; may be 0 when not given - lfs: int # ENDF level number of the final state - elfs: float # MF=8 excitation energy [eV]; None when absent - qm: float # [eV] - qi: float # [eV] - data: object # openmc.data.Tabulated1D - source: str # library label, e.g. 'ENDF/B-8.1' + zap: int # ZA of the product; may be 0 when not given + lfs: int # ENDF level number of the final state + elfs: float | None # MF=8 excitation energy [eV]; None when absent + qm: float # [eV] + qi: float # [eV] + data: Tabulated1D + source: str # library label, e.g. 'ENDF/B-8.1' @property def excitation_energy(self): @@ -63,8 +63,8 @@ class LevelAssignment: record: LevelRecord target: str status: str - liso: int = None - decay_elis: float = None + liso: int | None = None + decay_elis: float | None = None def library_label(evaluation): @@ -207,22 +207,26 @@ def assign_levels(records, ground_daughter, isomer_energies, elis_rtol): else: unmatched.append((i, record, z, a, target_elis)) - # Positional fallback per product nuclide - by_za = defaultdict(list) - for entry in unmatched: - by_za[entry[2], entry[3]].append(entry) - for (z, a), entries in by_za.items(): + # Positional fallback per product nuclide. Pairing is done per LEVEL, + # not per record, since one level may carry both an MF=9 and an MF=10 + # record and both must land on the same target. + by_za_level = defaultdict(lambda: defaultdict(list)) + for i, record, z, a, target_elis in unmatched: + by_za_level[z, a][record.lfs].append((i, record, target_elis)) + for (z, a), levels in by_za_level.items(): available = [ (liso, elis) for liso, elis in sorted(isomer_energies.get((z, a), [])) if liso not in matched_liso[z, a]] - entries.sort(key=lambda entry: entry[4]) - for (i, record, *_), (liso, elis) in zip(entries, available): - assignments[i] = LevelAssignment( - record, gnds_name(z, a, liso), POSITIONAL, liso, elis) - for i, record, *_ in entries[len(available):]: - assignments[i] = LevelAssignment( - record, gnds_name(z, a, 0), FOLDED) + ordered = sorted(levels.values(), key=lambda entries: entries[0][2]) + for entries, (liso, elis) in zip(ordered, available): + for i, record, _elis in entries: + assignments[i] = LevelAssignment( + record, gnds_name(z, a, liso), POSITIONAL, liso, elis) + for entries in ordered[len(available):]: + for i, record, _elis in entries: + assignments[i] = LevelAssignment( + record, gnds_name(z, a, 0), FOLDED) return assignments @@ -263,8 +267,47 @@ def group_records(records_and_targets): return dict(result) +def select_records(mts, primary, supplement): + """Choose the level records for one transmutation reaction. + + The primary library wins over the supplement. Within a library, only + the first MT with data (in ascending order, so a summary section is + preferred over any partials) is used so that a summary section and + its partials cannot be double counted; any further MTs carrying data + are returned so the caller can flag them. + + Parameters + ---------- + mts : set of int + Candidate ENDF MT numbers of the reaction + primary : dict + Mapping ``{mt: [LevelRecord]}`` from the primary evaluation + supplement : dict + Mapping ``{mt: [LevelRecord]}`` from supplementary evaluations + + Returns + ------- + records : list of LevelRecord + Records of the first MT with data, or an empty list + skipped_mts : list of int + Additional MTs that also carry data but were not used + + """ + available = primary if (mts & set(primary)) else supplement + with_data = sorted(mts & set(available)) + if not with_data: + return [], [] + return available[with_data[0]], with_data[1:] + + def _collapse(table, mode): - """Collapse an MF=9 yield table to a scalar for the requested mode.""" + """Collapse an MF=9 yield table to a scalar for the requested mode. + + Returns None when the thermal mode is requested but the table's grid + does not cover 0.0253 eV. In multigroup mode the group midpoints are + clamped into the tabulated range, so every group contributes using + the table's edge values. + """ x = table.data.x if mode == 'thermal': if x[0] > THERMAL_ENERGY or x[-1] < THERMAL_ENERGY: @@ -274,14 +317,14 @@ def _collapse(table, mode): energies, flux = mode energies = np.asarray(energies, dtype=float) flux = np.asarray(flux, dtype=float) - total = flux.sum() - if total <= 0.0: - return None + # Group representative energies: geometric midpoints of the group + # boundaries, clamped into the tabulated range so that groups + # straddling or outside the grid (e.g. a first group starting at + # 0 eV) use the table's edge values instead of contributing zero midpoints = np.sqrt(energies[:-1] * energies[1:]) - inside = (midpoints >= x[0]) & (midpoints <= x[-1]) - values = np.zeros_like(midpoints) - values[inside] = table.data(midpoints[inside]) - return float(values @ flux / total) + midpoints = np.clip(midpoints, x[0], x[-1]) + values = table.data(midpoints) + return float(values @ flux / flux.sum()) def compute_scalar_ratios(productions_by_target, ground_target, mode): @@ -290,14 +333,14 @@ def compute_scalar_ratios(productions_by_target, ground_target, mode): The default mode ``'none'`` reproduces the historical behavior: the full reaction rate goes to the ground-state target. The ``'thermal'`` mode evaluates the MF=9 yields at 0.0253 eV, and an ``(energies, - flux)`` tuple collapses them with a multigroup flux (evaluated at the - geometric midpoints of the group boundaries). The ground-state share - is always taken as one minus the metastable sum, which also covers - evaluations where the ground yield is implicit. When the requested - mode cannot be supported by the stored data (an MF=10-only state, a - grid not covering the spectrum, or metastable yields summing above - one), the ratios fall back to the ``'none'`` values and a note is - returned. + flux)`` tuple collapses them with a multigroup flux, evaluated at the + geometric midpoints of the group boundaries clamped into the + tabulated range. The ground-state share is always taken as one minus + the metastable sum, which also covers evaluations where the ground + yield is implicit. When the requested mode cannot be supported by + the stored data (an MF=10-only state, a thermal request outside the + grid, or metastable yields summing above one), the ratios fall back + to the ``'none'`` values and a note is returned. Parameters ---------- @@ -337,8 +380,8 @@ def compute_scalar_ratios(productions_by_target, ground_target, mode): collapsed = _collapse(table, mode) if collapsed is None: return default, ( - f'MF=9 grid for {target} does not cover the requested ' - 'spectrum; scalar branching ratios left at the default') + f'MF=9 grid for {target} does not cover 0.0253 eV; ' + 'scalar branching ratios left at the default') value += collapsed ratios[target] = value meta_sum += value diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 8601b87a7ca..dd60a2b7e77 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -25,7 +25,8 @@ from .nuclide import FissionYieldDistribution, Nuclide from ._isomeric import ( MATCHED, POSITIONAL, GROUND, FOLDED, IsomerMappingReport, assign_levels, - compute_scalar_ratios, extract_isomeric_production, group_records) + compute_scalar_ratios, extract_isomeric_production, group_records, + select_records) from .._xml import get_text from .._sparse_compat import csc_array, dok_array import openmc.data @@ -358,7 +359,8 @@ def from_endf(cls, decay_files, fpy_files, neutron_files, parent nuclides whose primary neutron file has none, e.g. TENDL files. Consulted in list order; only their MF=8/9/10 sections are read, so they cannot change which reactions exist. Each - stored table is tagged with the library it came from. + stored table is tagged with the library named in the + evaluation header. Requires ``isomeric_branching=True``. .. versionadded:: 0.15.4 scalar_branching : {'none', 'thermal'} or tuple, optional @@ -367,9 +369,12 @@ def from_endf(cls, decay_files, fpy_files, neutron_files, historical behavior (ground state gets 1.0, metastables 0.0). ``'thermal'`` evaluates the MF=9 yields at 0.0253 eV. A tuple of ``(energies, flux)`` collapses the MF=9 yields with a - multigroup flux given on group boundaries ``energies`` in [eV]. - Reactions whose data cannot support the requested mode fall - back to ``'none'`` with a note in the mapping report. + multigroup flux given on group boundaries ``energies`` in [eV], + evaluating each group at its geometric midpoint clamped into + the tabulated range. Reactions whose data cannot support the + requested mode fall back to ``'none'`` with a note in the + mapping report. Modes other than ``'none'`` require + ``isomeric_branching=True``. .. versionadded:: 0.15.4 elis_rtol : float, optional @@ -378,9 +383,9 @@ def from_endf(cls, decay_files, fpy_files, neutron_files, states in the decay data. .. versionadded:: 0.15.4 - isomer_mapping_log : str, optional + isomer_mapping_log : str or os.PathLike, optional Path to write a report of every level-to-isomer mapping - decision. + decision. Requires ``isomeric_branching=True``. .. versionadded:: 0.15.4 @@ -401,11 +406,32 @@ def from_endf(cls, decay_files, fpy_files, neutron_files, """ transmutation_reactions = reactions - if (scalar_branching not in ('none', 'thermal') - and not isinstance(scalar_branching, tuple)): + if isinstance(scalar_branching, tuple): + if len(scalar_branching) != 2: + raise ValueError( + "scalar_branching tuple must be (energies, flux)") + group_energies = np.asarray(scalar_branching[0], dtype=float) + group_flux = np.asarray(scalar_branching[1], dtype=float) + if (group_energies.ndim != 1 or group_flux.ndim != 1 + or len(group_energies) != len(group_flux) + 1): + raise ValueError( + "scalar_branching group boundaries must be a 1-D array " + "one longer than the flux vector") + if group_flux.sum() <= 0.0: + raise ValueError( + "scalar_branching flux must have a positive sum") + scalar_branching = (group_energies, group_flux) + elif scalar_branching not in ('none', 'thermal'): raise ValueError( "scalar_branching must be 'none', 'thermal', or a tuple " "of (energies, flux)") + check_greater_than('elis_rtol', elis_rtol, 0.0) + if not isomeric_branching and ( + branching_files or scalar_branching != 'none' + or isomer_mapping_log is not None): + raise ValueError( + "branching_files, scalar_branching, and isomer_mapping_log " + "have no effect unless isomeric_branching=True") report = IsomerMappingReport(elis_rtol) if isomeric_branching else None # Create dictionary mapping target to filename @@ -549,11 +575,16 @@ def from_endf(cls, decay_files, fpy_files, neutron_files, # over supplementary files records = [] if isomeric_branching and daughter is not None: - available = iso_records.get(parent, {}) - if not (mts & set(available)): - available = iso_supplement.get(parent, {}) - for mt in sorted(mts & set(available)): - records.extend(available[mt]) + records, skipped_mts = select_records( + mts, iso_records.get(parent, {}), + iso_supplement.get(parent, {})) + if skipped_mts: + report.add_note( + parent, name, + 'MF=9/10 data on MT ' + + ', '.join(map(str, skipped_mts)) + + ' not used (only the first MT with ' + 'data is read)') if records: cls._add_isomeric_reactions( @@ -1314,6 +1345,17 @@ def set_branch_ratios(self, branch_ratios, reaction="(n,gamma)", "with a sum outside tolerance of 1 +/- {:5.3e}:\n{}".format( reaction, tolerance, "\n".join(tail))) + def inferred_ground_target(parent_name, new_ratios): + """Ground target auto-added when only metastables are given.""" + if (not all("_m" in t for t in new_ratios) + or sums[parent_name] == 1.0): + return None + ground_target = grounds.get(parent_name) + if ground_target is None: + pz, pa, pm = zam(parent_name) + ground_target = gnds_name(pz, pa + 1, 0) + return ground_target + # Check up front whether the rewrite would remove targets that carry # energy-dependent isomeric production data, before any mutation @@ -1323,11 +1365,8 @@ def set_branch_ratios(self, branch_ratios, reaction="(n,gamma)", new_ratios = branch_ratios[parent_name] kept = set(new_ratios) # Account for the ground target that will be added automatically - if all("_m" in t for t in new_ratios) and sums[parent_name] != 1.0: - ground_target = grounds.get(parent_name) - if ground_target is None: - pz, pa, pm = zam(parent_name) - ground_target = gnds_name(pz, pa + 1, 0) + ground_target = inferred_ground_target(parent_name, new_ratios) + if ground_target is not None: kept.add(ground_target) for ix in rxn_index: target = parent.reactions[ix].target @@ -1368,19 +1407,14 @@ def set_branch_ratios(self, branch_ratios, reaction="(n,gamma)", saved_production[popped.target] = data # Add new reactions - all_meta = True + ground_target = inferred_ground_target(parent_name, new_ratios) for target, br in new_ratios.items(): - all_meta = all_meta and ("_m" in target) parent.add_reaction(reaction, target, rxn_Q, br) # If branching ratios don't add to unity, add reaction to ground # with remainder of branching ratio - if all_meta and sums[parent_name] != 1.0: + if ground_target is not None: ground_br = 1.0 - sums[parent_name] - ground_target = grounds.get(parent_name) - if ground_target is None: - pz, pa, pm = zam(parent_name) - ground_target = gnds_name(pz, pa + 1, 0) new_ratios[ground_target] = ground_br parent.add_reaction(reaction, ground_target, rxn_Q, ground_br) diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index 1d386ccebab..18f00b0bc3b 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -96,7 +96,8 @@ class ProductionTable: mt : int ENDF reaction number of the section that supplied the data. source : str - Source library identifier, e.g. 'ENDF/B-VIII.1' or 'TENDL-2025'. + Source library identifier as recorded in the evaluation header, + e.g. 'ENDF/B-8.1' or 'TENDL-2023.1'. QM : float Mass-difference Q value in [eV], verbatim from the TAB1 header. QI : float @@ -302,7 +303,9 @@ class Nuclide: Dictionary mapping ``(reaction type, target)`` tuples to lists of :class:`IsomericProduction` instances carrying verbatim energy-dependent isomeric production data from the source - evaluations. + evaluations. Entries whose key does not correspond to a reaction + present on this nuclide are flagged by :meth:`validate` and are + not written by :meth:`to_xml_element`. .. versionadded:: 0.15.4 """ @@ -528,6 +531,7 @@ def to_xml_element(self): elem.append(src_elem) elem.set('reactions', str(len(self.reactions))) + written_production = set() for rx, daughter, Q, br in self.reactions: rx_elem = ET.SubElement(elem, 'reaction') rx_elem.set('type', rx) @@ -536,8 +540,13 @@ def to_xml_element(self): rx_elem.set('target', daughter) if br != 1.0: rx_elem.set('branching_ratio', str(br)) - for iso in self.isomeric_production.get((rx, daughter), []): - rx_elem.append(iso.to_xml_element()) + # Production data is written once per (type, target) pair even + # if the pair appears on several reaction elements + key = (rx, daughter) + if key not in written_production: + written_production.add(key) + for iso in self.isomeric_production.get(key, []): + rx_elem.append(iso.to_xml_element()) if self.yield_data: fpy_elem = ET.SubElement(elem, 'neutron_fission_yields') diff --git a/tests/unit_tests/test_deplete_isomeric_extraction.py b/tests/unit_tests/test_deplete_isomeric_extraction.py index c1cb08cacc9..e2c93dec6d1 100644 --- a/tests/unit_tests/test_deplete_isomeric_extraction.py +++ b/tests/unit_tests/test_deplete_isomeric_extraction.py @@ -10,11 +10,11 @@ import numpy as np import pytest -import openmc.deplete from openmc.deplete import Chain from openmc.deplete._isomeric import ( - FOLDED, GROUND, MATCHED, POSITIONAL, assign_levels, - compute_scalar_ratios, extract_isomeric_production, group_records) + FOLDED, GROUND, MATCHED, POSITIONAL, IsomerMappingReport, assign_levels, + compute_scalar_ratios, extract_isomeric_production, group_records, + select_records) AM241_MF8_MT102 = """\ 9.524100+4 2.389860+2 0 0 2 09543 8102 1 @@ -264,11 +264,137 @@ def test_scalar_ratios_none_mode(am241_records): assert ratios == {'Am242': 1.0, 'Am242_m1': 0.0} +def test_assign_positional_keeps_level_records_together(): + # One physical level carrying both an MF=9 and an MF=10 record must + # land on a single target through the positional fallback, even when + # more than one metastable is available + records = extract_isomeric_production(StubEvaluation({ + (8, 102): AM241_MF8_MT102, + (9, 102): AM241_MF9_MT102, + }), source='test')[102] + duplicate = extract_isomeric_production(StubEvaluation({ + (8, 102): AM241_MF8_MT102, + (9, 102): AM241_MF9_MT102, + }), source='test')[102][1] + duplicate.mf = 10 + both = [records[1], duplicate] + + isomers = {(95, 242): [(1, 10000.0), (2, 2000000.0)]} + assignments = assign_levels(both, 'Am242', isomers, 0.01) + assert all(a.status == POSITIONAL for a in assignments) + assert {a.target for a in assignments} == {'Am242_m1'} + + # With no metastables both records of the level fold together + assignments = assign_levels(both, 'Am242', {}, 0.01) + assert all(a.status == FOLDED for a in assignments) + assert {a.target for a in assignments} == {'Am242'} + + +def test_select_records(): + primary = {102: ['p102'], 600: ['p600']} + supplement = {102: ['s102']} + + # Primary wins over the supplement, first MT wins within a library, + # and additionally populated MTs are reported back + records, skipped = select_records({102, 600}, primary, supplement) + assert records == ['p102'] + assert skipped == [600] + + # Supplement is used only when the primary has nothing for the MTs + records, skipped = select_records({102}, {}, supplement) + assert records == ['s102'] + assert skipped == [] + + records, skipped = select_records({16}, primary, supplement) + assert records == [] + assert skipped == [] + + +def test_scalar_ratios_zero_bounded_group_structure(am241_records): + # A first group starting at 0 eV has a geometric midpoint of 0, which + # must be clamped into the MF=9 grid rather than silently dropping + # the whole thermal group from the collapse + assignments = assign_levels( + am241_records[102], 'Am242', {(95, 242): [(1, 48600.0)]}, 0.5) + grouped = group_records( + [(a.record, a.target) for a in assignments]) + mode = ([0.0, 0.625, 2.0e7], [0.9, 0.1]) + ratios, note = compute_scalar_ratios(grouped, 'Am242', mode) + assert note is None + assert ratios['Am242_m1'] == pytest.approx(0.104007, rel=1e-3) + assert ratios['Am242'] == pytest.approx(1.0 - ratios['Am242_m1']) + + +def test_report_write(tmp_path, am241_records): + report = IsomerMappingReport(elis_rtol=0.5) + assignments = assign_levels( + am241_records[102], 'Am242', {(95, 242): [(1, 48600.0)]}, 0.5) + for assignment in assignments: + report.add('Am241', '(n,gamma)', assignment, + assignment.target, assignment.status) + report.add_note('Am241', '(n,gamma)', 'example note') + path = tmp_path / 'mapping.log' + report.write(path) + text = path.read_text() + assert 'Am242_m1' in text + assert 'matched: 1' in text + assert 'ground: 1' in text + assert 'example note' in text + + +def test_fold_when_target_missing_from_decay_data(am241_records): + # A matched metastable whose nuclide is absent from the decay data is + # folded onto the ground-state target via replace_missing + class FakeHalfLife: + nominal_value = 1000.0 + + class FakeDecay: + def __init__(self, stable): + self.nuclide = {'stable': stable} + self.half_life = FakeHalfLife() + + from openmc.deplete import Nuclide + decay_data = {'Am241': FakeDecay(False), 'Am242': FakeDecay(False)} + nuclide = Nuclide('Am241') + report = IsomerMappingReport(elis_rtol=0.5) + Chain._add_isomeric_reactions( + nuclide, '(n,gamma)', 'Am242', 5537755.0, am241_records[102], + decay_data, {(95, 242): [(1, 48600.0)]}, 0.5, 'none', report) + + targets = [rx.target for rx in nuclide.reactions] + assert targets == ['Am242'] + productions = nuclide.isomeric_production[('(n,gamma)', 'Am242')] + assert sorted(p.level for p in productions) == [0, 2] + statuses = {row['status'] for row in report.rows} + assert FOLDED in statuses + + def test_from_endf_invalid_scalar_mode(): with pytest.raises(ValueError, match='scalar_branching'): Chain.from_endf([], [], [], scalar_branching='banana') +def test_from_endf_scalar_tuple_validation(): + with pytest.raises(ValueError, match='one longer'): + Chain.from_endf([], [], [], isomeric_branching=True, + scalar_branching=([1.0, 2.0], [1.0, 2.0])) + with pytest.raises(ValueError, match='positive sum'): + Chain.from_endf([], [], [], isomeric_branching=True, + scalar_branching=([1.0, 2.0], [0.0])) + with pytest.raises(ValueError, match='energies, flux'): + Chain.from_endf([], [], [], isomeric_branching=True, + scalar_branching=(1.0, 2.0, 3.0)) + + +def test_from_endf_kwargs_require_flag(): + with pytest.raises(ValueError, match='isomeric_branching=True'): + Chain.from_endf([], [], [], branching_files=['some_file']) + with pytest.raises(ValueError, match='isomeric_branching=True'): + Chain.from_endf([], [], [], isomer_mapping_log='log.txt') + with pytest.raises(ValueError, match='isomeric_branching=True'): + Chain.from_endf([], [], [], scalar_branching='thermal') + + def _find_endf(directory, patterns): for pattern in patterns: matches = sorted(Path(directory).glob(pattern)) @@ -303,3 +429,28 @@ def test_from_endf_isomeric(): assert production[0].tables[0].mf == 9 assert production[0].tables[0].data(0.0253) == pytest.approx(0.1) assert chain.validate(strict=True) + + +@pytest.mark.skipif( + 'OPENMC_ENDF_DATA' not in os.environ, + reason='OPENMC_ENDF_DATA environment variable must be set') +def test_from_endf_isomeric_no_data_identical(tmp_path): + # With isomeric_branching=True but no MF=8/9/10 data for the + # requested reaction, the chain must be identical to the default one + endf_data = Path(os.environ['OPENMC_ENDF_DATA']) + neutron = _find_endf(endf_data / 'neutrons', ['*Am*241*']) + decay_files = sorted((endf_data / 'decay').glob('*.endf')) + fpy = _find_endf(endf_data / 'nfy', ['*U*235*']) + + # Am241 carries MF=9 only for (n,gamma); (n,2n) has no isomeric data + common = dict(decay_files=decay_files, fpy_files=[fpy], + neutron_files=[neutron], reactions=['(n,2n)'], + progress=False) + chain_default = Chain.from_endf(**common) + chain_isomeric = Chain.from_endf(**common, isomeric_branching=True) + + default_xml = tmp_path / 'default.xml' + isomeric_xml = tmp_path / 'isomeric.xml' + chain_default.export_to_xml(default_xml) + chain_isomeric.export_to_xml(isomeric_xml) + assert default_xml.read_text() == isomeric_xml.read_text() diff --git a/tests/unit_tests/test_deplete_isomeric_production.py b/tests/unit_tests/test_deplete_isomeric_production.py index 5c5e17b2517..4cf60f72973 100644 --- a/tests/unit_tests/test_deplete_isomeric_production.py +++ b/tests/unit_tests/test_deplete_isomeric_production.py @@ -258,6 +258,46 @@ def test_no_isomeric_data(): assert elem.find("reaction").find("isomeric_production") is None +def test_duplicate_reaction_pair_roundtrip(tmp_path): + # Production data attached to a (type, target) pair that appears on + # more than one reaction element must not multiply on round trip + nuc = Nuclide("A") + nuc.add_reaction("(n,gamma)", "B", 0.0, 0.5) + nuc.add_reaction("(n,gamma)", "B", 0.0, 0.5) + table = ProductionTable( + mf=9, mt=102, source="test", QM=0.0, QI=0.0, + data=Tabulated1D([1.0, 2.0], [0.5, 0.5])) + nuc.isomeric_production[("(n,gamma)", "B")] = [ + IsomericProduction(0, 0.0, [table])] + + elem = nuc.to_xml_element() + assert len(elem.findall("reaction/isomeric_production")) == 1 + reread = Nuclide.from_xml(elem) + assert len(reread.isomeric_production[("(n,gamma)", "B")]) == 1 + + +def test_orphan_entry_not_written(): + nuc = Nuclide("A") + nuc.add_reaction("(n,gamma)", "B", 0.0, 1.0) + table = ProductionTable( + mf=9, mt=102, source="test", QM=0.0, QI=0.0, + data=Tabulated1D([1.0, 2.0], [0.5, 0.5])) + nuc.isomeric_production[("(n,gamma)", "B_m1")] = [ + IsomericProduction(1, 100.0, [table])] + elem = nuc.to_xml_element() + assert elem.findall("reaction/isomeric_production") == [] + + +def test_validate_reaction_sum_message(): + # A nuclide with no decay modes and inconsistent reaction sums must + # report the reaction sum (this previously raised a NameError) + nuc = Nuclide("A") + nuc.add_reaction("(n,gamma)", "B", 0.0, 0.5) + nuc.add_reaction("(n,gamma)", "B_m1", 0.0, 0.25) + with pytest.raises(ValueError, match="sum to 0.75"): + nuc.validate(strict=True) + + def test_get_isomeric_production(chain): data = chain.get_isomeric_production("Am241", "(n,gamma)") assert set(data) == {"Am242", "Am242_m1"}