mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-26 21:25:36 -04:00
Refactor K-M slope calculation
This commit is contained in:
Paul Romano
parent
8f564052e5
commit
af579895dd
2 changed files with 80 additions and 300 deletions
|
|
@ -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)
|
||||
|
||||
|
|
|
|||
|
|
@ -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,
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue