From 0e7dff3e91cbf68b92822d6e07ee18c84eabda2a Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 11 Jun 2025 17:59:55 +0300 Subject: [PATCH 01/22] photonuclear physics part 1, data handling only --- docs/source/io_formats/cross_sections.rst | 2 +- docs/source/io_formats/nuclear_data.rst | 52 +- docs/source/pythonapi/data.rst | 10 + include/openmc/cross_sections.h | 1 + include/openmc/secondary_kalbach.h | 1 + openmc/data/__init__.py | 1 + openmc/data/data.py | 3 + openmc/data/energy_distribution.py | 3 +- openmc/data/kalbach_mann.py | 81 +- openmc/data/library.py | 4 +- openmc/data/njoy.py | 127 ++++ openmc/data/photonuclear.py | 820 +++++++++++++++++++++ openmc/data/reaction.py | 11 +- src/cross_sections.cpp | 2 + tests/unit_tests/test_data_kalbach_mann.py | 18 +- tools/ci/download-xs.sh | 4 +- 16 files changed, 1082 insertions(+), 58 deletions(-) create mode 100644 openmc/data/photonuclear.py diff --git a/docs/source/io_formats/cross_sections.rst b/docs/source/io_formats/cross_sections.rst index 9f0759a3ac2..a0d03940c3e 100644 --- a/docs/source/io_formats/cross_sections.rst +++ b/docs/source/io_formats/cross_sections.rst @@ -50,7 +50,7 @@ attributes: :type: The type of data contained in the file. Accepted values are 'neutron', - 'thermal', 'photon', and 'wmp'. + 'thermal', 'photon', 'photonuclear' and 'wmp'. .. _depletion_element: diff --git a/docs/source/io_formats/nuclear_data.rst b/docs/source/io_formats/nuclear_data.rst index 8174108e1ae..b376abc4438 100644 --- a/docs/source/io_formats/nuclear_data.rst +++ b/docs/source/io_formats/nuclear_data.rst @@ -256,6 +256,54 @@ temperature-dependent data set. For example, the data set corresponding to :Groups: - **distribution** -- Format for angle-energy distributions are detailed in :ref:`angle_energy`. + +-------------------------- +Incident Photonuclear Data +-------------------------- + +**/** + +:Attributes: - **filetype** (*char[]*) -- String indicating the type of file + - **version** (*int[2]*) -- Major and minor version of the data + +**//** + +:Attributes: - **Z** (*int*) -- Atomic number + - **A** (*int*) -- Mass number. For a natural element, A=0 is given. + - **metastable** (*int*) -- Metastable state (0=ground, 1=first + excited, etc.) + - **atomic_weight_ratio** (*double*) -- Mass in units of neutron masses + - **n_reaction** (*int*) -- Number of reactions + +:Datasets: + - **energy** (*double[]*) -- Energies in [eV] at which cross sections + are tabulated + +**//** + +**//reactions/reaction_/** + +:Attributes: - **mt** (*int*) -- ENDF MT reaction number + - **label** (*char[]*) -- Name of the reaction + - **Q_value** (*double*) -- Q value in eV + - **center_of_mass** (*int*) -- Whether the reference frame for + scattering is center-of-mass (1) or laboratory (0) + - **n_product** (*int*) -- Number of reaction products + - **redundant** (*int*) -- Whether reaction is redundant + +**//reactions/reaction_/** + +:Datasets: + - **xs** (*double[]*) -- Cross section values tabulated against the + nuclide energy grid + + :Attributes: + - **threshold_idx** (*int*) -- Index on the energy + grid that the reaction threshold corresponds to + +**//reactions/reaction_/product_/** + + Reaction product data is described in :ref:`product`. .. _product: @@ -415,6 +463,7 @@ Kalbach-Mann :Object type: Group :Attributes: - **type** (*char[]*) -- 'kalbach-mann' + - **is_photon** (*bool*) -- Whether the incident particle is a photon. :Datasets: - **energy** (*double[]*) -- Incoming energies at which distributions exist :Attributes: @@ -586,7 +635,8 @@ Level Inelastic :Attributes: - **type** (*char[]*) -- 'level' - **threshold** (*double*) -- Energy threshold in the laboratory system in eV - - **mass_ratio** (*double*) -- :math:`(A/(A + 1))^2` + - **mass_ratio** (*double*) -- for incident neutrons: :math:`(A/(A + 1))^2`, + while for incident photons: :math:`(A-1)/A` Continuous Tabular ------------------ diff --git a/docs/source/pythonapi/data.rst b/docs/source/pythonapi/data.rst index 1eaf90c9724..b8771c55786 100644 --- a/docs/source/pythonapi/data.rst +++ b/docs/source/pythonapi/data.rst @@ -50,6 +50,15 @@ The following classes are used for storing thermal neutron scattering data: CoherentElastic IncoherentElastic +The following classes are used for storing incident photonuclear interaction data: + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: myclass.rst + + IncidentPhotonuclear + PhotonuclearReaction Core Functions -------------- @@ -213,3 +222,4 @@ NJOY Interface njoy.make_pendf njoy.make_ace njoy.make_ace_thermal + njoy.make_ace_photonuclear diff --git a/include/openmc/cross_sections.h b/include/openmc/cross_sections.h index 06140a6a8cf..9620cd5e1e2 100644 --- a/include/openmc/cross_sections.h +++ b/include/openmc/cross_sections.h @@ -20,6 +20,7 @@ class Library { enum class Type { neutron = 1, photon = 3, + photonuclear = 6, thermal = 2, multigroup = 4, wmp = 5 diff --git a/include/openmc/secondary_kalbach.h b/include/openmc/secondary_kalbach.h index 83806d35248..a1a27dc85b9 100644 --- a/include/openmc/secondary_kalbach.h +++ b/include/openmc/secondary_kalbach.h @@ -44,6 +44,7 @@ class KalbachMann : public AngleEnergy { xt::xtensor a; //!< Parameterized function }; + bool is_photon_; //!< Whether the projectile is a photon int n_region_; //!< Number of interpolation regions vector breakpoints_; //!< Breakpoints between regions vector interpolation_; //!< Interpolation laws diff --git a/openmc/data/__init__.py b/openmc/data/__init__.py index c2d35565a8a..a6fdf1a35a2 100644 --- a/openmc/data/__init__.py +++ b/openmc/data/__init__.py @@ -12,6 +12,7 @@ from .data import * from .neutron import * from .photon import * +from .photonuclear import * from .decay import * from .reaction import * from . import ace diff --git a/openmc/data/data.py b/openmc/data/data.py index 2142a5dc905..cd5bc19f238 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -281,6 +281,9 @@ # Neutron mass in units of amu NEUTRON_MASS = 1.00866491595 +# Neutron mass in units of eV/c^2 +NEUTRON_MASS_EV = 939.56542052e6 + # Used in atomic_mass function as a cache _ATOMIC_MASS: dict[str, float] = {} diff --git a/openmc/data/energy_distribution.py b/openmc/data/energy_distribution.py index 069ab1b9beb..bb04e134097 100644 --- a/openmc/data/energy_distribution.py +++ b/openmc/data/energy_distribution.py @@ -1253,7 +1253,8 @@ def from_ace(cls, ace, idx, ldis): # Create continuous distribution eout_continuous = Tabular(data[0][n_discrete_lines:], data[1][n_discrete_lines:]/EV_PER_MEV, - INTERPOLATION_SCHEME[intt]) + INTERPOLATION_SCHEME[intt], + ignore_negative=True) eout_continuous.c = data[2][n_discrete_lines:] # If discrete lines are present, create a mixture distribution diff --git a/openmc/data/kalbach_mann.py b/openmc/data/kalbach_mann.py index d92bf9c213a..d5cc3f37f99 100644 --- a/openmc/data/kalbach_mann.py +++ b/openmc/data/kalbach_mann.py @@ -9,7 +9,7 @@ from openmc.stats import Tabular, Univariate, Discrete, Mixture from .function import Tabulated1D, INTERPOLATION_SCHEME from .angle_energy import AngleEnergy -from .data import EV_PER_MEV +from .data import EV_PER_MEV, NEUTRON_MASS_EV from .endf import get_list_record, get_tab2_record @@ -204,13 +204,14 @@ def kalbach_slope(energy_projectile, energy_emitted, za_projectile, Kalbach-Mann slope given with the same format as ACE file. """ - # TODO: develop for photons as projectile - # TODO: test for other particles than neutron - if za_projectile != 1: - raise NotImplementedError( - "Developed and tested for neutron projectile only." - ) + if za_projectile == 0: + # Calculate slope for photons using Eq. 3 in doi:10.1080/18811248.1995.9731830 + # or ENDF-6 Formats Manual section 6.2.3.2 + slope_n = kalbach_slope(energy_projectile, energy_emitted, 1, + za_emitted, za_target) + return slope_n * np.sqrt(0.5*energy_projectile/NEUTRON_MASS_EV)*np.minimum(4,np.maximum(1,9.3/np.sqrt(energy_emitted*1e-6))) + # Special handling of elemental carbon if za_emitted == 6000: za_emitted = 6012 @@ -268,6 +269,8 @@ class KalbachMann(AngleEnergy): slope : Iterable of openmc.data.Tabulated1D Kalbach-Chadwick angular distribution slope value 'a' as a function of outgoing energy for each incoming energy + is_photon : bool + Whether projectile is a photon, defaults to False Attributes ---------- @@ -285,14 +288,17 @@ class KalbachMann(AngleEnergy): slope : Iterable of openmc.data.Tabulated1D Kalbach-Chadwick angular distribution slope value 'a' as a function of outgoing energy for each incoming energy + is_photon : bool + Whether projectile particle is a photon """ def __init__(self, breakpoints, interpolation, energy, energy_out, - precompound, slope): + precompound, slope, is_photon=False): super().__init__() self.breakpoints = breakpoints self.interpolation = interpolation + self.is_photon = is_photon self.energy = energy self.energy_out = energy_out self.precompound = precompound @@ -317,6 +323,15 @@ def interpolation(self, interpolation): cv.check_type('Kalbach-Mann interpolation', interpolation, Iterable, Integral) self._interpolation = interpolation + + @property + def is_photon(self): + return self._is_photon + + @is_photon.setter + def is_photon(self, is_photon): + cv.check_type('Kalbach-Mann is_photon', is_photon, bool) + self._is_photon = is_photon @property def energy(self): @@ -367,6 +382,7 @@ def to_hdf5(self, group): """ group.attrs['type'] = np.bytes_('kalbach-mann') + group.attrs['is_photon'] = self.is_photon dset = group.create_dataset('energy', data=self.energy) dset.attrs['interpolation'] = np.vstack((self.breakpoints, @@ -436,6 +452,7 @@ def from_hdf5(cls, group): Kalbach-Mann energy distribution """ + is_photon = bool(group.attrs.get("is_photon", False)) interp_data = group['energy'].attrs['interpolation'] energy_breakpoints = interp_data[0, :] energy_interpolation = interp_data[1, :] @@ -491,7 +508,7 @@ def from_hdf5(cls, group): slope.append(km_a) return cls(energy_breakpoints, energy_interpolation, - energy, energy_out, precompound, slope) + energy, energy_out, precompound, slope, is_photon=is_photon) @classmethod def from_ace(cls, ace, idx, ldis): @@ -507,6 +524,8 @@ def from_ace(cls, ace, idx, ldis): ldis : int Index in XSS array of the start of the energy distribution block (e.g. JXS[11]) + is_photon : bool + Whether projectile is photon Returns ------- @@ -514,6 +533,7 @@ def from_ace(cls, ace, idx, ldis): Kalbach-Mann energy-angle distribution """ + is_photon = bool(ace.data_type.value == 'u') # Read number of interpolation regions and incoming energies n_regions = int(ace.xss[idx]) n_energy_in = int(ace.xss[idx + 1 + 2*n_regions]) @@ -586,10 +606,10 @@ def from_ace(cls, ace, idx, ldis): km_r.append(Tabulated1D(data[0], data[3])) km_a.append(Tabulated1D(data[0], data[4])) - return cls(breakpoints, interpolation, energy, energy_out, km_r, km_a) + return cls(breakpoints, interpolation, energy, energy_out, km_r, km_a, is_photon=is_photon) @classmethod - def from_endf(cls, file_obj, za_emitted, za_target, projectile_mass): + def from_endf(cls, file_obj, za_emitted, za_target, za_projectile): """Generate Kalbach-Mann distribution from an ENDF evaluation. If the projectile is a neutron, the slope is calculated when it is @@ -606,14 +626,8 @@ def from_endf(cls, file_obj, za_emitted, za_target, projectile_mass): ZA identifier of the emitted particle za_target : int ZA identifier of the target - projectile_mass : float - Mass of the projectile - - Warns - ----- - UserWarning - If the mass of the projectile is not equal to 1 (other than - a neutron), the slope is not calculated and set to 0 if missing. + za_projectile : int + ZA identifier of the projectile Returns ------- @@ -621,6 +635,7 @@ def from_endf(cls, file_obj, za_emitted, za_target, projectile_mass): Kalbach-Mann energy-angle distribution """ + is_photon = bool(za_projectile==0) params, tab2 = get_tab2_record(file_obj) lep = params[3] ne = params[5] @@ -654,31 +669,19 @@ def from_endf(cls, file_obj, za_emitted, za_target, projectile_mass): a_i = values[:, 3] calculated_slope.append(False) else: - # Check if the projectile is not a neutron - if not np.isclose(projectile_mass, 1.0, atol=1.0e-12, rtol=0.): - warn( - "Kalbach-Mann slope calculation is only available with " - "neutrons as projectile. Slope coefficients are set to 0." - ) - a_i = np.zeros_like(r_i) - calculated_slope.append(False) - - else: - # TODO: retrieve ZA of the projectile - za_projectile = 1 - a_i = [kalbach_slope(energy_projectile=energy[i], - energy_emitted=e, - za_projectile=za_projectile, - za_emitted=za_emitted, - za_target=za_target) - for e in eout_i] - calculated_slope.append(True) + a_i = [kalbach_slope(energy_projectile=energy[i], + energy_emitted=e, + za_projectile=za_projectile, + za_emitted=za_emitted, + za_target=za_target) + for e in eout_i] + calculated_slope.append(True) precompound.append(Tabulated1D(eout_i, r_i)) slope.append(Tabulated1D(eout_i, a_i)) km_distribution = cls(tab2.breakpoints, tab2.interpolation, energy, - energy_out, precompound, slope) + energy_out, precompound, slope, is_photon) # List of bool to indicate slope calculation by OpenMC km_distribution._calculated_slope = calculated_slope diff --git a/openmc/data/library.py b/openmc/data/library.py index a6ce1bbd32c..2bcd9877a99 100644 --- a/openmc/data/library.py +++ b/openmc/data/library.py @@ -37,7 +37,7 @@ def get_by_material(self, name, data_type='neutron'): name : str Name of material, e.g. 'Am241' data_type : str - Name of data type, e.g. 'neutron', 'photon', 'wmp', or 'thermal' + Name of data type, e.g. 'neutron', 'photon', 'photonuclear', 'wmp', or 'thermal' .. versionadded:: 0.12 @@ -61,7 +61,7 @@ def remove_by_material(self, name: str, data_type='neutron'): name : str Name of material, e.g. 'Am241' data_type : str - Name of data type, e.g. 'neutron', 'photon', 'wmp', or 'thermal' + Name of data type, e.g. 'neutron', 'photon', 'photonuclear', 'wmp', or 'thermal' """ library = self.get_by_material(name, data_type) diff --git a/openmc/data/njoy.py b/openmc/data/njoy.py index 4c538bd6e6f..ffa7389cb15 100644 --- a/openmc/data/njoy.py +++ b/openmc/data/njoy.py @@ -221,6 +221,13 @@ 222 64 {mt_elastic} {elastic_type} {data.nmix} {energy_max} {iwt}/ """ +_PHOTONUCLEAR_TEMPLATE_ACER = """ +acer / %%%%%%%%%%%%%%%%%%%%%%%% Write out in ACE format %%%%%%%%%%%%%%%%%%%%%%%% +{nendf} {nacer_in} 0 {nace} {ndir} +5 0 1 .{ext} / +'{library}: {zsymam}'/ +{mat} +""" def run(commands, tapein, tapeout, input_filename=None, stdout=False, njoy_exec='njoy'): @@ -673,3 +680,123 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None, for temperature in temperatures: (output_dir / f"ace_{temperature:.1f}").unlink() (output_dir / f"xsdir_{temperature:.1f}").unlink() + + +def make_ace_photonuclear(filename, acer=True, xsdir=None, + output_dir=None, pendf=False, error=0.001, + evaluation=None, **kwargs): + """Generate incident photonuclear ACE file from an ENDF file + + File names can be passed to + ``[acer, xsdir, pendf]`` + to specify the exact output for the given module. + Otherwise, the files will be writen to the current directory + or directory specified by ``output_dir``. Default file + names mirror the variable names, e.g. ``xsdir`` output + will be written to a file named ``xsdir`` unless otherwise + specified. + + Parameters + ---------- + filename : str + Path to ENDF file + acer : bool or str, optional + Flag indicating if acer should be run. If a string is give, write the + resulting ``ace`` file to this location. Path of ACE file to write. + Defaults to ``"ace"`` + xsdir : str, optional + Path of xsdir file to write. Defaults to ``"xsdir"`` in the same + directory as ``acer`` + output_dir : str, optional + Directory to write output for requested modules. If not provided + and at least one of ``[pendf, broadr, heatr, gaspr, purr, acer]`` + is ``True``, then write output files to current directory. If given, + must be a path to a directory. + pendf : str, optional + Path of pendf file to write. If omitted, the pendf file is not saved. + error : float, optional + Fractional error tolerance for NJOY processing + evaluation : openmc.data.endf.Evaluation, optional + If the ENDF file contains multiple material evaluations, this argument + indicates which evaluation should be used. + **kwargs + Keyword arguments passed to :func:`openmc.data.njoy.run` + + Raises + ------ + subprocess.CalledProcessError + If the NJOY process returns with a non-zero status + IOError + If ``output_dir`` does not point to a directory + + """ + if output_dir is None: + output_dir = Path() + else: + output_dir = Path(output_dir) + if not output_dir.is_dir(): + raise IOError(f"{output_dir} is not a directory") + + ev = evaluation if evaluation is not None else endf.Evaluation(filename) + mat = ev.material + zsymam = ev.target['zsymam'] + + # Determine name of library + library = '{}-{}.{}'.format(*ev.info['library']) + + # Create njoy commands by modules + commands = "" + + nendf, npendf = 20, 21 + tapein = {nendf: filename} + tapeout = {} + if pendf: + tapeout[npendf] = (output_dir / "pendf") if pendf is True else pendf + + # reconr + commands += _TEMPLATE_RECONR + nlast = npendf + + commands = commands.format(**locals()) + + # acer + if acer: + nacer_in = nlast + # Extend input with an ACER run + nace = nacer_in + 1 + ndir = nace + 1 + ext = f'{1:02}' + + commands += _PHOTONUCLEAR_TEMPLATE_ACER.format(**locals()) + + # Indicate tapes to save for each ACER run + tapeout[nace] = output_dir / f"ace_0.0" + tapeout[ndir] = output_dir / f"xsdir_0.0" + commands += 'stop\n' + run(commands, tapein, tapeout, **kwargs) + + if acer: + ace = (output_dir / "ace") if acer is True else Path(acer) + xsdir = (ace.parent / "xsdir") if xsdir is None else xsdir + with ace.open('w') as ace_file, xsdir.open('w') as xsdir_file: + # Get contents of ACE file + text = (output_dir / f"ace_0.0").read_text() + + # If the target is metastable, make sure that ZAID in the ACE + # file reflects this by adding 400 + if ev.target['isomeric_state'] > 0: + mass_first_digit = int(text[3]) + if mass_first_digit <= 2: + text = text[:3] + str(mass_first_digit + 4) + text[4:] + + # Concatenate into destination ACE file + ace_file.write(text) + + # Concatenate into destination xsdir file + xsdir_in = output_dir / f"xsdir_0.0" + xsdir_file.write(xsdir_in.read_text()) + + # Remove ACE/xsdir files + (output_dir / f"ace_0.0").unlink() + (output_dir / f"xsdir_0.0").unlink() + diff --git a/openmc/data/photonuclear.py b/openmc/data/photonuclear.py new file mode 100644 index 00000000000..810a01c07e5 --- /dev/null +++ b/openmc/data/photonuclear.py @@ -0,0 +1,820 @@ +from collections import OrderedDict +from collections.abc import Callable, Iterable +from io import StringIO +from numbers import Integral, Real +from warnings import warn +import os +import tempfile + +import h5py +import numpy as np + +from openmc.mixin import EqualityMixin +from openmc.stats import Uniform +import openmc.checkvalue as cv +from . import HDF5_VERSION, HDF5_VERSION_MAJOR +from .ace import Table, get_metadata, get_table, Library +from .angle_distribution import AngleDistribution +from .angle_energy import AngleEnergy +from .correlated import CorrelatedAngleEnergy +from .data import ATOMIC_SYMBOL, EV_PER_MEV +from .endf import Evaluation, get_head_record, get_tab1_record, get_cont_record +from .function import Tabulated1D +from .njoy import make_ace_photonuclear +from .reaction import Reaction, REACTION_NAME, FISSION_MTS, _get_products, _get_fission_products_endf, _get_photon_products_endf, _get_activation_products +from .product import Product +from .energy_distribution import EnergyDistribution, LevelInelastic, \ + DiscretePhoton +from .function import Tabulated1D, Polynomial +from .kalbach_mann import KalbachMann +from .laboratory import LaboratoryAngleEnergy +from .nbody import NBodyPhaseSpace +from .product import Product +from .uncorrelated import UncorrelatedAngleEnergy + +_REACTION_NAME = {key:value.replace("(n,","(gamma,") for key,value in REACTION_NAME.items()} + +class PhotonuclearReaction(EqualityMixin): + def __init__(self, mt): + self._center_of_mass = True + self._redundant = False + self._q_value = 0. + self._xs = None + self._products = [] + self._derived_products = [] + self.mt = mt + + def __repr__(self): + if self.mt in _REACTION_NAME: + return f"" + else: + return f"" + + @property + def center_of_mass(self): + return self._center_of_mass + + @center_of_mass.setter + def center_of_mass(self, center_of_mass): + cv.check_type('center of mass', center_of_mass, (bool, np.bool_)) + self._center_of_mass = center_of_mass + + @property + def redundant(self): + return self._redundant + + @redundant.setter + def redundant(self, redundant): + cv.check_type('redundant', redundant, (bool, np.bool_)) + self._redundant = redundant + + @property + def q_value(self): + return self._q_value + + @q_value.setter + def q_value(self, q_value): + cv.check_type('Q value', q_value, Real) + self._q_value = q_value + + @property + def products(self): + return self._products + + @products.setter + def products(self, products): + cv.check_type('reaction products', products, Iterable, Product) + self._products = products + + @property + def derived_products(self): + return self._derived_products + + @derived_products.setter + def derived_products(self, derived_products): + cv.check_type('reaction derived products', derived_products, + Iterable, Product) + self._derived_products = derived_products + + @property + def xs(self): + return self._xs + + @xs.setter + def xs(self, xs): + cv.check_type("reaction cross section", xs, Callable) + self._xs = xs + + @classmethod + def from_endf(cls, ev, mt): + """Generate a reaction from an ENDF evaluation + + Parameters + ---------- + ev : openmc.data.endf.Evaluation + ENDF evaluation + mt : int + The MT value of the reaction to get data for + + Returns + ------- + rx : openmc.data.Reaction + Reaction data + + """ + rx = cls(mt) + + # Integrated cross section + if (3, mt) in ev.section: + file_obj = StringIO(ev.section[3, mt]) + get_head_record(file_obj) + params, rx.xs = get_tab1_record(file_obj) + rx.q_value = params[1] + + # Get fission product yields (nu) as well as delayed neutron energy + # distributions + if mt in FISSION_MTS: + rx.products, rx.derived_products = _get_fission_products_endf(ev) + + if (6, mt) in ev.section: + # Product angle-energy distribution + for product in _get_products(ev, mt): + if mt in FISSION_MTS and product.particle == 'neutron': + rx.products[0].applicability = product.applicability + rx.products[0].distribution = product.distribution + else: + rx.products.append(product) + elif (4, mt) in ev.section or (5, mt) in ev.section: + # Uncorrelated angle-energy distribution + + # Note that the energy distribution for MT=455 is read in + # _get_fission_products_endf rather than here + if (5, mt) in ev.section: + product = Product('photon') + file_obj = StringIO(ev.section[5, mt]) + items = get_head_record(file_obj) + nk = items[4] + for i in range(nk): + params, applicability = get_tab1_record(file_obj) + dist = UncorrelatedAngleEnergy() + dist.energy = EnergyDistribution.from_endf(file_obj, params) + + product.applicability.append(applicability) + product.distribution.append(dist) + elif mt == 2: + # Elastic scattering -- no energy distribution is given since it + # can be calulcated analytically + product = Product('photon') + dist = UncorrelatedAngleEnergy() + product.distribution.append(dist) + elif mt >= 50 and mt < 91: + # Level inelastic scattering -- no energy distribution is given + # since it can be calculated analytically. Here we determine the + # necessary parameters to create a LevelInelastic object + product = Product('neutron') + dist = UncorrelatedAngleEnergy() + + A = ev.target['mass'] + threshold = abs(rx.q_value) + mass_ratio = (A-1)/A + dist.energy = LevelInelastic(threshold, mass_ratio) + + product.distribution.append(dist) + + if (4, mt) in ev.section: + for dist in product.distribution: + dist.angle = AngleDistribution.from_endf(ev, mt) + + if mt in FISSION_MTS and (5, mt) in ev.section: + # For fission reactions, + rx.products[0].applicability = product.applicability + rx.products[0].distribution = product.distribution + else: + rx.products.append(product) + + if (8, mt) in ev.section: + rx.products += _get_activation_products(ev, rx) + + if (12, mt) in ev.section or (13, mt) in ev.section: + rx.products += _get_photon_products_endf(ev, rx) + return rx + + def to_hdf5(self, group): + """Write reaction to an HDF5 group + + Parameters + ---------- + group : h5py.Group + HDF5 group to write to + + """ + + group.attrs['mt'] = self.mt + if self.mt in _REACTION_NAME: + group.attrs['label'] = np.bytes_(_REACTION_NAME[self.mt]) + else: + group.attrs['label'] = np.bytes_(self.mt) + group.attrs['Q_value'] = self.q_value + group.attrs['center_of_mass'] = 1 if self.center_of_mass else 0 + group.attrs['redundant'] = 1 if self.redundant else 0 + + dset = group.create_dataset('xs', data=self.xs.y) + threshold_idx = getattr(self.xs, '_threshold_idx', 0) + dset.attrs['threshold_idx'] = threshold_idx + + for i, p in enumerate(self.products): + pgroup = group.create_group(f'product_{i}') + p.to_hdf5(pgroup) + + @classmethod + def from_hdf5(cls, group, energy): + """Generate reaction from an HDF5 group + + Parameters + ---------- + group : h5py.Group + HDF5 group to read from + energy : array + array of energies at which cross sections are tabulated at. + + Returns + ------- + openmc.data.Reaction + Reaction data + + """ + + mt = group.attrs['mt'] + rx = cls(mt) + rx.q_value = group.attrs['Q_value'] + rx.center_of_mass = bool(group.attrs['center_of_mass']) + rx.redundant = bool(group.attrs.get('redundant', False)) + xs = group['xs'][()] + threshold_idx = group['xs'].attrs['threshold_idx'] + tabulated_xs = Tabulated1D(energy[threshold_idx:], xs) + tabulated_xs._threshold_idx = threshold_idx + rx.xs = tabulated_xs + + # Determine number of products + n_product = 0 + for name in group: + if name.startswith('product_'): + n_product += 1 + + # Read reaction products + for i in range(n_product): + pgroup = group[f'product_{i}'] + rx.products.append(Product.from_hdf5(pgroup)) + + return rx + + @classmethod + def from_ace(cls, ace, i_reaction): + # Get nuclide energy grid + n_grid = ace.nxs[3] + grid = ace.xss[ace.jxs[1] : ace.jxs[1] + n_grid] * EV_PER_MEV + + if i_reaction > 0: + mt = int(ace.xss[ace.jxs[6] + i_reaction - 1]) + rx = cls(mt) + + # Get Q-value of reaction + rx.q_value = ace.xss[ace.jxs[7] + i_reaction - 1]*EV_PER_MEV + + # ================================================================== + # CROSS SECTION + + # Get locator for cross-section data + loc = int(ace.xss[ace.jxs[8] + i_reaction - 1]) + + # Determine starting index on energy grid + threshold_idx = int(ace.xss[ace.jxs[9] + loc - 1]) - 1 + + # Determine number of energies in reaction + n_energy = int(ace.xss[ace.jxs[9] + loc]) + energy = grid[threshold_idx : threshold_idx + n_energy] + + # Read reaction cross section + xs = ace.xss[ace.jxs[9] + loc + 1 : ace.jxs[9] + loc + 1 + n_energy] + + # Fix negatives -- known issue for Y89 in JEFF 3.2 + if np.any(xs < 0.0): + warn( + "Negative cross sections found for MT={} in {}. Setting " + "to zero.".format(rx.mt, ace.name) + ) + xs[xs < 0.0] = 0.0 + + tabulated_xs = Tabulated1D(energy, xs) + tabulated_xs._threshold_idx = threshold_idx + rx.xs = tabulated_xs + + # ================================================================== + # YIELD AND ANGLE-ENERGY DISTRIBUTION + + for i_typ in range(ace.nxs[5]): + loc = ace.jxs[10]+i_typ*ace.nxs[7]-1 + mts = ace.xss[int(ace.xss[loc+5]):int(ace.xss[loc+5])+int(ace.xss[loc+2])].astype(int) + try: + i_mtr = mts.tolist().index(mt) + except ValueError: + #skip products for other reactions + continue + match int(ace.xss[loc+1]): + case 1: + particle=Product('neutron') + case 2: + particle=Product('photon') + case _: + #TODO: support more product particles + # for now skip particle if it is not a neutron or photon + continue + + # Determine reference frame + ty = int(ace.xss[int(ace.xss[loc+6])+i_mtr]) + rx.center_of_mass = (ty < 0) + + # Determine multiplicity + idx = int(ace.xss[loc+8])+int(ace.xss[int(ace.xss[loc+7])+i_mtr])-1 + match int(ace.xss[idx]): + case 6 | 16 | 12: + assert int(ace.xss[idx+1]) == mt + yield_ = Tabulated1D.from_ace(ace, idx+2) + case _: + raise NotImplementedError('partial yields not implemented yet') + particle.yield_ = yield_ + + + # Determine locator for energy distribution + idx = int(ace.xss[int(ace.xss[loc+11])+i_mtr]) + distribution = AngleEnergy.from_ace(ace, int(ace.xss[loc+12]), idx) + + # Determine locator for angular distribution + idx = int(ace.xss[int(ace.xss[loc+9])+i_mtr]) + if idx==0: + # No angular distribution data are given for this reaction, + # isotropic scattering is asssumed in LAB + energy = np.array([particle.yield_.x[0], particle.yield_.x[-1]]) + mu_isotropic = Uniform(-1., 1.) + distribution.angle = AngleDistribution( + energy, [mu_isotropic, mu_isotropic]) + elif idx==-1: + pass + else: + assert idx>=0 + distribution.angle = AngleDistribution.from_ace(ace, int(ace.xss[loc+10]), idx) + + particle.distribution.append(distribution) + rx.products.append(particle) + else: + # elastic cross section + mt = 2 + rx = cls(mt) + + # Get elastic cross section values + elastic_xs = ace.xss[ace.jxs[4] : ace.jxs[4] + ace.nxs[3]] + + # Fix negatives -- known issue for Ti46,49,50 in JEFF 3.2 + if np.any(elastic_xs < 0.0): + warn( + "Negative elastic scattering cross section found for {}. " + "Setting to zero.".format(ace.name) + ) + elastic_xs[elastic_xs < 0.0] = 0.0 + + tabulated_xs = Tabulated1D(grid, elastic_xs) + tabulated_xs._threshold_idx = 0 + rx.xs = tabulated_xs + + return rx + +class IncidentPhotonuclear(EqualityMixin): + """photo-nuclear interaction data. + + This class stores data derived from an ENDF-6 format photo-nuclear interaction + sublibrary. Instances of this class are not normally instantiated by the + user but rather created using the factory methods + :meth:`Photonuclear.from_hdf5`, :meth:`Photonuclear.from_ace`, and + :meth:`Photonuclear.from_endf`. + + Parameters + ---------- + name : str + Name of the nuclide using the GND naming convention + atomic_number : int + Number of photo-nuclears in the target nucleus + atomic_number : int + Number of photo-nuclears in the target nucleus + mass_number : int + Number of nucleons in the target nucleus + metastable : int + Metastable state of the target nucleus. A value of zero indicates ground + state. + atomic_weight_ratio : float + Atomic weight ratio of the target nuclide. + + Attributes + ---------- + atomic_number : int + Number of photo-nuclears in the target nucleus + atomic_symbol : str + Atomic symbol of the nuclide, e.g., 'Zr' + atomic_weight_ratio : float + Atomic weight ratio of the target nuclide. + mass_number : int + Number of nucleons in the target nucleus + metastable : int + Metastable state of the target nucleus. A value of zero indicates ground + state. + name : str + Name of the nuclide using the GND naming convention + reactions : collections.OrderedDict + Contains the cross sections, secondary angle and energy distributions, + and other associated data for each reaction. The keys are the MT values + and the values are Reaction objects. + + """ + + def __init__( + self, name, atomic_number, mass_number, metastable, atomic_weight_ratio + ): + self.name = name + self.atomic_number = atomic_number + self.mass_number = mass_number + self.reactions = OrderedDict() + self.energy = [] + self.metastable = metastable + self.atomic_weight_ratio = atomic_weight_ratio + + def __contains__(self, mt): + return mt in self.reactions + + def __getitem__(self, mt): + if mt in self.reactions: + return self.reactions[mt] + else: + raise KeyError("No reaction with MT={}.".format(mt)) + + def __repr__(self): + return "".format(self.name) + + def __iter__(self): + return iter(self.reactions.values()) + + @property + def atomic_number(self): + return self._atomic_number + + @property + def name(self): + return self._name + + @property + def mass_number(self): + return self._mass_number + + @property + def metastable(self): + return self._metastable + + @property + def atomic_weight_ratio(self): + return self._atomic_weight_ratio + + @property + def atomic_symbol(self): + return ATOMIC_SYMBOL[self.atomic_number] + + @name.setter + def name(self, name): + cv.check_type("name", name, str) + self._name = name + + @atomic_number.setter + def atomic_number(self, atomic_number): + cv.check_type("atomic number", atomic_number, Integral) + cv.check_greater_than("atomic number", atomic_number, 0, True) + self._atomic_number = atomic_number + + @mass_number.setter + def mass_number(self, mass_number): + cv.check_type("mass number", mass_number, Integral) + cv.check_greater_than("mass number", mass_number, 0, True) + self._mass_number = mass_number + + @metastable.setter + def metastable(self, metastable): + cv.check_type("metastable", metastable, Integral) + cv.check_greater_than("metastable", metastable, 0, True) + self._metastable = metastable + + @atomic_weight_ratio.setter + def atomic_weight_ratio(self, atomic_weight_ratio): + cv.check_type("atomic weight ratio", atomic_weight_ratio, Real) + cv.check_greater_than("atomic weight ratio", atomic_weight_ratio, 0.0) + self._atomic_weight_ratio = atomic_weight_ratio + + @classmethod + def from_endf(cls, ev_or_filename): + """Generate incident photo-nuclear data from an ENDF evaluation + + Parameters + ---------- + ev_or_filename : openmc.data.endf.Evaluation or str + ENDF evaluation to read from. If given as a string, it is assumed to + be the filename for the ENDF file. + + + Returns + ------- + openmc.data.Photonuclear + photo-nuclear interaction data + + """ + + if isinstance(ev_or_filename, Evaluation): + ev = ev_or_filename + else: + ev = Evaluation(ev_or_filename) + + atomic_number = ev.target["atomic_number"] + mass_number = ev.target["mass_number"] + metastable = ev.target["isomeric_state"] + atomic_weight_ratio = ev.target["mass"] + + element = ATOMIC_SYMBOL[atomic_number] + name = "{}{}".format(element, mass_number) + + data = cls(name, atomic_number, mass_number, metastable, atomic_weight_ratio) + + # Read each reaction + for mf, mt, nc, mod in ev.reaction_list: + if mf == 3: + data.reactions[mt] = PhotonuclearReaction.from_endf(ev, mt) + + data._evaluation = ev + return data + + @classmethod + def from_ace(cls, ace_or_filename, metastable_scheme="nndc"): + """Generate incident photo-nuclear continuous-energy data from an ACE table + + Parameters + ---------- + ace_or_filename : openmc.data.ace.Table or str + ACE table to read from. If the value is a string, it is assumed to + be the filename for the ACE file. + metastable_scheme : {'nndc', 'mcnp'} + Determine how ZAID identifiers are to be interpreted in the case of + a metastable nuclide. Because the normal ZAID (=1000*Z + A) does not + encode metastable information, different conventions are used among + different libraries. In MCNP libraries, the convention is to add 400 + for a metastable nuclide except for Am242m, for which 95242 is + metastable and 95642 (or 1095242 in newer libraries) is the ground + state. For NNDC libraries, ZAID is given as 1000*Z + A + 100*m. + + Returns + ------- + openmc.data.Photonuclear + Incident photo-nuclear continuous-energy data + + """ + + # First obtain the data for the first provided ACE table/file + if isinstance(ace_or_filename, Table): + ace = ace_or_filename + else: + ace = get_table(ace_or_filename) + + # If mass number hasn't been specified, make an educated guess + zaid, xs = ace.name.split(".") + if not xs.endswith("u"): + raise TypeError( + "{} is not a continuous-energy photo-nuclear ACE table.".format(ace) + ) + name, element, Z, mass_number, metastable = get_metadata( + int(zaid), metastable_scheme + ) + + data = cls(name, Z, mass_number, metastable, ace.atomic_weight_ratio) + + # Read energy grid + n_energy = ace.nxs[3] + i = ace.jxs[1] + energy = ace.xss[i : i + n_energy] * EV_PER_MEV + data.energy = energy + + total_xs = ace.xss[ace.jxs[2] : ace.jxs[2] + n_energy] + nonelastic_xs = ace.xss[ace.jxs[3] : ace.jxs[3] + n_energy] + elastic_xs = total_xs-nonelastic_xs + heating_number = ace.xss[ace.jxs[5] : ace.jxs[5] + n_energy] * EV_PER_MEV + + # Create redundant reaction for total (MT=1) + total = PhotonuclearReaction(1) + total.xs = Tabulated1D(energy, total_xs) + total.redundant = True + data.reactions[1] = total + + # Create redundant reaction for nonelastic (MT=3) + nonelastic = PhotonuclearReaction(3) + nonelastic.xs = Tabulated1D(energy, nonelastic_xs) + nonelastic.redundant = True + data.reactions[3] = nonelastic + + # Create redundant reaction for heating (MT=301) + heating = PhotonuclearReaction(301) + heating.xs = Tabulated1D(energy, heating_number*total_xs) + heating.redundant = True + data.reactions[301] = heating + + # Read each reaction + n_reaction = ace.nxs[4] + 1 + for i in range(n_reaction): + if i==0: + if ace.jxs[4]==0: + assert np.allclose(total_xs,nonelastic_xs) + else: + raise NotImplementedError('photonuclear elastic scattering is not supported.') + else: + rx = PhotonuclearReaction.from_ace(ace, i) + data.reactions[rx.mt] = rx + + return data + + def export_to_hdf5(self, path, mode="a", libver="earliest"): + """Export incident photo-nuclear data to an HDF5 file. + + Parameters + ---------- + path : str + Path to write HDF5 file to + mode : {'r', 'r+', 'w', 'x', 'a'} + Mode that is used to open the HDF5 file. This is the second argument + to the :class:`h5py.File` constructor. + libver : {'earliest', 'latest'} + Compatibility mode for the HDF5 file. 'latest' will produce files + that are less backwards compatible but have performance benefits. + + """ + + # Open file and write version + f = h5py.File(str(path), mode, libver=libver) + f.attrs["filetype"] = np.bytes_("data_photonuclear") + if "version" not in f.attrs: + f.attrs["version"] = np.array(HDF5_VERSION) + + # If data come from ENDF, don't allow exporting to HDF5 + if hasattr(self, '_evaluation'): + raise NotImplementedError('Cannot export incident photonuclear data that ' + 'originated from an ENDF file.') + # Write basic data + g = f.create_group(self.name) + g.attrs['Z'] = self.atomic_number + g.attrs['A'] = self.mass_number + g.attrs['metastable'] = self.metastable + g.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio + + # Determine union energy grid + union_grid = np.array([]) + for rx in self: + union_grid = np.union1d(union_grid, rx.xs.x) + for product in rx.products: + union_grid = np.union1d(union_grid, product.yield_.x) + g.create_dataset("energy", data=union_grid) + + # Write reaction data + rxs_group = g.create_group('reactions') + for rx in self.reactions.values(): + # Skip writing redundant reaction if it doesn't have neutron + # production or is a summed transmutation reaction. + # Also write heating. + if rx.redundant: + neutron_rx = any(p.particle == 'neutron' for p in rx.products) + keep_mts = (301,) + if not (neutron_rx or rx.mt in keep_mts): + continue + + rx_group = rxs_group.create_group(f'reaction_{rx.mt:03}') + rx.to_hdf5(rx_group) + + f.close() + + @classmethod + def from_hdf5(cls, group_or_filename): + """Generate photo-nuclear interaction data from HDF5 group + + Parameters + ---------- + group_or_filename : h5py.Group or str + HDF5 group containing interaction data. If given as a string, it is + assumed to be the filename for the HDF5 file, and the first group is + used to read from. + + Returns + ------- + openmc.data.Photonuclear + photo-nuclear interaction data + + """ + if isinstance(group_or_filename, h5py.Group): + group = group_or_filename + else: + h5file = h5py.File(str(group_or_filename), "r") + + # Make sure version matches + if "version" in h5file.attrs: + major, minor = h5file.attrs["version"] + # For now all versions of HDF5 data can be read + else: + raise IOError( + "HDF5 data does not indicate a version. Your installation of " + "the OpenMC Python API expects version {}.x data.".format( + HDF5_VERSION_MAJOR + ) + ) + + group = list(h5file.values())[0] + + name = group.name[1:] + atomic_number = group.attrs["Z"] + mass_number = group.attrs["A"] + metastable = group.attrs["metastable"] + atomic_weight_ratio = group.attrs["atomic_weight_ratio"] + + data = cls(name, atomic_number, mass_number, metastable, atomic_weight_ratio) + + # Read energy grid + data.energy = group["energy"][()] + + # Read reaction data + rxs_group = group["reactions"] + for name, obj in sorted(rxs_group.items()): + if name.startswith("reaction_"): + rx = PhotonuclearReaction.from_hdf5(obj, data.energy) + data.reactions[rx.mt] = rx + + return data + + @classmethod + def from_njoy(cls, filename, evaluation=None, **kwargs): + """Generate incident photo-nuclear data by running NJOY. + + Parameters + ---------- + filename : str + Path to ENDF file + evaluation : openmc.data.endf.Evaluation, optional + If the ENDF file contains multiple material evaluations, this + argument indicates which evaluation to use. + **kwargs + Keyword arguments passed to :func:`openmc.data.njoy.make_ace_photonuclear` + + Returns + ------- + data : openmc.data.Photonuclear + Incident photo-nuclear continuous-energy data + + """ + with tempfile.TemporaryDirectory() as tmpdir: + # Run NJOY to create an ACE library + kwargs.setdefault("output_dir", tmpdir) + for key in ("acer", "pendf"): + kwargs.setdefault(key, os.path.join(kwargs["output_dir"], key)) + kwargs["evaluation"] = evaluation + make_ace_photonuclear(filename, [0.0], **kwargs) + + # Create instance from ACE tables within library + lib = Library(kwargs["acer"]) + data = cls.from_ace(lib.tables[0]) + + return data + + def _get_redundant_reaction(self, mt, mts): + """Create redundant reaction from its components + + Parameters + ---------- + mt : int + MT value of the desired reaction + mts : iterable of int + MT values of its components + + Returns + ------- + openmc.Reaction + Redundant reaction + + """ + + rx = PhotonuclearReaction(mt) + # Get energy grid + energy = self.energy + xss = [self.reactions[mt_i].xs for mt_i in mts] + idx = min([xs._threshold_idx if hasattr(xs, '_threshold_idx') + else 0 for xs in xss]) + rx.xs = Tabulated1D(energy[idx:], Sum(xss)(energy[idx:])) + rx.xs._threshold_idx = idx + + rx.redundant = True + + return rx diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index 65b59582cf8..553bda5323d 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -87,6 +87,8 @@ def _get_products(ev, mt): is not defined in the 'center-of-mass' system. The breakup logic is not implemented which can lead to this error being raised while the definition of the product is correct. + NotImplementedError + When the projectile is not a neutron and not a photon. Returns ------- @@ -163,11 +165,16 @@ def _get_products(ev, mt): ) zat = ev.target["atomic_number"] * 1000 + ev.target["mass_number"] - projectile_mass = ev.projectile["mass"] + if ev.projectile['mass'] == 0.0: + projectile_za = 0 + elif np.isclose(ev.projectile['mass'], 1.0, atol=1.0e-12, rtol=0.): + projectile_za = 1 + else: + raise NotImplementedError('Unknown projectile') p.distribution = [KalbachMann.from_endf(file_obj, za, zat, - projectile_mass)] + projectile_za)] elif law == 2: # Discrete two-body scattering diff --git a/src/cross_sections.cpp b/src/cross_sections.cpp index b1bfde03d13..511628948e6 100644 --- a/src/cross_sections.cpp +++ b/src/cross_sections.cpp @@ -53,6 +53,8 @@ Library::Library(pugi::xml_node node, const std::string& directory) type_ = Type::thermal; } else if (type == "photon") { type_ = Type::photon; + } else if (type == "photonuclear") { + type_ = Type::photonuclear; } else if (type == "wmp") { type_ = Type::wmp; } else { diff --git a/tests/unit_tests/test_data_kalbach_mann.py b/tests/unit_tests/test_data_kalbach_mann.py index 3b837af7da9..df281195541 100644 --- a/tests/unit_tests/test_data_kalbach_mann.py +++ b/tests/unit_tests/test_data_kalbach_mann.py @@ -102,16 +102,14 @@ def test_kalbach_slope(): energy_projectile = 10.2 # [eV] energy_emitted = 5.4 # [eV] - # Check that NotImplementedError is raised if the projectile is not - # a neutron - with pytest.raises(NotImplementedError): - kalbach_slope( - energy_projectile=energy_projectile, - energy_emitted=energy_emitted, - za_projectile=1000, - za_emitted=1, - za_target=6012 - ) + + kalbach_slope( + energy_projectile=energy_projectile, + energy_emitted=energy_emitted, + za_projectile=0, + za_emitted=1, + za_target=6012 + ) assert kalbach_slope( energy_projectile=energy_projectile, diff --git a/tools/ci/download-xs.sh b/tools/ci/download-xs.sh index e831d8215ad..20e0785f6d6 100755 --- a/tools/ci/download-xs.sh +++ b/tools/ci/download-xs.sh @@ -3,11 +3,11 @@ set -ex # Download HDF5 data if [[ ! -e $HOME/nndc_hdf5/cross_sections.xml ]]; then - wget -q -O - https://anl.box.com/shared/static/teaup95cqv8s9nn56hfn7ku8mmelr95p.xz | tar -C $HOME -xJ + wget -q -O - https://github.com/GuySten/data/releases/download/v1.0.0/nndc_hdf5_test.tar.xz | tar -C $HOME -xJ fi # Download ENDF/B-VII.1 distribution ENDF=$HOME/endf-b-vii.1 -if [[ ! -d $ENDF/neutrons || ! -d $ENDF/photoat || ! -d $ENDF/atomic_relax ]]; then +if [[ ! -d $ENDF/neutrons || ! -d $ENDF/photoat || ! -d $ENDF/atomic_relax || ! -d $ENDF/gammas ]]; then wget -q -O - https://anl.box.com/shared/static/4kd2gxnf4gtk4w1c8eua5fsua22kvgjb.xz | tar -C $HOME -xJ fi From 222ec3995a327eb49dce3994dba7b953a545b772 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 11 Jun 2025 19:59:20 +0300 Subject: [PATCH 02/22] added various fixes --- openmc/data/photonuclear.py | 68 ++++++++++++++++++++++++++++++++++--- 1 file changed, 63 insertions(+), 5 deletions(-) diff --git a/openmc/data/photonuclear.py b/openmc/data/photonuclear.py index 810a01c07e5..5e3df98b628 100644 --- a/openmc/data/photonuclear.py +++ b/openmc/data/photonuclear.py @@ -18,7 +18,7 @@ from .angle_energy import AngleEnergy from .correlated import CorrelatedAngleEnergy from .data import ATOMIC_SYMBOL, EV_PER_MEV -from .endf import Evaluation, get_head_record, get_tab1_record, get_cont_record +from .endf import Evaluation, SUM_RULES, get_head_record, get_tab1_record, get_cont_record from .function import Tabulated1D from .njoy import make_ace_photonuclear from .reaction import Reaction, REACTION_NAME, FISSION_MTS, _get_products, _get_fission_products_endf, _get_photon_products_endf, _get_activation_products @@ -118,7 +118,7 @@ def from_endf(cls, ev, mt): Returns ------- - rx : openmc.data.Reaction + rx : openmc.data.PhotonuclearReaction Reaction data """ @@ -353,7 +353,7 @@ def from_ace(cls, ace, i_reaction): idx = int(ace.xss[int(ace.xss[loc+9])+i_mtr]) if idx==0: # No angular distribution data are given for this reaction, - # isotropic scattering is asssumed in LAB + # isotropic scattering is assumed in LAB energy = np.array([particle.yield_.x[0], particle.yield_.x[-1]]) mu_isotropic = Uniform(-1., 1.) distribution.angle = AngleDistribution( @@ -421,6 +421,9 @@ class IncidentPhotonuclear(EqualityMixin): Atomic symbol of the nuclide, e.g., 'Zr' atomic_weight_ratio : float Atomic weight ratio of the target nuclide. + fission_energy : None or openmc.data.FissionEnergyRelease + The energy released by fission, tabulated by component (e.g. prompt + neutrons or beta particles) and dependent on incident neutron energy mass_number : int Number of nucleons in the target nucleus metastable : int @@ -443,6 +446,7 @@ def __init__( self.mass_number = mass_number self.reactions = OrderedDict() self.energy = [] + self._fission_energy = None self.metastable = metastable self.atomic_weight_ratio = atomic_weight_ratio @@ -513,6 +517,16 @@ def atomic_weight_ratio(self, atomic_weight_ratio): cv.check_type("atomic weight ratio", atomic_weight_ratio, Real) cv.check_greater_than("atomic weight ratio", atomic_weight_ratio, 0.0) self._atomic_weight_ratio = atomic_weight_ratio + + @property + def fission_energy(self): + return self._fission_energy + + @fission_energy.setter + def fission_energy(self, fission_energy): + cv.check_type('fission energy release', fission_energy, + FissionEnergyRelease) + self._fission_energy = fission_energy @classmethod def from_endf(cls, ev_or_filename): @@ -551,6 +565,11 @@ def from_endf(cls, ev_or_filename): for mf, mt, nc, mod in ev.reaction_list: if mf == 3: data.reactions[mt] = PhotonuclearReaction.from_endf(ev, mt) + + # Read fission energy release (requires that we already know nu for + # fission) + if (1, 458) in ev.section: + data.fission_energy = FissionEnergyRelease.from_endf(ev, data) data._evaluation = ev return data @@ -696,7 +715,12 @@ def export_to_hdf5(self, path, mode="a", libver="earliest"): rx_group = rxs_group.create_group(f'reaction_{rx.mt:03}') rx.to_hdf5(rx_group) - + + # Write fission energy release data + if self.fission_energy is not None: + fer_group = g.create_group('fission_energy_release') + self.fission_energy.to_hdf5(fer_group) + f.close() @classmethod @@ -753,6 +777,11 @@ def from_hdf5(cls, group_or_filename): rx = PhotonuclearReaction.from_hdf5(obj, data.energy) data.reactions[rx.mt] = rx + # Read fission energy release data + if 'fission_energy_release' in group: + fer_group = group['fission_energy_release'] + data.fission_energy = FissionEnergyRelease.from_hdf5(fer_group) + return data @classmethod @@ -781,14 +810,43 @@ def from_njoy(cls, filename, evaluation=None, **kwargs): for key in ("acer", "pendf"): kwargs.setdefault(key, os.path.join(kwargs["output_dir"], key)) kwargs["evaluation"] = evaluation - make_ace_photonuclear(filename, [0.0], **kwargs) + make_ace_photonuclear(filename, **kwargs) # Create instance from ACE tables within library lib = Library(kwargs["acer"]) data = cls.from_ace(lib.tables[0]) + + # Add fission energy release data + ev = evaluation if evaluation is not None else Evaluation(filename) + if (1, 458) in ev.section: + data.fission_energy = FissionEnergyRelease.from_endf(ev, data) return data + def get_reaction_components(self, mt): + """Determine what reactions make up redundant reaction. + + Parameters + ---------- + mt : int + ENDF MT number of the reaction to find components of. + + Returns + ------- + mts : list of int + ENDF MT numbers of reactions that make up the redundant reaction and + have cross sections provided. + + """ + mts = [] + if mt in SUM_RULES: + for mt_i in SUM_RULES[mt]: + mts += self.get_reaction_components(mt_i) + if mts: + return mts + else: + return [mt] if mt in self else [] + def _get_redundant_reaction(self, mt, mts): """Create redundant reaction from its components From cc8d1a1f6c9fd0ff35d93cf4143da16c26354332 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 11 Jun 2025 20:51:25 +0300 Subject: [PATCH 03/22] added missing import --- openmc/data/photonuclear.py | 1 + 1 file changed, 1 insertion(+) diff --git a/openmc/data/photonuclear.py b/openmc/data/photonuclear.py index 5e3df98b628..c40e5281ac6 100644 --- a/openmc/data/photonuclear.py +++ b/openmc/data/photonuclear.py @@ -19,6 +19,7 @@ from .correlated import CorrelatedAngleEnergy from .data import ATOMIC_SYMBOL, EV_PER_MEV from .endf import Evaluation, SUM_RULES, get_head_record, get_tab1_record, get_cont_record +from .fission_energy import FissionEnergyRelease from .function import Tabulated1D from .njoy import make_ace_photonuclear from .reaction import Reaction, REACTION_NAME, FISSION_MTS, _get_products, _get_fission_products_endf, _get_photon_products_endf, _get_activation_products From bf738c4eb7a09d7a45214928fb558ca4d8e3e617 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 11 Jun 2025 20:21:15 +0300 Subject: [PATCH 04/22] added tests --- tests/unit_tests/test_data_photonuclear.py | 142 +++++++++++++++++++++ 1 file changed, 142 insertions(+) create mode 100644 tests/unit_tests/test_data_photonuclear.py diff --git a/tests/unit_tests/test_data_photonuclear.py b/tests/unit_tests/test_data_photonuclear.py new file mode 100644 index 00000000000..755b061f2f2 --- /dev/null +++ b/tests/unit_tests/test_data_photonuclear.py @@ -0,0 +1,142 @@ +from collections.abc import Mapping, Callable +import os + +import numpy as np +import pandas as pd +import pytest +import openmc.data + +from . import needs_njoy + + +@pytest.fixture(scope='module') +def pu239(): + """Pu239 HDF5 data.""" + directory = os.path.dirname(os.environ['OPENMC_CROSS_SECTIONS']) + filename = os.path.join(directory, 'photonuclear', 'Pu239.h5') + return openmc.data.IncidentPhotonuclear.from_hdf5(filename) + +@pytest.fixture(scope='module') +def u235(): + endf_data = os.environ['OPENMC_ENDF_DATA'] + endf_file = os.path.join(endf_data, 'gammas', 'g-092_U_235.endf') + return openmc.data.IncidentPhotonuclear.from_njoy(endf_file) + +@pytest.fixture(scope='module') +def be9(): + """Be9 ENDF data (contains laboratory angle-energy distribution).""" + endf_data = os.environ['OPENMC_ENDF_DATA'] + filename = os.path.join(endf_data, 'gammas', 'g-004_Be_009.endf') + return openmc.data.IncidentPhotonuclear.from_endf(filename) + + +@pytest.fixture(scope='module') +def h2(): + endf_data = os.environ['OPENMC_ENDF_DATA'] + endf_file = os.path.join(endf_data, 'gammas', 'g-001_H_002.endf') + return openmc.data.IncidentPhotonuclear.from_njoy(endf_file) + + +def test_attributes(pu239): + assert pu239.name == 'Pu239' + assert pu239.mass_number == 239 + assert pu239.metastable == 0 + assert pu239.atomic_symbol == 'Pu' + assert pu239.atomic_weight_ratio == pytest.approx(236.9986) + + +def test_fission_energy(pu239): + fer = pu239.fission_energy + assert isinstance(fer, openmc.data.FissionEnergyRelease) + components = ['betas', 'delayed_neutrons', 'delayed_photons', 'fragments', + 'neutrinos', 'prompt_neutrons', 'prompt_photons', 'recoverable', + 'total', 'q_prompt', 'q_recoverable', 'q_total'] + for c in components: + assert isinstance(getattr(fer, c), Callable) + +@needs_njoy +def test_fission_energy_endf(u235): + fer = u235.fission_energy + assert isinstance(fer, openmc.data.FissionEnergyRelease) + components = ['betas', 'delayed_neutrons', 'delayed_photons', 'fragments', + 'neutrinos', 'prompt_neutrons', 'prompt_photons', 'recoverable', + 'total', 'q_prompt', 'q_recoverable', 'q_total'] + for c in components: + assert isinstance(getattr(fer, c), Callable) + + +def test_energy_grid(pu239): + grid = pu239.energy + assert np.all(np.diff(grid) >= 0.0) + + +def test_reactions(pu239): + assert 18 in pu239.reactions + assert isinstance(pu239.reactions[18], openmc.data.PhotonuclearReaction) + with pytest.raises(KeyError): + pu239.reactions[2] + +def test_fission(pu239): + fission = pu239.reactions[18] + assert not fission.center_of_mass + assert fission.q_value == pytest.approx(197380000.0) + assert fission.mt == 18 + assert len(fission.products) == 1 + prompt = fission.products[0] + assert prompt.particle == 'neutron' + assert prompt.yield_(1.0e-5) == pytest.approx(1.74559) + + +@needs_njoy +def test_kerma(run_in_tmpdir, h2): + assert 301 in h2 + h2.export_to_hdf5("H2.h5") + read_in = openmc.data.IncidentPhotonuclear.from_hdf5("H2.h5") + assert 301 in read_in + assert np.all(read_in[301].xs.y == h2[301].xs.y) + + +@needs_njoy +def test_get_reaction_components(h2): + assert h2.get_reaction_components(1) == [50] + assert h2.get_reaction_components(51) == [] + + +def test_export_to_hdf5(tmpdir, pu239, be9): + filename = str(tmpdir.join('pu239.h5')) + pu239.export_to_hdf5(filename) + assert os.path.exists(filename) + with pytest.raises(NotImplementedError): + be9.export_to_hdf5('be9.h5') + + +@needs_njoy +def test_ace_convert(run_in_tmpdir): + endf_data = os.environ['OPENMC_ENDF_DATA'] + filename = os.path.join(endf_data, 'gammas', 'g-001_H_002.endf') + ace_ascii = 'ace_ascii' + ace_binary = 'ace_binary' + openmc.data.njoy.make_ace_photonuclear(filename, acer=ace_ascii) + + # Convert to binary + openmc.data.ace.ascii_to_binary(ace_ascii, ace_binary) + + # Make sure conversion worked + lib_ascii = openmc.data.ace.Library(ace_ascii) + lib_binary = openmc.data.ace.Library(ace_binary) + for tab_a, tab_b in zip(lib_ascii.tables, lib_binary.tables): + assert tab_a.name == tab_b.name + assert tab_a.atomic_weight_ratio == pytest.approx(tab_b.atomic_weight_ratio) + assert tab_a.temperature == pytest.approx(tab_b.temperature) + assert np.all(tab_a.nxs == tab_b.nxs) + assert np.all(tab_a.jxs == tab_b.jxs) + assert tab_a.zaid == tab_b.zaid + assert tab_a.data_type == tab_b.data_type + + +def test_ace_table_types(): + TT = openmc.data.ace.TableType + assert TT.from_suffix('u') == TT.PHOTONUCLEAR + assert TT.from_suffix('80u') == TT.PHOTONUCLEAR + + From 0778b17e811ab00e5b33c118a6ee07433ec72da6 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 11 Jun 2025 22:57:59 +0300 Subject: [PATCH 05/22] applied makeclean's suggestions --- openmc/data/data.py | 3 ++- openmc/data/kalbach_mann.py | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/openmc/data/data.py b/openmc/data/data.py index cd5bc19f238..60d06e4e70c 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -274,6 +274,7 @@ # Unit conversions EV_PER_MEV = 1.0e6 JOULE_PER_EV = 1.602176634e-19 +EV_PER_AMU = 931.49410372e6 # eV/c^2 per amu # Avogadro's constant AVOGADRO = 6.02214076e23 @@ -282,7 +283,7 @@ NEUTRON_MASS = 1.00866491595 # Neutron mass in units of eV/c^2 -NEUTRON_MASS_EV = 939.56542052e6 +NEUTRON_MASS_EV = NEUTRON_MASS*EV_PER_AMU # Used in atomic_mass function as a cache _ATOMIC_MASS: dict[str, float] = {} diff --git a/openmc/data/kalbach_mann.py b/openmc/data/kalbach_mann.py index d5cc3f37f99..9eefc80c6e7 100644 --- a/openmc/data/kalbach_mann.py +++ b/openmc/data/kalbach_mann.py @@ -210,7 +210,7 @@ def kalbach_slope(energy_projectile, energy_emitted, za_projectile, # or ENDF-6 Formats Manual section 6.2.3.2 slope_n = kalbach_slope(energy_projectile, energy_emitted, 1, za_emitted, za_target) - return slope_n * np.sqrt(0.5*energy_projectile/NEUTRON_MASS_EV)*np.minimum(4,np.maximum(1,9.3/np.sqrt(energy_emitted*1e-6))) + return slope_n * np.sqrt(0.5*energy_projectile/NEUTRON_MASS_EV)*np.minimum(4,np.maximum(1,9.3/np.sqrt(energy_emitted/EV_PER_MEV))) # Special handling of elemental carbon if za_emitted == 6000: From faf005fd91331c2ef3d12a60fc491711a070270e Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 11 Jun 2025 23:17:32 +0300 Subject: [PATCH 06/22] add fission energy release to nuclear data rst docs --- docs/source/io_formats/nuclear_data.rst | 28 +++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/docs/source/io_formats/nuclear_data.rst b/docs/source/io_formats/nuclear_data.rst index b376abc4438..3e4ed49ae85 100644 --- a/docs/source/io_formats/nuclear_data.rst +++ b/docs/source/io_formats/nuclear_data.rst @@ -304,6 +304,34 @@ Incident Photonuclear Data **//reactions/reaction_/product_/** Reaction product data is described in :ref:`product`. + +**//fission_energy_release/** + +:Datasets: - **fragments** (:ref:`function <1d_functions>`) -- Energy + released in the form of fragments as a function of incident + neutron energy. + - **prompt_neutrons** (:ref:`function <1d_functions>`) -- Energy + released in the form of prompt neutrons as a function of incident + neutron energy. + - **delayed_neutrons** (:ref:`function <1d_functions>`) -- Energy + released in the form of delayed neutrons as a function of incident + neutron energy. + - **prompt_photons** (:ref:`function <1d_functions>`) -- Energy + released in the form of prompt photons as a function of incident + neutron energy. + - **delayed_photons** (:ref:`function <1d_functions>`) -- Energy + released in the form of delayed photons as a function of incident + neutron energy. + - **betas** (:ref:`function <1d_functions>`) -- Energy released in + the form of betas as a function of incident neutron energy. + - **neutrinos** (:ref:`function <1d_functions>`) -- Energy released + in the form of neutrinos as a function of incident neutron energy. + - **q_prompt** (:ref:`function <1d_functions>`) -- The prompt fission + Q-value (fragments + prompt neutrons + prompt photons - incident + energy) + - **q_recoverable** (:ref:`function <1d_functions>`) -- The + recoverable fission Q-value (Q_prompt + delayed neutrons + delayed + photons + betas) .. _product: From ad62d0884fe058c0774260ab5b6bd49c001769a1 Mon Sep 17 00:00:00 2001 From: GuySten <62616591+GuySten@users.noreply.github.com> Date: Thu, 12 Jun 2025 18:44:16 +0300 Subject: [PATCH 07/22] Update ci.yml --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 38a49b60219..0f9b980df8b 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -153,7 +153,7 @@ jobs: path: | ~/nndc_hdf5 ~/endf-b-vii.1 - key: ${{ runner.os }}-build-xs-cache + key: ${{ runner.os }}-build-xs-cache-${{ hashFiles("$GITHUB_WORKSPACE/tools/ci/download-xs.sh")}} - name: before shell: bash From 57ee46928bed55e0c5cc78685a9dfe62c1bcf7ad Mon Sep 17 00:00:00 2001 From: GuySten <62616591+GuySten@users.noreply.github.com> Date: Thu, 12 Jun 2025 18:57:29 +0300 Subject: [PATCH 08/22] Update ci.yml --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 0f9b980df8b..fd7745ec8d1 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -153,7 +153,7 @@ jobs: path: | ~/nndc_hdf5 ~/endf-b-vii.1 - key: ${{ runner.os }}-build-xs-cache-${{ hashFiles("$GITHUB_WORKSPACE/tools/ci/download-xs.sh")}} + key: ${{ runner.os }}-build-xs-cache-${{ hashFiles(format('{0}/tools/ci/download-xs.sh', github.workspace)) }} - name: before shell: bash From 3b5599733a528ae800e180ced2de924db97b5c04 Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 12 Jun 2025 21:16:35 +0300 Subject: [PATCH 09/22] now do not fail if diff is only in line endings --- tests/regression_tests/model_xml/test.py | 15 ++++++++------ tests/regression_tests/ncrystal/test.py | 11 ++++++---- tests/testing_harness.py | 26 +++++++++++++++--------- tests/unit_tests/mesh_to_vtk/test.py | 6 ++++-- 4 files changed, 36 insertions(+), 22 deletions(-) diff --git a/tests/regression_tests/model_xml/test.py b/tests/regression_tests/model_xml/test.py index c67a72ed372..7716102265e 100644 --- a/tests/regression_tests/model_xml/test.py +++ b/tests/regression_tests/model_xml/test.py @@ -38,11 +38,14 @@ def _compare_results(self): if not compare: expected = open(self.results_true).readlines() actual = open('results_test.dat').readlines() - diff = unified_diff(expected, actual, self.results_true, - 'results_test.dat') - print('Result differences:') - print(''.join(colorize(diff))) - os.rename('results_test.dat', 'results_error.dat') + diff = list(unified_diff(expected, actual, self.results_true, + 'results_test.dat')) + if diff: + print('Result differences:') + print(''.join(colorize(diff))) + os.rename('results_test.dat', 'results_error.dat') + else: + compare = True assert compare, 'Results do not agree' def _cleanup(self): @@ -97,4 +100,4 @@ def test_input_arg(run_in_tmpdir): # even in the presence of other, valid XML files pincell.export_to_model_xml() with pytest.raises(RuntimeError, match='ex-em-ell.xml'): - openmc.run(path_input='ex-em-ell.xml') \ No newline at end of file + openmc.run(path_input='ex-em-ell.xml') diff --git a/tests/regression_tests/ncrystal/test.py b/tests/regression_tests/ncrystal/test.py index 8da05e1cfd7..22683cffef2 100644 --- a/tests/regression_tests/ncrystal/test.py +++ b/tests/regression_tests/ncrystal/test.py @@ -99,8 +99,11 @@ def test_cfg_from_xml(): actual = open('materials.xml.after', 'r').readlines() compare = filecmp.cmp('materials.xml.orig','materials.xml.after') if not compare: - diff = unified_diff(expected, actual, 'materials.xml.orig', - 'materials.xml.after') - print('Input differences:') - print(''.join(diff)) + diff = list(unified_diff(expected, actual, 'materials.xml.orig', + 'materials.xml.after')) + if diff: + print('Input differences:') + print(''.join(diff)) + else: + compare = True assert compare, 'Materials not read correctly from XML' diff --git a/tests/testing_harness.py b/tests/testing_harness.py index 6de475e6e04..1162e727fd4 100644 --- a/tests/testing_harness.py +++ b/tests/testing_harness.py @@ -133,11 +133,14 @@ def _compare_results(self): if not compare: expected = open('results_true.dat').readlines() actual = open('results_test.dat').readlines() - diff = unified_diff(expected, actual, 'results_true.dat', - 'results_test.dat') - print('Result differences:') - print(''.join(colorize(diff))) - os.rename('results_test.dat', 'results_error.dat') + diff = list(unified_diff(expected, actual, 'results_true.dat', + 'results_test.dat')) + if diff: + print('Result differences:') + print(''.join(colorize(diff))) + os.rename('results_test.dat', 'results_error.dat') + else: + compare = True assert compare, 'Results do not agree' def _cleanup(self): @@ -360,11 +363,14 @@ def _compare_inputs(self): if not compare: expected = open(self.inputs_true, 'r').readlines() actual = open('inputs_test.dat', 'r').readlines() - diff = unified_diff(expected, actual, self.inputs_true, - 'inputs_test.dat') - print('Input differences:') - print(''.join(colorize(diff))) - os.rename('inputs_test.dat', 'inputs_error.dat') + diff = list(unified_diff(expected, actual, self.inputs_true, + 'inputs_test.dat')) + if diff: + print('Input differences:') + print(''.join(colorize(diff))) + os.rename('inputs_test.dat', 'inputs_error.dat') + else: + compare = True assert compare, 'Input files are broken.' def _cleanup(self): diff --git a/tests/unit_tests/mesh_to_vtk/test.py b/tests/unit_tests/mesh_to_vtk/test.py index f00aa462612..eadba3979f6 100644 --- a/tests/unit_tests/mesh_to_vtk/test.py +++ b/tests/unit_tests/mesh_to_vtk/test.py @@ -82,7 +82,8 @@ def test_mesh_write_vtk(mesh_params, run_in_tmpdir): assert filecmp.cmp(test_data, filename) except AssertionError as e: diff = diff_file(test_data, filename) - raise AssertionError(diff) from e + if diff: + raise AssertionError(diff) from e # check data writing def test_mesh_write_vtk_data(run_in_tmpdir): @@ -97,7 +98,8 @@ def test_mesh_write_vtk_data(run_in_tmpdir): assert filecmp.cmp(filename, filename_expected) except AssertionError as e: diff = diff_file(filename_expected, filename) - raise AssertionError(diff) from e + if diff: + raise AssertionError(diff) from e From 057fc3df57829993a5bf1ea3e5d210023217fd78 Mon Sep 17 00:00:00 2001 From: GuySten <62616591+GuySten@users.noreply.github.com> Date: Thu, 12 Jun 2025 23:44:03 +0300 Subject: [PATCH 10/22] Update download-xs.sh --- tools/ci/download-xs.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/ci/download-xs.sh b/tools/ci/download-xs.sh index 20e0785f6d6..0c534c0e6e0 100755 --- a/tools/ci/download-xs.sh +++ b/tools/ci/download-xs.sh @@ -9,5 +9,5 @@ fi # Download ENDF/B-VII.1 distribution ENDF=$HOME/endf-b-vii.1 if [[ ! -d $ENDF/neutrons || ! -d $ENDF/photoat || ! -d $ENDF/atomic_relax || ! -d $ENDF/gammas ]]; then - wget -q -O - https://anl.box.com/shared/static/4kd2gxnf4gtk4w1c8eua5fsua22kvgjb.xz | tar -C $HOME -xJ + wget -q -O - https://github.com/GuySten/data/releases/download/v1.0.0/endf71.tar.xz | tar -C $HOME -xJ fi From 049ea70eb9931a69844c55cb52e398dc360cd293 Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 13 Jun 2025 01:00:42 +0300 Subject: [PATCH 11/22] reverted path --- tools/ci/download-xs.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/ci/download-xs.sh b/tools/ci/download-xs.sh index 0c534c0e6e0..20e0785f6d6 100755 --- a/tools/ci/download-xs.sh +++ b/tools/ci/download-xs.sh @@ -9,5 +9,5 @@ fi # Download ENDF/B-VII.1 distribution ENDF=$HOME/endf-b-vii.1 if [[ ! -d $ENDF/neutrons || ! -d $ENDF/photoat || ! -d $ENDF/atomic_relax || ! -d $ENDF/gammas ]]; then - wget -q -O - https://github.com/GuySten/data/releases/download/v1.0.0/endf71.tar.xz | tar -C $HOME -xJ + wget -q -O - https://anl.box.com/shared/static/4kd2gxnf4gtk4w1c8eua5fsua22kvgjb.xz | tar -C $HOME -xJ fi From 09dae65fd42176f734065ff8c49b03ce886cd984 Mon Sep 17 00:00:00 2001 From: GuySten <62616591+GuySten@users.noreply.github.com> Date: Sat, 14 Jun 2025 14:25:55 +0300 Subject: [PATCH 12/22] Apply suggestions from code review Co-authored-by: Ahnaf Tahmid Chowdhury --- openmc/data/photonuclear.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/openmc/data/photonuclear.py b/openmc/data/photonuclear.py index c40e5281ac6..6422a4c6f28 100644 --- a/openmc/data/photonuclear.py +++ b/openmc/data/photonuclear.py @@ -329,6 +329,8 @@ def from_ace(cls, ace, i_reaction): case _: #TODO: support more product particles # for now skip particle if it is not a neutron or photon + warn(f"Unsupported secondary particle type in {ace.name} for MT={mt}. " + "This product will be skipped.") continue # Determine reference frame @@ -759,6 +761,7 @@ def from_hdf5(cls, group_or_filename): ) group = list(h5file.values())[0] + h5file.close() name = group.name[1:] atomic_number = group.attrs["Z"] From 7a1fe592bc158cd51c0352df777b79cfb36cafb2 Mon Sep 17 00:00:00 2001 From: GuySten <62616591+GuySten@users.noreply.github.com> Date: Sat, 14 Jun 2025 14:41:12 +0300 Subject: [PATCH 13/22] Update IncidentPhotonuclear to support sum rules like IncidentNeutron --- openmc/data/photonuclear.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/openmc/data/photonuclear.py b/openmc/data/photonuclear.py index 6422a4c6f28..4acc2d6c025 100644 --- a/openmc/data/photonuclear.py +++ b/openmc/data/photonuclear.py @@ -460,8 +460,13 @@ def __getitem__(self, mt): if mt in self.reactions: return self.reactions[mt] else: - raise KeyError("No reaction with MT={}.".format(mt)) - + # Try to create a redundant cross section + mts = self.get_reaction_components(mt) + if len(mts) > 0: + return self._get_redundant_reaction(mt, mts) + else: + raise KeyError(f'No reaction with MT={mt}.') + def __repr__(self): return "".format(self.name) From c0c6802919b54be7aa45802d20eb62bb41c0ba87 Mon Sep 17 00:00:00 2001 From: GuySten <62616591+GuySten@users.noreply.github.com> Date: Sat, 14 Jun 2025 14:43:53 +0300 Subject: [PATCH 14/22] Update photonuclear.py --- openmc/data/photonuclear.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/openmc/data/photonuclear.py b/openmc/data/photonuclear.py index 4acc2d6c025..12045f0555e 100644 --- a/openmc/data/photonuclear.py +++ b/openmc/data/photonuclear.py @@ -565,7 +565,10 @@ def from_endf(cls, ev_or_filename): atomic_weight_ratio = ev.target["mass"] element = ATOMIC_SYMBOL[atomic_number] - name = "{}{}".format(element, mass_number) + if metastable > 0: + name = f'{element}{mass_number}_m{metastable}' + else: + name = f'{element}{mass_number}' data = cls(name, atomic_number, mass_number, metastable, atomic_weight_ratio) From 497a1fc1a59c6005d25e19a8e04a7e942c0a7a03 Mon Sep 17 00:00:00 2001 From: GuySten <62616591+GuySten@users.noreply.github.com> Date: Sat, 14 Jun 2025 22:30:23 +0300 Subject: [PATCH 15/22] close h5file when done --- openmc/data/photonuclear.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/openmc/data/photonuclear.py b/openmc/data/photonuclear.py index 12045f0555e..d3fb0328976 100644 --- a/openmc/data/photonuclear.py +++ b/openmc/data/photonuclear.py @@ -769,7 +769,6 @@ def from_hdf5(cls, group_or_filename): ) group = list(h5file.values())[0] - h5file.close() name = group.name[1:] atomic_number = group.attrs["Z"] @@ -794,6 +793,10 @@ def from_hdf5(cls, group_or_filename): fer_group = group['fission_energy_release'] data.fission_energy = FissionEnergyRelease.from_hdf5(fer_group) + # Close h5file when done if was opened + if not isinstance(group_or_filename, h5py.Group): + h5file.close() + return data @classmethod From 060ca6563b63f4623eb51f9436b4c05c0acf87c0 Mon Sep 17 00:00:00 2001 From: GuySten <62616591+GuySten@users.noreply.github.com> Date: Sat, 14 Jun 2025 22:42:35 +0300 Subject: [PATCH 16/22] fix test --- tests/unit_tests/test_data_photonuclear.py | 14 ++++---------- 1 file changed, 4 insertions(+), 10 deletions(-) diff --git a/tests/unit_tests/test_data_photonuclear.py b/tests/unit_tests/test_data_photonuclear.py index 755b061f2f2..eda507140f6 100644 --- a/tests/unit_tests/test_data_photonuclear.py +++ b/tests/unit_tests/test_data_photonuclear.py @@ -16,12 +16,14 @@ def pu239(): filename = os.path.join(directory, 'photonuclear', 'Pu239.h5') return openmc.data.IncidentPhotonuclear.from_hdf5(filename) + @pytest.fixture(scope='module') def u235(): endf_data = os.environ['OPENMC_ENDF_DATA'] endf_file = os.path.join(endf_data, 'gammas', 'g-092_U_235.endf') return openmc.data.IncidentPhotonuclear.from_njoy(endf_file) + @pytest.fixture(scope='module') def be9(): """Be9 ENDF data (contains laboratory angle-energy distribution).""" @@ -45,17 +47,8 @@ def test_attributes(pu239): assert pu239.atomic_weight_ratio == pytest.approx(236.9986) -def test_fission_energy(pu239): - fer = pu239.fission_energy - assert isinstance(fer, openmc.data.FissionEnergyRelease) - components = ['betas', 'delayed_neutrons', 'delayed_photons', 'fragments', - 'neutrinos', 'prompt_neutrons', 'prompt_photons', 'recoverable', - 'total', 'q_prompt', 'q_recoverable', 'q_total'] - for c in components: - assert isinstance(getattr(fer, c), Callable) - @needs_njoy -def test_fission_energy_endf(u235): +def test_fission_energy(u235): fer = u235.fission_energy assert isinstance(fer, openmc.data.FissionEnergyRelease) components = ['betas', 'delayed_neutrons', 'delayed_photons', 'fragments', @@ -76,6 +69,7 @@ def test_reactions(pu239): with pytest.raises(KeyError): pu239.reactions[2] + def test_fission(pu239): fission = pu239.reactions[18] assert not fission.center_of_mass From 54ecab6a36e99ee8d1bda58f9037ca3075272a6c Mon Sep 17 00:00:00 2001 From: GuySten <62616591+GuySten@users.noreply.github.com> Date: Sun, 15 Jun 2025 00:13:17 +0300 Subject: [PATCH 17/22] Ignore photonuclear data when expanding elements --- openmc/element.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/element.py b/openmc/element.py index f1d8f5f24d5..c0fb978e9a3 100644 --- a/openmc/element.py +++ b/openmc/element.py @@ -142,7 +142,7 @@ def expand(self, percent, percent_type, enrichment=None, library_nuclides = set() tree = ET.parse(cross_sections) root = tree.getroot() - for child in root.findall('library'): + for child in root.findall("library[not(@type = 'photonuclear')]"): nuclide = child.attrib['materials'] if re.match(r'{}\d+'.format(self), nuclide): library_nuclides.add(nuclide) From b802e4aeef4908af58d7368653aab17395305f82 Mon Sep 17 00:00:00 2001 From: GuySten <62616591+GuySten@users.noreply.github.com> Date: Sun, 15 Jun 2025 00:34:03 +0300 Subject: [PATCH 18/22] fixed bug --- openmc/element.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/openmc/element.py b/openmc/element.py index c0fb978e9a3..855c453064b 100644 --- a/openmc/element.py +++ b/openmc/element.py @@ -142,7 +142,9 @@ def expand(self, percent, percent_type, enrichment=None, library_nuclides = set() tree = ET.parse(cross_sections) root = tree.getroot() - for child in root.findall("library[not(@type = 'photonuclear')]"): + for child in root.findall('library'): + if child.attrib['type'] == 'photonuclear': + continue nuclide = child.attrib['materials'] if re.match(r'{}\d+'.format(self), nuclide): library_nuclides.add(nuclide) From 0e6e5c95c9d1bfae2aab8e819b916853c5df8fa4 Mon Sep 17 00:00:00 2001 From: GuySten Date: Mon, 21 Jul 2025 11:45:51 +0300 Subject: [PATCH 19/22] added mt=50 to REACTION_NAME dict --- openmc/data/photonuclear.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/openmc/data/photonuclear.py b/openmc/data/photonuclear.py index d3fb0328976..77ac5032490 100644 --- a/openmc/data/photonuclear.py +++ b/openmc/data/photonuclear.py @@ -22,7 +22,7 @@ from .fission_energy import FissionEnergyRelease from .function import Tabulated1D from .njoy import make_ace_photonuclear -from .reaction import Reaction, REACTION_NAME, FISSION_MTS, _get_products, _get_fission_products_endf, _get_photon_products_endf, _get_activation_products +from .reaction import Reaction, REACTION_NAME as _REACTION_NAME, FISSION_MTS, _get_products, _get_fission_products_endf, _get_photon_products_endf, _get_activation_products from .product import Product from .energy_distribution import EnergyDistribution, LevelInelastic, \ DiscretePhoton @@ -33,7 +33,8 @@ from .product import Product from .uncorrelated import UncorrelatedAngleEnergy -_REACTION_NAME = {key:value.replace("(n,","(gamma,") for key,value in REACTION_NAME.items()} +REACTION_NAME = {50 : '(gamma,n0)'} +REACTION_NAME.update({key:value.replace("(n,","(gamma,") for key,value in REACTION_NAME_.items()}) class PhotonuclearReaction(EqualityMixin): def __init__(self, mt): @@ -47,7 +48,7 @@ def __init__(self, mt): def __repr__(self): if self.mt in _REACTION_NAME: - return f"" + return f"" else: return f"" @@ -211,8 +212,8 @@ def to_hdf5(self, group): """ group.attrs['mt'] = self.mt - if self.mt in _REACTION_NAME: - group.attrs['label'] = np.bytes_(_REACTION_NAME[self.mt]) + if self.mt in REACTION_NAME: + group.attrs['label'] = np.bytes_(REACTION_NAME[self.mt]) else: group.attrs['label'] = np.bytes_(self.mt) group.attrs['Q_value'] = self.q_value From c13c6a21b911026a302c16b3174ca13a198e31c3 Mon Sep 17 00:00:00 2001 From: GuySten Date: Mon, 21 Jul 2025 12:30:59 +0300 Subject: [PATCH 20/22] fixed union grid logic --- openmc/data/photonuclear.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/openmc/data/photonuclear.py b/openmc/data/photonuclear.py index 77ac5032490..af329ecf328 100644 --- a/openmc/data/photonuclear.py +++ b/openmc/data/photonuclear.py @@ -34,7 +34,7 @@ from .uncorrelated import UncorrelatedAngleEnergy REACTION_NAME = {50 : '(gamma,n0)'} -REACTION_NAME.update({key:value.replace("(n,","(gamma,") for key,value in REACTION_NAME_.items()}) +REACTION_NAME.update({key:value.replace("(n,","(gamma,") for key,value in _REACTION_NAME.items()}) class PhotonuclearReaction(EqualityMixin): def __init__(self, mt): @@ -713,6 +713,17 @@ def export_to_hdf5(self, path, mode="a", libver="earliest"): union_grid = np.union1d(union_grid, product.yield_.x) g.create_dataset("energy", data=union_grid) + # Update reactions according to union energy grid + for rx in self: + if tuple(rx.xs.interpolation) != (2,): + raise NotImplementedError('Only linear-linear interpolable reactions are supported.') + + threshold_idx, = np.flatnonzero(union_grid == rx.xs.x[0]) + xs = rx.xs(union_grid[threshold_idx:]) + tab1d = Tabulated1D(union_grid[threshold_idx:], xs) + tab1d._threshold_idx = threshold_idx + rx.xs = tab1d + # Write reaction data rxs_group = g.create_group('reactions') for rx in self.reactions.values(): From b55a61fd2105f32c2974a9eb7e096b614ea884c0 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 23 Jul 2025 21:30:06 +0300 Subject: [PATCH 21/22] invalidate test data xs cache --- tools/ci/download-xs.sh | 1 + 1 file changed, 1 insertion(+) diff --git a/tools/ci/download-xs.sh b/tools/ci/download-xs.sh index 20e0785f6d6..dd7775e9e6e 100755 --- a/tools/ci/download-xs.sh +++ b/tools/ci/download-xs.sh @@ -1,4 +1,5 @@ #!/bin/bash +# touching this script invalidate the test data xs cache set -ex # Download HDF5 data From f5682d62df7564621a59dc19a9bbb101550b2ce1 Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 29 Aug 2025 15:46:48 +0300 Subject: [PATCH 22/22] simplified code --- tests/unit_tests/test_data_photonuclear.py | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/tests/unit_tests/test_data_photonuclear.py b/tests/unit_tests/test_data_photonuclear.py index eda507140f6..abb318d8ad4 100644 --- a/tests/unit_tests/test_data_photonuclear.py +++ b/tests/unit_tests/test_data_photonuclear.py @@ -5,6 +5,7 @@ import pandas as pd import pytest import openmc.data +import openmc from . import needs_njoy @@ -12,29 +13,26 @@ @pytest.fixture(scope='module') def pu239(): """Pu239 HDF5 data.""" - directory = os.path.dirname(os.environ['OPENMC_CROSS_SECTIONS']) + directory = os.path.dirname(openmc.config.get('cross_sections')) filename = os.path.join(directory, 'photonuclear', 'Pu239.h5') return openmc.data.IncidentPhotonuclear.from_hdf5(filename) @pytest.fixture(scope='module') -def u235(): - endf_data = os.environ['OPENMC_ENDF_DATA'] +def u235(endf_data): endf_file = os.path.join(endf_data, 'gammas', 'g-092_U_235.endf') return openmc.data.IncidentPhotonuclear.from_njoy(endf_file) @pytest.fixture(scope='module') -def be9(): +def be9(endf_data): """Be9 ENDF data (contains laboratory angle-energy distribution).""" - endf_data = os.environ['OPENMC_ENDF_DATA'] filename = os.path.join(endf_data, 'gammas', 'g-004_Be_009.endf') return openmc.data.IncidentPhotonuclear.from_endf(filename) @pytest.fixture(scope='module') -def h2(): - endf_data = os.environ['OPENMC_ENDF_DATA'] +def h2(endf_data): endf_file = os.path.join(endf_data, 'gammas', 'g-001_H_002.endf') return openmc.data.IncidentPhotonuclear.from_njoy(endf_file) @@ -105,8 +103,7 @@ def test_export_to_hdf5(tmpdir, pu239, be9): @needs_njoy -def test_ace_convert(run_in_tmpdir): - endf_data = os.environ['OPENMC_ENDF_DATA'] +def test_ace_convert(endf_data, run_in_tmpdir): filename = os.path.join(endf_data, 'gammas', 'g-001_H_002.endf') ace_ascii = 'ace_ascii' ace_binary = 'ace_binary'