mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 22:26:08 -04:00
Merge KM refactoring
Refactor of KM slope branch
This commit is contained in:
Joffrey Dorville
GitHub
commit
64b7889cec
4 changed files with 171 additions and 507 deletions
|
|
@ -24,7 +24,6 @@ and product yields.
|
|||
FissionProductYields
|
||||
WindowedMultipole
|
||||
ProbabilityTables
|
||||
AtomicRepresentation
|
||||
|
||||
The following classes are used for storing atomic data (incident photon cross
|
||||
sections, atomic relaxation):
|
||||
|
|
@ -65,11 +64,11 @@ Core Functions
|
|||
dose_coefficients
|
||||
gnd_name
|
||||
isotopes
|
||||
kalbach_slope
|
||||
linearize
|
||||
thin
|
||||
water_density
|
||||
zam
|
||||
return_kalbach_slope
|
||||
|
||||
One-dimensional Functions
|
||||
-------------------------
|
||||
|
|
|
|||
|
|
@ -13,138 +13,66 @@ 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):
|
||||
class _AtomicRepresentation(EqualityMixin):
|
||||
"""Atomic representation of an isotope or a particle.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
z: int
|
||||
z : int
|
||||
Number of protons (atomic number)
|
||||
a: int
|
||||
a : int
|
||||
Number of nucleons (mass number)
|
||||
|
||||
Raises
|
||||
------
|
||||
IOError:
|
||||
ValueError
|
||||
When the number of protons (z) declared is higher than the number
|
||||
of nucleons (a)
|
||||
|
||||
Attributes
|
||||
----------
|
||||
z: int
|
||||
z : int
|
||||
Number of protons (atomic number)
|
||||
a: int
|
||||
a : int
|
||||
Number of nucleons (mass number)
|
||||
n: int
|
||||
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
|
||||
za : int
|
||||
ZA identifier, 1000*Z + A, where Z is the atomic number and A the mass
|
||||
number
|
||||
|
||||
"""
|
||||
def __init__(self, z, a):
|
||||
self._consistency_check(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
|
||||
self.z = self._z
|
||||
self.a = self._a
|
||||
|
||||
def __add__(self, other):
|
||||
"""Adds two AtomicRepresentations.
|
||||
|
||||
"""
|
||||
"""Add two _AtomicRepresentations"""
|
||||
z = self.z + other.z
|
||||
a = self.a + other.a
|
||||
return AtomicRepresentation(z=z, a=a)
|
||||
return _AtomicRepresentation(z=z, a=a)
|
||||
|
||||
def __sub__(self, other):
|
||||
"""Substracts two AtomicRepresentations.
|
||||
|
||||
"""
|
||||
"""Substract two _AtomicRepresentations"""
|
||||
z = self.z - other.z
|
||||
a = self.a - other.a
|
||||
return AtomicRepresentation(z=z, a=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
|
||||
|
|
@ -152,98 +80,30 @@ class AtomicRepresentation(EqualityMixin):
|
|||
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)
|
||||
)
|
||||
def za(self):
|
||||
return self.z * 1000 + self.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.
|
||||
def from_za(cls, za):
|
||||
"""Instantiate an _AtomicRepresentation from a ZA identifier.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
iza: int
|
||||
ZA identifier
|
||||
za : int
|
||||
ZA identifier, 1000*Z + A, where Z is the atomic number and A the
|
||||
mass number
|
||||
|
||||
Returns
|
||||
-------
|
||||
AtomicRepresentation
|
||||
_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)
|
||||
|
||||
z, a = divmod(za, 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 +111,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) 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]
|
||||
|
||||
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_a
|
||||
)
|
||||
|
||||
|
||||
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, za_projectile,
|
||||
za_emitted, za_target):
|
||||
"""Returns Kalbach-Mann slope from calculations.
|
||||
|
||||
|
||||
The associated reaction is defined as:
|
||||
A + a -> C -> B + b
|
||||
|
||||
|
|
@ -414,10 +172,6 @@ def return_kalbach_slope(energy_projectile,
|
|||
- 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
|
||||
|
|
@ -426,53 +180,71 @@ def return_kalbach_slope(energy_projectile,
|
|||
|
||||
Parameters
|
||||
----------
|
||||
energy_projectile: float
|
||||
energy_projectile : float
|
||||
Energy of the projectile in the laboratory system in eV
|
||||
energy_emitted: float
|
||||
energy_emitted : float
|
||||
Energy of the emitted particle in the center of mass system in eV
|
||||
iza_projectile: int
|
||||
za_projectile : int
|
||||
ZA identifier of the projectile
|
||||
iza_emitted: int
|
||||
za_emitted : int
|
||||
ZA identifier of the emitted particle
|
||||
iza_target: int
|
||||
za_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).
|
||||
NotImplementedError
|
||||
When the projectile is not a neutron
|
||||
|
||||
Returns
|
||||
-------
|
||||
slope: float
|
||||
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:
|
||||
if za_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)
|
||||
# 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
|
||||
|
||||
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)
|
||||
|
||||
return float("%7e" % slope)
|
||||
# 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):
|
||||
|
|
@ -815,7 +587,7 @@ class KalbachMann(AngleEnergy):
|
|||
return cls(breakpoints, interpolation, energy, energy_out, km_r, km_a)
|
||||
|
||||
@classmethod
|
||||
def from_endf(cls, file_obj, iza_emitted, iza_target, projectile_mass):
|
||||
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
|
||||
|
|
@ -825,11 +597,11 @@ class KalbachMann(AngleEnergy):
|
|||
----------
|
||||
file_obj : file-like object
|
||||
ENDF file positioned at the start of the Kalbach-Mann distribution
|
||||
iza_emitted : int
|
||||
za_emitted : int
|
||||
ZA identifier of the emitted particle
|
||||
iza_target : int
|
||||
za_target : int
|
||||
ZA identifier of the target
|
||||
projectile_mass: float
|
||||
projectile_mass : float
|
||||
Mass of the projectile
|
||||
|
||||
Warns
|
||||
|
|
@ -887,13 +659,13 @@ class KalbachMann(AngleEnergy):
|
|||
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)
|
||||
# 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)
|
||||
|
||||
|
|
|
|||
|
|
@ -300,7 +300,7 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None,
|
|||
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.
|
||||
If the smoothing option (ACER card 6) is on (True) or off (False).
|
||||
**kwargs
|
||||
Keyword arguments passed to :func:`openmc.data.njoy.run`
|
||||
|
||||
|
|
@ -383,10 +383,7 @@ def make_ace(filename, temperatures=None, acer=True, xsdir=None,
|
|||
|
||||
# acer
|
||||
if acer:
|
||||
if smoothing is True:
|
||||
ismoothing = 1
|
||||
else:
|
||||
ismoothing = 0
|
||||
ismoothing = int(smoothing)
|
||||
nacer_in = nlast
|
||||
for i, temperature in enumerate(temperatures):
|
||||
# Extend input with an ACER run for each temperature
|
||||
|
|
|
|||
|
|
@ -2,19 +2,14 @@
|
|||
retrieved from ENDF files."""
|
||||
|
||||
import os
|
||||
from pathlib import Path
|
||||
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.kalbach_mann import _separation_energy, _AtomicRepresentation
|
||||
from openmc.data import kalbach_slope
|
||||
from openmc.data import KalbachMann
|
||||
|
||||
from . import needs_njoy
|
||||
|
|
@ -23,46 +18,43 @@ from . import needs_njoy
|
|||
@pytest.fixture(scope='module')
|
||||
def neutron():
|
||||
"""Neutron AtomicRepresentation."""
|
||||
return AtomicRepresentation(z=0, a=1)
|
||||
return _AtomicRepresentation(z=0, a=1)
|
||||
|
||||
|
||||
@pytest.fixture(scope='module')
|
||||
def triton():
|
||||
"""Triton AtomicRepresentation."""
|
||||
return AtomicRepresentation(z=1, a=3)
|
||||
return _AtomicRepresentation(z=1, a=3)
|
||||
|
||||
|
||||
@pytest.fixture(scope='module')
|
||||
def b10():
|
||||
"""B10 AtomicRepresentation."""
|
||||
return AtomicRepresentation(z=5, a=10)
|
||||
return _AtomicRepresentation(z=5, a=10)
|
||||
|
||||
|
||||
@pytest.fixture(scope='module')
|
||||
def c12():
|
||||
"""C12 AtomicRepresentation."""
|
||||
return AtomicRepresentation(z=6, a=12)
|
||||
return _AtomicRepresentation(z=6, a=12)
|
||||
|
||||
|
||||
@pytest.fixture(scope='module')
|
||||
def c13():
|
||||
"""C13 AtomicRepresentation."""
|
||||
return AtomicRepresentation(z=6, a=13)
|
||||
return _AtomicRepresentation(z=6, a=13)
|
||||
|
||||
|
||||
@pytest.fixture(scope='module')
|
||||
def na23():
|
||||
"""Na23 AtomicRepresentation."""
|
||||
return AtomicRepresentation(z=11, a=23)
|
||||
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 the _AtomicRepresentation class."""
|
||||
# Test instantiation from_za
|
||||
assert b10 == _AtomicRepresentation.from_za(5010)
|
||||
|
||||
# Test addition
|
||||
assert c13 + b10 == na23
|
||||
|
|
@ -75,77 +67,37 @@ 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
|
||||
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.breaking_energy == _BREAKING_ENERGY_TRITON
|
||||
assert triton.M == _M_TRITON
|
||||
assert triton.m == _SM_TRITON
|
||||
assert triton.iza == 1003
|
||||
assert triton.za == 1003
|
||||
|
||||
# Test instanciation errors
|
||||
with pytest.raises(IOError):
|
||||
AtomicRepresentation(z=5, a=1)
|
||||
# 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):
|
||||
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):
|
||||
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,36 +105,36 @@ 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,
|
||||
iza_emitted=1,
|
||||
iza_target=6012
|
||||
za_projectile=1000,
|
||||
za_emitted=1,
|
||||
za_target=6012
|
||||
)
|
||||
|
||||
assert return_kalbach_slope(
|
||||
assert kalbach_slope(
|
||||
energy_projectile=energy_projectile,
|
||||
energy_emitted=energy_emitted,
|
||||
iza_projectile=1,
|
||||
iza_emitted=1003,
|
||||
iza_target=6012
|
||||
za_projectile=1,
|
||||
za_emitted=1003,
|
||||
za_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')
|
||||
"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_type, endf_filename):
|
||||
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 (ie. Kalbach-Mann systematics) and the slope is not given
|
||||
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
|
||||
|
|
@ -195,15 +147,14 @@ def test_comparison_slope_hdf5(hdf5_filename, endf_type, endf_filename):
|
|||
|
||||
"""
|
||||
# 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_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 = os.environ['OPENMC_ENDF_DATA']
|
||||
endf_path = os.path.join(endf_directory, endf_type, endf_filename)
|
||||
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]
|
||||
|
|
@ -216,65 +167,10 @@ def test_comparison_slope_hdf5(hdf5_filename, endf_type, endf_filename):
|
|||
|
||||
# Results check
|
||||
for i, hdf5_slope in enumerate(hdf5_distribution.slope):
|
||||
|
||||
assert endf_distribution._calculated_slope[i] is True
|
||||
assert endf_distribution._calculated_slope[i]
|
||||
|
||||
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
|
||||
decimal=5
|
||||
)
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue