diff --git a/docs/source/io_formats/depletion_chain.rst b/docs/source/io_formats/depletion_chain.rst index 74413e7b610..4f354106aad 100644 --- a/docs/source/io_formats/depletion_chain.rst +++ b/docs/source/io_formats/depletion_chain.rst @@ -96,6 +96,85 @@ 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 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 + + :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/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 new file mode 100644 index 00000000000..bbc8eddac50 --- /dev/null +++ b/openmc/deplete/_isomeric.py @@ -0,0 +1,473 @@ +"""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 Tabulated1D, 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 | 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): + """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 = None + decay_elis: float | None = 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. 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]] + 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 + + +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 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. + + 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: + return None + return float(table.data(THERMAL_ENERGY)) + + energies, flux = mode + energies = np.asarray(energies, dtype=float) + flux = np.asarray(flux, dtype=float) + # 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:]) + 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): + """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 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 + ---------- + 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 0.0253 eV; ' + '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 a2ff6dea0b6..dd60a2b7e77 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 @@ -22,6 +23,10 @@ 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, + select_records) from .._xml import get_text from .._sparse_compat import csc_array, dok_array import openmc.data @@ -313,14 +318,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 ---------- @@ -340,6 +347,47 @@ 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 named in the + evaluation header. Requires ``isomeric_branching=True``. + + .. 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], + 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 + 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 or os.PathLike, optional + Path to write a report of every level-to-isomer mapping + decision. Requires ``isomeric_branching=True``. + + .. versionadded:: 0.15.4 Returns ------- @@ -358,10 +406,39 @@ def from_endf(cls, decay_files, fpy_files, neutron_files, """ transmutation_reactions = reactions + 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 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 @@ -372,17 +449,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...') @@ -462,7 +570,29 @@ 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: + 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( + 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] @@ -532,8 +662,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. @@ -963,8 +1161,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 +1215,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 +1237,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 +1345,46 @@ 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 + + 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 + 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 + 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,27 +1396,33 @@ 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 + 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) + # 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 +1547,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 +1567,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 +1579,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/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index 95881483474..18f00b0bc3b 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,195 @@ 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 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 + 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 +299,15 @@ 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. 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 """ def __init__(self, name=None): @@ -122,6 +322,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 +474,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 @@ -322,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) @@ -330,6 +540,13 @@ def to_xml_element(self): rx_elem.set('target', daughter) if br != 1.0: rx_elem.set('branching_ratio', str(br)) + # 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') @@ -353,6 +570,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 +634,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 +643,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 +670,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_extraction.py b/tests/unit_tests/test_deplete_isomeric_extraction.py new file mode 100644 index 00000000000..e2c93dec6d1 --- /dev/null +++ b/tests/unit_tests/test_deplete_isomeric_extraction.py @@ -0,0 +1,456 @@ +"""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 + +from openmc.deplete import Chain +from openmc.deplete._isomeric import ( + 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 + 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_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)) + 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) + + +@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 new file mode 100644 index 00000000000..4cf60f72973 --- /dev/null +++ b/tests/unit_tests/test_deplete_isomeric_production.py @@ -0,0 +1,364 @@ +"""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 + + +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"} + 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