diff --git a/openmc/data/kalbach_mann.py b/openmc/data/kalbach_mann.py index 11e3acf40..5c0423f30 100644 --- a/openmc/data/kalbach_mann.py +++ b/openmc/data/kalbach_mann.py @@ -13,70 +13,6 @@ 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. @@ -151,27 +87,6 @@ class AtomicRepresentation(EqualityMixin): 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 @@ -234,16 +149,11 @@ class AtomicRepresentation(EqualityMixin): 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) - + z, a = divmod(iza, 1000) return cls(z, a) -def _calculate_separation_energy(compound, nucleus, particle): +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 @@ -251,158 +161,56 @@ def _calculate_separation_energy(compound, nucleus, particle): Parameters ---------- - compound: AtomicRepresentation + compound : AtomicRepresentation Atomic representation of the compound (C) - nucleus: AtomicRepresentation + nucleus : AtomicRepresentation Atomic representation of the nucleus (A or B) - particle: AtomicRepresentation + particle : AtomicRepresentation Atomic representation of the particle (a or b) Returns ------- - separation_energy: float + 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)) + # 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 - separation_energy = coef_1 - coef_2 - coef_3 + coef_4 \ - - coef_5 + coef_6 - particle.breaking_energy + # Determine breakup energy of incident particle (ENDF-6 Formats Manual, + # Appendix H, Table 3) + iza_to_breaking_energy = { + 1: 0.0, + 1001: 0.0, + 1002: 2.224566, + 1003: 8.481798, + 2003: 7.718043, + 2004: 28.29566 + } + I_b = iza_to_breaking_energy[particle.iza] - return separation_energy + # 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_b + ) -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): +def 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 @@ -456,21 +264,42 @@ def return_kalbach_slope(energy_projectile, "Developed and tested for neutron projectile only." ) + # Special handling of elemental carbon + if iza_emitted == 6000: + iza_emitted = 6012 + if iza_target == 6000: + iza_target = 6012 + 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 - ) + # 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 + iza_to_M = {1: 1.0, 1001: 1.0, 1002: 1.0, 2004: 0.0} + iza_to_m = {1: 0.5, 1001: 1.0, 1002: 1.0, 1003: 1.0, 2003: 1.0, 2004: 2.0} + M = iza_to_M[projectile.iza] + m = iza_to_m[emitted.iza] + 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 + slope = 0.04 * x_1 + 1.8e-6 * x_1**3 + 6.7e-7 * M * m * x_3**4 return float("%7e" % slope) @@ -889,11 +718,11 @@ class KalbachMann(AngleEnergy): 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) + a_i = [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) diff --git a/tests/unit_tests/test_data_kalbach_mann.py b/tests/unit_tests/test_data_kalbach_mann.py index 99869808d..643fffcfc 100644 --- a/tests/unit_tests/test_data_kalbach_mann.py +++ b/tests/unit_tests/test_data_kalbach_mann.py @@ -7,14 +7,8 @@ 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.kalbach_mann import _separation_energy +from openmc.data import kalbach_slope from openmc.data import KalbachMann from . import needs_njoy @@ -61,9 +55,6 @@ def test_atomic_representation(neutron, triton, b10, c12, c13, na23): # 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 @@ -75,18 +66,12 @@ def test_atomic_representation(neutron, triton, b10, c12, c13, na23): 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 @@ -102,50 +87,16 @@ def test_atomic_representation(neutron, triton, b10, c12, c13, na23): neutron - triton -def test__calculate_separation_energy(triton, b10, c13): +def test_separation_energy(triton, b10, c13): """Comparison to hand-calculations on a simple example.""" - assert _calculate_separation_energy( + assert _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(): +def test_kalbach_slope(): """Comparison to hand-calculations for n + c12 -> c13 -> triton + b10.""" energy_projectile = 10.2 # [eV] energy_emitted = 5.4 # [eV] @@ -153,7 +104,7 @@ def test_return_kalbach_slope(): # Check that NotImplementedError is raised if the projectile is not # a neutron with pytest.raises(NotImplementedError): - return_kalbach_slope( + kalbach_slope( energy_projectile=energy_projectile, energy_emitted=energy_emitted, iza_projectile=1000, @@ -161,7 +112,7 @@ def test_return_kalbach_slope(): iza_target=6012 ) - assert return_kalbach_slope( + assert kalbach_slope( energy_projectile=energy_projectile, energy_emitted=energy_emitted, iza_projectile=1,