diff --git a/docs/source/pythonapi/data.rst b/docs/source/pythonapi/data.rst index 454409f1f..c7c65e14a 100644 --- a/docs/source/pythonapi/data.rst +++ b/docs/source/pythonapi/data.rst @@ -67,6 +67,7 @@ Core Functions gnd_name half_life isotopes + kalbach_slope linearize thin water_density diff --git a/openmc/data/kalbach_mann.py b/openmc/data/kalbach_mann.py index f4c914d90..a036da694 100644 --- a/openmc/data/kalbach_mann.py +++ b/openmc/data/kalbach_mann.py @@ -5,6 +5,7 @@ from warnings import warn import numpy as np import openmc.checkvalue as cv +from openmc.mixin import EqualityMixin from openmc.stats import Tabular, Univariate, Discrete, Mixture from .function import Tabulated1D, INTERPOLATION_SCHEME from .angle_energy import AngleEnergy @@ -12,6 +13,240 @@ from .data import EV_PER_MEV from .endf import get_list_record, get_tab2_record +class _AtomicRepresentation(EqualityMixin): + """Atomic representation of an isotope or a particle. + + Parameters + ---------- + z : int + Number of protons (atomic number) + a : int + Number of nucleons (mass number) + + Raises + ------ + ValueError + When the number of protons (z) declared is higher than the number + of nucleons (a) + + Attributes + ---------- + z : int + Number of protons (atomic number) + a : int + Number of nucleons (mass number) + n : int + Number of neutrons + za : int + ZA identifier, 1000*Z + A, where Z is the atomic number and A the mass + number + + """ + def __init__(self, z, a): + # Sanity checks on values + cv.check_type('z', z, Integral) + cv.check_greater_than('z', z, 0, equality=True) + cv.check_type('a', a, Integral) + cv.check_greater_than('a', a, 0, equality=True) + if z > a: + raise ValueError(f"Number of protons ({z}) must be less than or " + f"equal to number of nucleons ({a}).") + + self._z = z + self._a = a + + def __add__(self, other): + """Add two _AtomicRepresentations""" + z = self.z + other.z + a = self.a + other.a + return _AtomicRepresentation(z=z, a=a) + + def __sub__(self, other): + """Substract two _AtomicRepresentations""" + z = self.z - other.z + a = self.a - other.a + return _AtomicRepresentation(z=z, a=a) + + @property + def a(self): + return self._a + + @property + def z(self): + return self._z + + @property + def n(self): + return self.a - self.z + + @property + def za(self): + return self.z * 1000 + self.a + + @classmethod + def from_za(cls, za): + """Instantiate an _AtomicRepresentation from a ZA identifier. + + Parameters + ---------- + za : int + ZA identifier, 1000*Z + A, where Z is the atomic number and A the + mass number + + Returns + ------- + _AtomicRepresentation + Atomic representation of the isotope/particle + + """ + z, a = divmod(za, 1000) + return cls(z, a) + + +def _separation_energy(compound, nucleus, particle): + """Calculates the separation energy as defined in ENDF-6 manual + BNL-203218-2018-INRE, Revision 215, File 6 description for LAW=1 + and LANG=2. This function can be used for the incident or emitted + particle of the following reaction: A + a -> C -> B + b + + Parameters + ---------- + compound : _AtomicRepresentation + Atomic representation of the compound (C) + nucleus : _AtomicRepresentation + Atomic representation of the nucleus (A or B) + particle : _AtomicRepresentation + Atomic representation of the particle (a or b) + + Returns + ------- + separation_energy : float + Separation energy in MeV + + """ + # Determine A, Z, and N for compound and nucleus + A_c = compound.a + Z_c = compound.z + N_c = compound.n + A_a = nucleus.a + Z_a = nucleus.z + N_a = nucleus.n + + # Determine breakup energy of incident particle (ENDF-6 Formats Manual, + # Appendix H, Table 3) in MeV + za_to_breaking_energy = { + 1: 0.0, + 1001: 0.0, + 1002: 2.224566, + 1003: 8.481798, + 2003: 7.718043, + 2004: 28.29566 + } + I_a = za_to_breaking_energy[particle.za] + + # Eq. 4 in in doi:10.1103/PhysRevC.37.2350 or ENDF-6 Formats Manual section + # 6.2.3.2 + return ( + 15.68 * (A_c - A_a) - + 28.07 * ((N_c - Z_c)**2 / A_c - (N_a - Z_a)**2 / A_a) - + 18.56 * (A_c**(2./3.) - A_a**(2./3.)) + + 33.22 * ((N_c - Z_c)**2 / A_c**(4./3.) - (N_a - Z_a)**2 / A_a**(4./3.)) - + 0.717 * (Z_c**2 / A_c**(1./3.) - Z_a**2 / A_a**(1./3.)) + + 1.211 * (Z_c**2 / A_c - Z_a**2 / A_a) - + I_a + ) + + +def kalbach_slope(energy_projectile, energy_emitted, za_projectile, + za_emitted, za_target): + """Returns Kalbach-Mann slope from calculations. + + The associated reaction is defined as: + A + a -> C -> B + b + + Where: + + - A is the targeted nucleus, + - a is the projectile, + - C is the compound, + - B is the residual nucleus, + - b is the emitted particle. + + The Kalbach-Mann slope calculation is done as defined in ENDF-6 manual + BNL-203218-2018-INRE, Revision 215, File 6 description for LAW=1 and + LANG=2. One exception to this, is that the entrance and emission channel + energies are not calculated with the AWR number, but approximated with + the number of mass instead. + + Parameters + ---------- + energy_projectile : float + Energy of the projectile in the laboratory system in eV + energy_emitted : float + Energy of the emitted particle in the center of mass system in eV + za_projectile : int + ZA identifier of the projectile + za_emitted : int + ZA identifier of the emitted particle + za_target : int + ZA identifier of the targeted nucleus + + Raises + ------ + NotImplementedError + When the projectile is not a neutron + + Returns + ------- + slope : float + 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." + ) + + # Special handling of elemental carbon + if za_emitted == 6000: + za_emitted = 6012 + if za_target == 6000: + za_target = 6012 + + projectile = _AtomicRepresentation.from_za(za_projectile) + emitted = _AtomicRepresentation.from_za(za_emitted) + target = _AtomicRepresentation.from_za(za_target) + compound = projectile + target + residual = compound - emitted + + # Calculate entrance and emission channel energy in MeV, defined in section + # 6.2.3.2 in the ENDF-6 Formats Manual + epsilon_a = energy_projectile * target.a / (target.a + projectile.a) / EV_PER_MEV + epsilon_b = energy_emitted * (residual.a + emitted.a) \ + / (residual.a * EV_PER_MEV) + + # Calculate separation energies using Eq. 4 in doi:10.1103/PhysRevC.37.2350 + # or ENDF-6 Formats Manual section 6.2.3.2 + s_a = _separation_energy(compound, target, projectile) + s_b = _separation_energy(compound, residual, emitted) + + # See Eq. 10 in doi:10.1103/PhysRevC.37.2350 or section 6.2.3.2 in the + # ENDF-6 Formats Manual + za_to_M = {1: 1.0, 1001: 1.0, 1002: 1.0, 2004: 0.0} + za_to_m = {1: 0.5, 1001: 1.0, 1002: 1.0, 1003: 1.0, 2003: 1.0, 2004: 2.0} + M = za_to_M[projectile.za] + m = za_to_m[emitted.za] + e_a = epsilon_a + s_a + e_b = epsilon_b + s_b + r_1 = min(e_a, 130.) + r_3 = min(e_a, 41.) + x_1 = r_1 * e_b / e_a + x_3 = r_3 * e_b / e_a + return 0.04 * x_1 + 1.8e-6 * x_1**3 + 6.7e-7 * M * m * x_3**4 + + class KalbachMann(AngleEnergy): """Kalbach-Mann distribution @@ -319,7 +554,7 @@ class KalbachMann(AngleEnergy): n_energy_out = int(ace.xss[idx + 1]) data = ace.xss[idx + 2:idx + 2 + 5*n_energy_out].copy() data.shape = (5, n_energy_out) - data[0,:] *= EV_PER_MEV + data[0, :] *= EV_PER_MEV # Create continuous distribution eout_continuous = Tabular(data[0][n_discrete_lines:], @@ -352,13 +587,28 @@ class KalbachMann(AngleEnergy): return cls(breakpoints, interpolation, energy, energy_out, km_r, km_a) @classmethod - def from_endf(cls, file_obj): - """Generate Kalbach-Mann distribution from an ENDF evaluation + def from_endf(cls, file_obj, za_emitted, za_target, projectile_mass): + """Generate Kalbach-Mann distribution from an ENDF evaluation. + + If the projectile is a neutron, the slope is calculated when it is + not given explicitly. Parameters ---------- file_obj : file-like object ENDF file positioned at the start of the Kalbach-Mann distribution + za_emitted : int + 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. Returns ------- @@ -374,6 +624,7 @@ class KalbachMann(AngleEnergy): energy_out = [] precompound = [] slope = [] + calculated_slope = [] for i in range(ne): items, values = get_list_record(file_obj) energy[i] = items[1] @@ -385,19 +636,46 @@ class KalbachMann(AngleEnergy): values.shape = (n_energy_out, n_angle + 2) # Outgoing energy distribution at the i-th incoming energy - eout_i = values[:,0] - eout_p_i = values[:,1] + eout_i = values[:, 0] + eout_p_i = values[:, 1] energy_out_i = Tabular(eout_i, eout_p_i, INTERPOLATION_SCHEME[lep]) energy_out.append(energy_out_i) - # Precompound and slope factors for Kalbach-Mann - r_i = values[:,2] + # Precompound factors for Kalbach-Mann + r_i = values[:, 2] + + # Slope factors for Kalbach-Mann if n_angle == 2: - a_i = values[:,3] + a_i = values[:, 3] + calculated_slope.append(False) else: - a_i = np.zeros_like(r_i) + # 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) + precompound.append(Tabulated1D(eout_i, r_i)) slope.append(Tabulated1D(eout_i, a_i)) - return cls(tab2.breakpoints, tab2.interpolation, energy, - energy_out, precompound, slope) + km_distribution = cls(tab2.breakpoints, tab2.interpolation, energy, + energy_out, precompound, slope) + + # List of bool to indicate slope calculation by OpenMC + km_distribution._calculated_slope = calculated_slope + + return km_distribution diff --git a/openmc/data/njoy.py b/openmc/data/njoy.py index 305edaf4f..d8ecaae18 100644 --- a/openmc/data/njoy.py +++ b/openmc/data/njoy.py @@ -127,7 +127,7 @@ acer / %%%%%%%%%%%%%%%%%%%%%%%% Write out in ACE format %%%%%%%%%%%%%%%%%%%%%%%% 1 0 1 .{ext} / '{library}: {zsymam} at {temperature}'/ {mat} {temperature} -1 1/ +1 1 {ismooth}/ / """ @@ -248,7 +248,8 @@ def make_pendf(filename, pendf='pendf', error=0.001, stdout=False): def make_ace(filename, temperatures=None, acer=True, xsdir=None, output_dir=None, pendf=False, error=0.001, broadr=True, - heatr=True, gaspr=True, purr=True, evaluation=None, **kwargs): + heatr=True, gaspr=True, purr=True, evaluation=None, + smoothing=True, **kwargs): """Generate incident neutron ACE file from an ENDF file File names can be passed to @@ -298,6 +299,8 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None, evaluation : openmc.data.endf.Evaluation, optional If the ENDF file contains multiple material evaluations, this argument indicates which evaluation should be used. + smoothing : bool, optional + If the smoothing option (ACER card 6) is on (True) or off (False). **kwargs Keyword arguments passed to :func:`openmc.data.njoy.run` @@ -380,6 +383,7 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None, # acer if acer: + ismooth = int(smoothing) nacer_in = nlast for i, temperature in enumerate(temperatures): # Extend input with an ACER run for each temperature diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index 5e4287f16..a2431cf1c 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -80,6 +80,14 @@ def _get_products(ev, mt): mt : int The MT value of the reaction to get products for + Raises + ------ + IOError + When the Kalbach-Mann systematics is used, but the product + 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. + Returns ------- products : list of openmc.data.Product @@ -141,7 +149,26 @@ def _get_products(ev, mt): if lang == 1: p.distribution = [CorrelatedAngleEnergy.from_endf(file_obj)] elif lang == 2: - p.distribution = [KalbachMann.from_endf(file_obj)] + # Products need to be described in the center-of-mass system + product_center_of_mass = False + if reference_frame == 'center-of-mass': + product_center_of_mass = True + elif reference_frame == 'light-heavy': + product_center_of_mass = (awr <= 4.0) + # TODO: 'breakup' logic not implemented + + if product_center_of_mass is False: + raise IOError( + "Kalbach-Mann representation must be defined in the " + "'center-of-mass' system" + ) + + zat = ev.target["atomic_number"] * 1000 + ev.target["mass_number"] + projectile_mass = ev.projectile["mass"] + p.distribution = [KalbachMann.from_endf(file_obj, + za, + zat, + projectile_mass)] elif law == 2: # Discrete two-body scattering diff --git a/tests/unit_tests/test_data_kalbach_mann.py b/tests/unit_tests/test_data_kalbach_mann.py new file mode 100644 index 000000000..3b837af7d --- /dev/null +++ b/tests/unit_tests/test_data_kalbach_mann.py @@ -0,0 +1,176 @@ +"""Test of the Kalbach-Mann slope calculation when data are +retrieved from ENDF files.""" + +import os +from pathlib import Path +import pytest + +import numpy as np + +from openmc.data import IncidentNeutron +from openmc.data.kalbach_mann import _separation_energy, _AtomicRepresentation +from openmc.data import kalbach_slope +from openmc.data import KalbachMann + +from . import needs_njoy + + +@pytest.fixture(scope='module') +def neutron(): + """Neutron AtomicRepresentation.""" + return _AtomicRepresentation(z=0, a=1) + + +@pytest.fixture(scope='module') +def triton(): + """Triton AtomicRepresentation.""" + return _AtomicRepresentation(z=1, a=3) + + +@pytest.fixture(scope='module') +def b10(): + """B10 AtomicRepresentation.""" + return _AtomicRepresentation(z=5, a=10) + + +@pytest.fixture(scope='module') +def c12(): + """C12 AtomicRepresentation.""" + return _AtomicRepresentation(z=6, a=12) + + +@pytest.fixture(scope='module') +def c13(): + """C13 AtomicRepresentation.""" + return _AtomicRepresentation(z=6, a=13) + + +@pytest.fixture(scope='module') +def na23(): + """Na23 AtomicRepresentation.""" + return _AtomicRepresentation(z=11, a=23) + + +def test_atomic_representation(neutron, triton, b10, c12, c13, na23): + """Test the _AtomicRepresentation class.""" + # Test instantiation from_za + assert b10 == _AtomicRepresentation.from_za(5010) + + # Test addition + assert c13 + b10 == na23 + + # Test substraction + assert c13 - c12 == neutron + assert c13 - b10 == triton + + # Test properties when no information for Kalbach-Mann are given + assert c13.a == 13 + assert c13.z == 6 + assert c13.n == 7 + assert c13.za == 6013 + + # Test properties when information for Kalbach-Mann are given + assert triton.a == 3 + assert triton.z == 1 + assert triton.n == 2 + assert triton.za == 1003 + + # Test instantiation errors + with pytest.raises(ValueError): + _AtomicRepresentation(z=5, a=1) + with pytest.raises(ValueError): + _AtomicRepresentation(z=-1, a=1) + with pytest.raises(ValueError): + _AtomicRepresentation(z=5, a=0) + with pytest.raises(ValueError): + _AtomicRepresentation(z=5, a=-2) + with pytest.raises(ValueError): + neutron - triton + + +def test_separation_energy(triton, b10, c13): + """Comparison to hand-calculations on a simple example.""" + assert _separation_energy( + compound=c13, + nucleus=b10, + particle=triton + ) == pytest.approx(18.6880713) + + +def test_kalbach_slope(): + """Comparison to hand-calculations for n + c12 -> c13 -> triton + b10.""" + 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 + ) + + assert kalbach_slope( + energy_projectile=energy_projectile, + energy_emitted=energy_emitted, + za_projectile=1, + za_emitted=1003, + za_target=6012 + ) == pytest.approx(0.8409921475) + + +@pytest.mark.parametrize( + "hdf5_filename, endf_filename", [ + ('O16.h5', 'n-008_O_016.endf'), + ('Ca46.h5', 'n-020_Ca_046.endf'), + ('Hg204.h5', 'n-080_Hg_204.endf') + ] +) +def test_comparison_slope_hdf5(hdf5_filename, endf_filename): + """Test the calculation of the Kalbach-Mann slope done by OpenMC + by comparing it to HDF5 data. The test is based on the first product + of MT=5 (neutron). The isotopes tested have been selected because the + corresponding products in ENDF/B-VII.1 are described using MF=6, LAW=1, + LANG=2 (i.e., Kalbach-Mann systematics) and the slope is not given + explicitly. + + If an error occurs during the "validity check", this means that + the nuclear data evaluation has evolved and the distribution might + no longer be described using Kalbach-Mann systematics. Another + isotope needs to be identified and tested. + + Warning: This test is valid as long as ENDF files are not directly + used to generate the HDF5 files used in the tests. + + """ + # HDF5 data + hdf5_directory = Path(os.environ['OPENMC_CROSS_SECTIONS']).parent + hdf5_data = IncidentNeutron.from_hdf5(hdf5_directory / hdf5_filename) + hdf5_product = hdf5_data[5].products[0] + hdf5_distribution = hdf5_product.distribution[0] + + # ENDF data + endf_directory = Path(os.environ['OPENMC_ENDF_DATA']) + endf_path = endf_directory / 'neutrons' / endf_filename + endf_data = IncidentNeutron.from_endf(endf_path) + endf_product = endf_data[5].products[0] + endf_distribution = endf_product.distribution[0] + + # Validity check + assert isinstance(endf_distribution, KalbachMann) + assert isinstance(hdf5_distribution, KalbachMann) + assert endf_product.particle == hdf5_product.particle + assert len(endf_distribution.slope) == len(hdf5_distribution.slope) + + # Results check + for i, hdf5_slope in enumerate(hdf5_distribution.slope): + assert endf_distribution._calculated_slope[i] + + np.testing.assert_array_almost_equal( + endf_distribution.slope[i].y, + hdf5_slope.y, + decimal=5 + )