From cbe9cc8c290d675b515fa40507fd4179c74a9b5b Mon Sep 17 00:00:00 2001 From: Joffrey Dorville Date: Fri, 1 Apr 2022 17:31:16 -0600 Subject: [PATCH] Kalbach-Mann slope calculation for ENDF files when the projectile is a neutron --- openmc/data/kalbach_mann.py | 528 ++++++++++++++++++++- openmc/data/njoy.py | 11 +- openmc/data/reaction.py | 29 +- tests/unit_tests/test_data_kalbach_mann.py | 280 +++++++++++ 4 files changed, 834 insertions(+), 14 deletions(-) create mode 100644 tests/unit_tests/test_data_kalbach_mann.py diff --git a/openmc/data/kalbach_mann.py b/openmc/data/kalbach_mann.py index f4c914d909..11e3acf40c 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,468 @@ from .data import EV_PER_MEV from .endf import get_list_record, get_tab2_record +# Kalbach-Mann constants as defined in ENDF-6 manual BNL-203218-2018-INRE, +# Revision 215, File 6 description for LAW=1 and LANG=2. +_C1 = 0.04 # [1/MeV] +_C2 = 1.8E-6 # [1/MeV^3] +_C3 = 6.7E-7 # [1/MeV^4] +_ET1 = 130. # [MeV] +_ET3 = 41. # [MeV] +_M_NEUTRON = 1. +_M_PROTON = 1. +_M_DEUTERON = 1. +_M_TRITON = None +_M_3HE = None +_M_ALPHA = 0. +_SM_NEUTRON = 1/2. +_SM_PROTON = 1. +_SM_DEUTERON = 1. +_SM_TRITON = 1. +_SM_3HE = 1. +_SM_ALPHA = 2. + +# Kalbach-Mann M coefficients +_TABULATED_PARTICLE_M = { + 1: _M_NEUTRON, + 1001: _M_PROTON, + 1002: _M_DEUTERON, + 1003: _M_TRITON, + 2003: _M_3HE, + 2004: _M_ALPHA +} + +# Kalbach-Mann m coefficients +_TABULATED_PARTICLE_SM = { + 1: _SM_NEUTRON, + 1001: _SM_PROTON, + 1002: _SM_DEUTERON, + 1003: _SM_TRITON, + 2003: _SM_3HE, + 2004: _SM_ALPHA +} + +# Breaking energy as defined in ENDF-6 manual BNL-203218-2018-INRE, +# Revision 215, Appendix H, Table 3. +_BREAKING_ENERGY_NEUTRON = 0. +_BREAKING_ENERGY_PROTON = 0. +_BREAKING_ENERGY_DEUTERON = 2.224566 # [MeV] +_BREAKING_ENERGY_TRITON = 8.481798 # [MeV] +_BREAKING_ENERGY_3HE = 7.718043 # [MeV] +_BREAKING_ENERGY_ALPHA = 28.29566 # [MeV] + +_TABULATED_BREAKING_ENERGY = { + 1: _BREAKING_ENERGY_NEUTRON, + 1001: _BREAKING_ENERGY_PROTON, + 1002: _BREAKING_ENERGY_DEUTERON, + 1003: _BREAKING_ENERGY_TRITON, + 2003: _BREAKING_ENERGY_3HE, + 2004: _BREAKING_ENERGY_ALPHA +} + +# Abundant IZA translation in merged library +_IZA_TRANSLATION = { + 6000: 6012, +} + + +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 + ------ + IOError: + 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 + breaking_energy: float + Energy required to break the isotope or particle into their + constituent nucleons from tabulated values + M: float + Kalbach-Mann M coefficient + m: float + Kalbach-Mann m coefficient + iza: int + ZA identifier defined as: + iza = Z x 1000 + A, + where Z is the number of protons and A the number of nucleons + + """ + def __init__(self, z, a): + self._consistency_check(z, a) + self._z = z + self._a = a + self.z = self._z + self.a = self._a + + def __add__(self, other): + """Adds two AtomicRepresentations. + + """ + z = self.z + other.z + a = self.a + other.a + return AtomicRepresentation(z=z, a=a) + + def __sub__(self, other): + """Substracts two AtomicRepresentations. + + """ + z = self.z - other.z + a = self.a - other.a + return AtomicRepresentation(z=z, a=a) + + @property + def a(self): + self._consistency_check(self._z, self._a) + return self._a + + @property + def z(self): + self._consistency_check(self._z, self._a) + return self._z + + @property + def n(self): + return self.a - self.z + + @property + def breaking_energy(self): + breaking_energy = None + if self.iza in _TABULATED_BREAKING_ENERGY: + breaking_energy = _TABULATED_BREAKING_ENERGY[self.iza] + return breaking_energy + + @property + def M(self): + M = None + if self.iza in _TABULATED_PARTICLE_M: + M = _TABULATED_PARTICLE_M[self.iza] + return M + + @property + def m(self): + m = None + if self.iza in _TABULATED_PARTICLE_SM: + m = _TABULATED_PARTICLE_SM[self.iza] + return m + + @property + def iza(self): + iza = self.z * 1000 + self.a + return iza + + @a.setter + def a(self, an): + cv.check_type('a', an, Integral) + cv.check_greater_than('a', an, 0, equality=True) + self._consistency_check(self._z, an) + self._a = an + + @z.setter + def z(self, zn): + cv.check_type('z', zn, Integral) + cv.check_greater_than('z', zn, 0, equality=True) + self._consistency_check(zn, self._a) + self._z = zn + + @staticmethod + def _consistency_check(z, a): + """Simple consistency check. + + Parameters + ---------- + z: int + Number of protons (atomic number) + a: int + Number of nucleons (mass number) + + Raises + ------ + IOError: + When the number of protons (z) declared is higher than the number + of nucleons (a) + + """ + if z > a: + raise IOError( + "Number of protons (%i) incompatible with number of " + "nucleons (%i)" % (z, a) + ) + + @classmethod + def from_iza(cls, iza): + """Instantiates an AtomicRepresentation from a ZA identifier. + + The ZA identifier is defined as: + iza = Z x 1000 + A, + where Z is the number of protons and A the number of nucleons. + + Parameters + ---------- + iza: int + ZA identifier + + Returns + ------- + AtomicRepresentation + Atomic representation of the isotope/particle + + """ + if iza in _IZA_TRANSLATION.keys(): + iza = _IZA_TRANSLATION[iza] + + z = int(iza/1000) + a = np.mod(iza, 1000) + + return cls(z, a) + + +def _calculate_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 + + """ + coef_1 = 15.68 * (compound.a - nucleus.a) + coef_2 = 28.07 * ((compound.n - compound.z)**2 / float(compound.a) - \ + (nucleus.n - nucleus.z)**2 / float(nucleus.a)) + coef_3 = 18.56 * (compound.a**(2./3.) - nucleus.a**(2./3.)) + coef_4 = 33.22 * ((compound.n - compound.z)**2 / float(compound.a)**(4./3.) - \ + (nucleus.n - nucleus.z)**2 / float(nucleus.a)**(4./3.)) + coef_5 = 0.717 * (compound.z**2 / float(compound.a)**(1./3.) - \ + nucleus.z**2 / float(nucleus.a)**(1./3.)) + coef_6 = 1.211 * (compound.z**2 / float(compound.a) - \ + nucleus.z**2 / float(nucleus.a)) + + separation_energy = coef_1 - coef_2 - coef_3 + coef_4 \ + - coef_5 + coef_6 - particle.breaking_energy + + return separation_energy + + +def _return_entrance_channel_energy(e_p, awr_t, awr_p): + """Returns the entrance channel energy as defined in ENDF-6 manual + BNL-203218-2018-INRE, Revision 215, File 6 description for LAW=1 + and LANG=2. + + Parameters + ---------- + e_p: float + Energy of the incident projectile in the laboratory system in eV + awr_t: float + Atomic weight ratio of the target + awr_p: float + Atomic weight ratio of the projectile + + Returns + ------- + epsilon_p: float + Entrance channel energy in eV + + """ + epsilon_p = e_p * awr_t / (awr_t + awr_p) + return epsilon_p + + +def _return_emission_channel_energy(e_e, awr_r, awr_e): + """Returns the emission channel energy as defined in ENDF-6 manual + BNL-203218-2018-INRE, Revision 215, File 6 description for LAW=1 + and LANG=2. + + Parameters + ---------- + e_e: float + Energy of the emitted particle in the center of mass system in eV + awr_r: float + Atomic weight ratio of the residual nucleus + awr_e: float + Atomic weight ratio of the emitted particle + + Returns + ------- + epsilon_e: float + Emission channel energy in eV + + """ + epsilon_e = e_e * (awr_r + awr_e) / awr_r + return epsilon_e + + +def _calculate_kalbach_slope(energy_projectile, + energy_emitted, + projectile, + target, + compound, + emitted, + residual): + """Calculate the Kalbach slope for projectiles other than photons + as defined in ENDF-6 manual BNL-203218-2018-INRE, Revision 215, + File 6 description for LAW=1 and LANG=2. + + 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 + projectile: AtomicRepresentation + Atomic representation of the projectile + target: AtomicRepresentation + Atomic representation of the target + compound: AtomicRepresentation + Atomic representation of the compound + emitted: AtomicRepresentation + Atomic representation of the emitted particle + residual: AtomicRepresentation + Atomic representation of the residual nucleus + + Returns + ------- + slope: float + Kalbach-Mann slope + + """ + epsilon_a = _return_entrance_channel_energy( + energy_projectile, + target.a, + projectile.a + ) / EV_PER_MEV + epsilon_b = _return_emission_channel_energy( + energy_emitted, + residual.a, + emitted.a + ) / EV_PER_MEV + + s_a = _calculate_separation_energy(compound, target, projectile) + s_b = _calculate_separation_energy(compound, residual, emitted) + + e_a = epsilon_a + s_a + e_b = epsilon_b + s_b + + r_1 = min(e_a, _ET1) + r_3 = min(e_a, _ET3) + + x_1 = r_1 * e_b / e_a + x_3 = r_3 * e_b / e_a + + slope = _C1 * x_1 \ + + _C2 * x_1**3 \ + + _C3 * projectile.M * emitted.m * x_3**4 + + return slope + + +def return_kalbach_slope(energy_projectile, + energy_emitted, + iza_projectile, + iza_emitted, + iza_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. + + This function uses the concept of ZA identifier defined as: + iza = Z x 1000 + A, + where Z is the number of protons and A the number of nucleons. + + 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 + iza_projectile: int + ZA identifier of the projectile + iza_emitted: int + ZA identifier of the emitted particle + iza_target: int + ZA identifier of the targeted nucleus + + Raises + ------ + NotImplementedError: + When the ZA identifier of the projectile is not equal to 1 + (ie. other than 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 iza_projectile != 1: + raise NotImplementedError( + "Developed and tested for neutron projectile only." + ) + + projectile = AtomicRepresentation.from_iza(iza_projectile) + emitted = AtomicRepresentation.from_iza(iza_emitted) + target = AtomicRepresentation.from_iza(iza_target) + compound = projectile + target + residual = compound - emitted + + slope = _calculate_kalbach_slope( + energy_projectile, + energy_emitted, + projectile, + target, + compound, + emitted, + residual + ) + + return float("%7e" % slope) + + class KalbachMann(AngleEnergy): """Kalbach-Mann distribution @@ -319,7 +782,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 +815,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, iza_emitted, iza_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 + iza_emitted : int + ZA identifier of the emitted particle + iza_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 +852,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 +864,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 IZA of the projectile + iza_projectile = 1 + a_i = [return_kalbach_slope(energy_projectile=energy[i], + energy_emitted=e, + iza_projectile=iza_projectile, + iza_emitted=iza_emitted, + iza_target=iza_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 305edaf4f2..fbb26be05e 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 {ismoothing}/ / """ @@ -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 in ACER is on (1) or off (0) in the card 6. **kwargs Keyword arguments passed to :func:`openmc.data.njoy.run` @@ -380,6 +383,10 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None, # acer if acer: + if smoothing is True: + ismoothing = 1 + else: + ismoothing = 0 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 5e4287f16e..42fc576da2 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 0000000000..99869808de --- /dev/null +++ b/tests/unit_tests/test_data_kalbach_mann.py @@ -0,0 +1,280 @@ +"""Test of the Kalbach-Mann slope calculation when data are +retrieved from ENDF files.""" + +import os +import pytest +import numpy as np + +from openmc.data import IncidentNeutron +from openmc.data import AtomicRepresentation +from openmc.data.kalbach_mann import _BREAKING_ENERGY_TRITON +from openmc.data.kalbach_mann import _M_TRITON +from openmc.data.kalbach_mann import _SM_TRITON +from openmc.data.kalbach_mann import _calculate_separation_energy +from openmc.data.kalbach_mann import _return_emission_channel_energy +from openmc.data.kalbach_mann import _return_entrance_channel_energy +from openmc.data.kalbach_mann import _calculate_kalbach_slope +from openmc.data import return_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 instanciation from_iza + assert b10 == AtomicRepresentation.from_iza(5010) + + # Test instanciation from_iza using IZA translation + assert c12 == AtomicRepresentation.from_iza(6000) + + # 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.breaking_energy is None + assert c13.M is None + assert c13.m is None + assert c13.iza == 6013 + + # Test properties when information for Kalbach-Mann are given + assert triton.a == 3 + assert triton.z == 1 + assert triton.n == 2 + assert triton.breaking_energy == _BREAKING_ENERGY_TRITON + assert triton.M == _M_TRITON + assert triton.m == _SM_TRITON + assert triton.iza == 1003 + + # Test instanciation errors + with pytest.raises(IOError): + AtomicRepresentation(z=5, a=1) + with pytest.raises(ValueError): + AtomicRepresentation(z=-1, a=1) + with pytest.raises(IOError): + AtomicRepresentation(z=5, a=0) + with pytest.raises(IOError): + AtomicRepresentation(z=5, a=-2) + with pytest.raises(OSError): + neutron - triton + + +def test__calculate_separation_energy(triton, b10, c13): + """Comparison to hand-calculations on a simple example.""" + assert _calculate_separation_energy( + compound=c13, + nucleus=b10, + particle=triton + ) == pytest.approx(18.6880713) + + +def test__return_entrance_channel_energy(): + """Comparison to hand-calculations on a simple example.""" + assert _return_entrance_channel_energy( + e_p=5.2, + awr_t=13.7, + awr_p=7.2 + ) == pytest.approx(3.4086124) + + +def test__return_emission_channel_energy(): + """Comparison to hand-calculations on a simple example.""" + assert _return_emission_channel_energy( + e_e=6.8, + awr_r=49.2, + awr_e=5.4 + ) == pytest.approx(7.5463415) + + +def test__calculate_kalbach_slope(neutron, triton, b10, c12, c13): + """Comparison to hand-calculations for n + c12 -> c13 -> triton + b10.""" + energy_projectile = 10.2 # [eV] + energy_emitted = 5.4 # [eV] + + assert _calculate_kalbach_slope( + energy_projectile=energy_projectile, + energy_emitted=energy_emitted, + projectile=neutron, + target=c12, + compound=c13, + emitted=triton, + residual=b10 + ) == pytest.approx(0.8409921475) + + +def test_return_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): + return_kalbach_slope( + energy_projectile=energy_projectile, + energy_emitted=energy_emitted, + iza_projectile=1000, + iza_emitted=1, + iza_target=6012 + ) + + assert return_kalbach_slope( + energy_projectile=energy_projectile, + energy_emitted=energy_emitted, + iza_projectile=1, + iza_emitted=1003, + iza_target=6012 + ) == pytest.approx(0.8409921475) + + +@pytest.mark.parametrize( + "hdf5_filename, endf_type, endf_filename", [ + ('O16.h5', 'neutrons', 'n-008_O_016.endf'), + ('Ca46.h5', 'neutrons', 'n-020_Ca_046.endf'), + ('Hg204.h5', 'neutrons', 'n-080_Hg_204.endf') + ] +) +def test_comparison_slope_hdf5(hdf5_filename, endf_type, 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 (ie. 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 = os.path.dirname(os.environ['OPENMC_CROSS_SECTIONS']) + hdf5_path = os.path.join(hdf5_directory, hdf5_filename) + hdf5_data = IncidentNeutron.from_hdf5(hdf5_path) + hdf5_product = hdf5_data[5].products[0] + hdf5_distribution = hdf5_product.distribution[0] + + # ENDF data + endf_directory = os.environ['OPENMC_ENDF_DATA'] + endf_path = os.path.join(endf_directory, endf_type, 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] is True + + np.testing.assert_array_almost_equal( + endf_distribution.slope[i].y, + hdf5_slope.y, + decimal=6 + ) + + +@needs_njoy +@pytest.mark.parametrize( + "endf_type, endf_filename", [ + ('neutrons', 'n-008_O_016.endf'), + ('neutrons', 'n-020_Ca_046.endf'), + ('neutrons', 'n-080_Hg_204.endf') + ] +) +def test_comparison_slope_njoy(endf_type, endf_filename): + """Test the calculation of the Kalbach-Mann slope done by OpenMC + by comparing it to an NJOY calculation. 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 (ie. 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. + + """ + endf_directory = os.environ['OPENMC_ENDF_DATA'] + endf_path = os.path.join(endf_directory, endf_type, endf_filename) + + # ENDF data + endf_data = IncidentNeutron.from_endf(endf_path) + endf_product = endf_data[5].products[0] + endf_distribution = endf_product.distribution[0] + + # NJOY data + njoy_data = IncidentNeutron.from_njoy(endf_path, heatr=False, gaspr=False, + purr=False, smoothing=False) + njoy_product = njoy_data[5].products[0] + njoy_distribution = njoy_product.distribution[0] + + # Validity check + assert isinstance(endf_distribution, KalbachMann) + assert isinstance(njoy_distribution, KalbachMann) + assert endf_product.particle == njoy_product.particle + assert len(endf_distribution.slope) == len(njoy_distribution.slope) + + # Results check + for i, njoy_slope in enumerate(njoy_distribution.slope): + + assert endf_distribution._calculated_slope[i] is True + + np.testing.assert_array_almost_equal( + endf_distribution.slope[i].y, + njoy_slope.y, + decimal=6 + )