crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 05:35:49 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge pull request #2022 from JoffreyDorville/kalbach_mann_slope

Browse source 
Kalbach-Mann slope calculation for ENDF files
This commit is contained in:
Paul Romano 2022-08-05 11:51:00 -05:00 • committed by GitHub
commit 2adf34b9a6
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
5 changed files with 500 additions and 14 deletions
Download patch file Download diff file Expand all files Collapse all files

1
docs/source/pythonapi/data.rst
View file

@ -67,6 +67,7 @@ Core Functions
gnd_name
half_life
isotopes
kalbach_slope
linearize
thin
water_density

300
openmc/data/kalbach_mann.py
View file

@ -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

8
openmc/data/njoy.py
View file

@ -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

29
openmc/data/reaction.py
View file

@ -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

176
tests/unit_tests/test_data_kalbach_mann.py Normal file
View file

@ -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
)