This commit is contained in:
GuySten 2026-07-21 01:07:03 +06:00 committed by GitHub
commit f1104afd6d
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
18 changed files with 1332 additions and 58 deletions

View file

@ -50,7 +50,7 @@ attributes:
:type:
The type of data contained in the file. Accepted values are 'neutron',
'thermal', 'photon', and 'wmp'.
'thermal', 'photon', 'photonuclear' and 'wmp'.
.. _depletion_element:

View file

@ -256,6 +256,82 @@ temperature-dependent data set. For example, the data set corresponding to
:Groups:
- **distribution** -- Format for angle-energy distributions are
detailed in :ref:`angle_energy`.
--------------------------
Incident Photonuclear Data
--------------------------
**/**
:Attributes: - **filetype** (*char[]*) -- String indicating the type of file
- **version** (*int[2]*) -- Major and minor version of the data
**/<nuclide name>/**
:Attributes: - **Z** (*int*) -- Atomic number
- **A** (*int*) -- Mass number. For a natural element, A=0 is given.
- **metastable** (*int*) -- Metastable state (0=ground, 1=first
excited, etc.)
- **atomic_weight_ratio** (*double*) -- Mass in units of neutron masses
- **n_reaction** (*int*) -- Number of reactions
:Datasets:
- **energy** (*double[]*) -- Energies in [eV] at which cross sections
are tabulated
**/<nuclide name>/**
**/<nuclide name>/reactions/reaction_<mt>/**
:Attributes: - **mt** (*int*) -- ENDF MT reaction number
- **label** (*char[]*) -- Name of the reaction
- **Q_value** (*double*) -- Q value in eV
- **center_of_mass** (*int*) -- Whether the reference frame for
scattering is center-of-mass (1) or laboratory (0)
- **n_product** (*int*) -- Number of reaction products
- **redundant** (*int*) -- Whether reaction is redundant
**/<nuclide name>/reactions/reaction_<mt>/**
:Datasets:
- **xs** (*double[]*) -- Cross section values tabulated against the
nuclide energy grid
:Attributes:
- **threshold_idx** (*int*) -- Index on the energy
grid that the reaction threshold corresponds to
**/<nuclide name>/reactions/reaction_<mt>/product_<j>/**
Reaction product data is described in :ref:`product`.
**/<nuclide name>/fission_energy_release/**
:Datasets: - **fragments** (:ref:`function <1d_functions>`) -- Energy
released in the form of fragments as a function of incident
neutron energy.
- **prompt_neutrons** (:ref:`function <1d_functions>`) -- Energy
released in the form of prompt neutrons as a function of incident
neutron energy.
- **delayed_neutrons** (:ref:`function <1d_functions>`) -- Energy
released in the form of delayed neutrons as a function of incident
neutron energy.
- **prompt_photons** (:ref:`function <1d_functions>`) -- Energy
released in the form of prompt photons as a function of incident
neutron energy.
- **delayed_photons** (:ref:`function <1d_functions>`) -- Energy
released in the form of delayed photons as a function of incident
neutron energy.
- **betas** (:ref:`function <1d_functions>`) -- Energy released in
the form of betas as a function of incident neutron energy.
- **neutrinos** (:ref:`function <1d_functions>`) -- Energy released
in the form of neutrinos as a function of incident neutron energy.
- **q_prompt** (:ref:`function <1d_functions>`) -- The prompt fission
Q-value (fragments + prompt neutrons + prompt photons - incident
energy)
- **q_recoverable** (:ref:`function <1d_functions>`) -- The
recoverable fission Q-value (Q_prompt + delayed neutrons + delayed
photons + betas)
.. _product:
@ -415,6 +491,7 @@ Kalbach-Mann
:Object type: Group
:Attributes: - **type** (*char[]*) -- 'kalbach-mann'
- **is_photon** (*bool*) -- Whether the incident particle is a photon.
:Datasets: - **energy** (*double[]*) -- Incoming energies at which distributions exist
:Attributes:
@ -586,7 +663,8 @@ Level Inelastic
:Attributes: - **type** (*char[]*) -- 'level'
- **threshold** (*double*) -- Energy threshold in the laboratory
system in eV
- **mass_ratio** (*double*) -- :math:`(A/(A + 1))^2`
- **mass_ratio** (*double*) -- for incident neutrons: :math:`(A/(A + 1))^2`,
while for incident photons: :math:`(A-1)/A`
Continuous Tabular
------------------

View file

@ -50,6 +50,15 @@ The following classes are used for storing thermal neutron scattering data:
CoherentElastic
IncoherentElastic
The following classes are used for storing incident photonuclear interaction data:
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
IncidentPhotonuclear
PhotonuclearReaction
Core Functions
--------------
@ -215,3 +224,4 @@ NJOY Interface
njoy.make_pendf
njoy.make_ace
njoy.make_ace_thermal
njoy.make_ace_photonuclear

View file

@ -20,6 +20,7 @@ public:
enum class Type {
neutron = 1,
photon = 3,
photonuclear = 6,
thermal = 2,
multigroup = 4,
wmp = 5

View file

@ -62,6 +62,7 @@ private:
tensor::Tensor<double> a; //!< Parameterized function
};
bool is_photon_; //!< Whether the projectile is a photon
int n_region_; //!< Number of interpolation regions
vector<int> breakpoints_; //!< Breakpoints between regions
vector<Interpolation> interpolation_; //!< Interpolation laws

View file

@ -12,6 +12,7 @@ WMP_VERSION = (WMP_VERSION_MAJOR, WMP_VERSION_MINOR)
from .data import *
from .neutron import *
from .photon import *
from .photonuclear import *
from .decay import *
from .reaction import *
from . import ace

View file

@ -219,6 +219,7 @@ DADZ = {
# Unit conversions
JOULE_PER_EV = 1.602176634e-19
EV_PER_AMU = 931.49410372e6 # eV/c^2 per amu
# Avogadro's constant
AVOGADRO = 6.02214076e23
@ -226,6 +227,9 @@ AVOGADRO = 6.02214076e23
# Neutron mass in units of amu
NEUTRON_MASS = 1.00866491595
# Neutron mass in units of eV/c^2
NEUTRON_MASS_EV = NEUTRON_MASS*EV_PER_AMU
# Used in atomic_mass function as a cache
_ATOMIC_MASS: dict[str, float] = {}

View file

@ -1253,7 +1253,8 @@ class ContinuousTabular(EnergyDistribution):
# Create continuous distribution
eout_continuous = Tabular(data[0][n_discrete_lines:],
data[1][n_discrete_lines:]/EV_PER_MEV,
INTERPOLATION_SCHEME[intt])
INTERPOLATION_SCHEME[intt],
ignore_negative=True)
eout_continuous.c = data[2][n_discrete_lines:]
# If discrete lines are present, create a mixture distribution

View file

@ -9,7 +9,7 @@ from openmc.mixin import EqualityMixin
from openmc.stats import Tabular, Univariate, Discrete, Mixture
from .function import Tabulated1D, INTERPOLATION_SCHEME
from .angle_energy import AngleEnergy
from .data import EV_PER_MEV
from .data import EV_PER_MEV, NEUTRON_MASS_EV
from .endf import get_list_record, get_tab2_record
@ -204,13 +204,14 @@ def kalbach_slope(energy_projectile, energy_emitted, za_projectile,
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."
)
if za_projectile == 0:
# Calculate slope for photons using Eq. 3 in doi:10.1080/18811248.1995.9731830
# or ENDF-6 Formats Manual section 6.2.3.2
slope_n = kalbach_slope(energy_projectile, energy_emitted, 1,
za_emitted, za_target)
return slope_n * np.sqrt(0.5*energy_projectile/NEUTRON_MASS_EV)*np.minimum(4,np.maximum(1,9.3/np.sqrt(energy_emitted/EV_PER_MEV)))
# Special handling of elemental carbon
if za_emitted == 6000:
za_emitted = 6012
@ -268,6 +269,8 @@ class KalbachMann(AngleEnergy):
slope : Iterable of openmc.data.Tabulated1D
Kalbach-Chadwick angular distribution slope value 'a' as a function of
outgoing energy for each incoming energy
is_photon : bool
Whether projectile is a photon, defaults to False
Attributes
----------
@ -285,14 +288,17 @@ class KalbachMann(AngleEnergy):
slope : Iterable of openmc.data.Tabulated1D
Kalbach-Chadwick angular distribution slope value 'a' as a function of
outgoing energy for each incoming energy
is_photon : bool
Whether projectile particle is a photon
"""
def __init__(self, breakpoints, interpolation, energy, energy_out,
precompound, slope):
precompound, slope, is_photon=False):
super().__init__()
self.breakpoints = breakpoints
self.interpolation = interpolation
self.is_photon = is_photon
self.energy = energy
self.energy_out = energy_out
self.precompound = precompound
@ -317,6 +323,15 @@ class KalbachMann(AngleEnergy):
cv.check_type('Kalbach-Mann interpolation', interpolation,
Iterable, Integral)
self._interpolation = interpolation
@property
def is_photon(self):
return self._is_photon
@is_photon.setter
def is_photon(self, is_photon):
cv.check_type('Kalbach-Mann is_photon', is_photon, bool)
self._is_photon = is_photon
@property
def energy(self):
@ -367,6 +382,7 @@ class KalbachMann(AngleEnergy):
"""
group.attrs['type'] = np.bytes_('kalbach-mann')
group.attrs['is_photon'] = self.is_photon
dset = group.create_dataset('energy', data=self.energy)
dset.attrs['interpolation'] = np.vstack((self.breakpoints,
@ -436,6 +452,7 @@ class KalbachMann(AngleEnergy):
Kalbach-Mann energy distribution
"""
is_photon = bool(group.attrs.get("is_photon", False))
interp_data = group['energy'].attrs['interpolation']
energy_breakpoints = interp_data[0, :]
energy_interpolation = interp_data[1, :]
@ -491,7 +508,7 @@ class KalbachMann(AngleEnergy):
slope.append(km_a)
return cls(energy_breakpoints, energy_interpolation,
energy, energy_out, precompound, slope)
energy, energy_out, precompound, slope, is_photon=is_photon)
@classmethod
def from_ace(cls, ace, idx, ldis):
@ -507,6 +524,8 @@ class KalbachMann(AngleEnergy):
ldis : int
Index in XSS array of the start of the energy distribution block
(e.g. JXS[11])
is_photon : bool
Whether projectile is photon
Returns
-------
@ -514,6 +533,7 @@ class KalbachMann(AngleEnergy):
Kalbach-Mann energy-angle distribution
"""
is_photon = bool(ace.data_type.value == 'u')
# Read number of interpolation regions and incoming energies
n_regions = int(ace.xss[idx])
n_energy_in = int(ace.xss[idx + 1 + 2*n_regions])
@ -586,10 +606,10 @@ class KalbachMann(AngleEnergy):
km_r.append(Tabulated1D(data[0], data[3]))
km_a.append(Tabulated1D(data[0], data[4]))
return cls(breakpoints, interpolation, energy, energy_out, km_r, km_a)
return cls(breakpoints, interpolation, energy, energy_out, km_r, km_a, is_photon=is_photon)
@classmethod
def from_endf(cls, file_obj, za_emitted, za_target, projectile_mass):
def from_endf(cls, file_obj, za_emitted, za_target, za_projectile):
"""Generate Kalbach-Mann distribution from an ENDF evaluation.
If the projectile is a neutron, the slope is calculated when it is
@ -606,14 +626,8 @@ class KalbachMann(AngleEnergy):
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.
za_projectile : int
ZA identifier of the projectile
Returns
-------
@ -621,6 +635,7 @@ class KalbachMann(AngleEnergy):
Kalbach-Mann energy-angle distribution
"""
is_photon = bool(za_projectile==0)
params, tab2 = get_tab2_record(file_obj)
lep = params[3]
ne = params[5]
@ -654,31 +669,19 @@ class KalbachMann(AngleEnergy):
a_i = values[:, 3]
calculated_slope.append(False)
else:
# 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)
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))
km_distribution = cls(tab2.breakpoints, tab2.interpolation, energy,
energy_out, precompound, slope)
energy_out, precompound, slope, is_photon)
# List of bool to indicate slope calculation by OpenMC
km_distribution._calculated_slope = calculated_slope

View file

@ -37,7 +37,7 @@ class DataLibrary(list):
name : str
Name of material, e.g. 'Am241'
data_type : str
Name of data type, e.g. 'neutron', 'photon', 'wmp', or 'thermal'
Name of data type, e.g. 'neutron', 'photon', 'photonuclear', 'wmp', or 'thermal'
.. versionadded:: 0.12
@ -61,7 +61,7 @@ class DataLibrary(list):
name : str
Name of material, e.g. 'Am241'
data_type : str
Name of data type, e.g. 'neutron', 'photon', 'wmp', or 'thermal'
Name of data type, e.g. 'neutron', 'photon', 'photonuclear', 'wmp', or 'thermal'
"""
library = self.get_by_material(name, data_type)

View file

@ -317,6 +317,13 @@ acer / %%%%%%%%%%%%%%%%%%%%%%%% Write out in ACE format %%%%%%%%%%%%%%%%%%%%%%%%
222 64 {mt_elastic} {elastic_type} {data.nmix} {energy_max} {iwt}/
"""
_PHOTONUCLEAR_TEMPLATE_ACER = """
acer / %%%%%%%%%%%%%%%%%%%%%%%% Write out in ACE format %%%%%%%%%%%%%%%%%%%%%%%%
{nendf} {nacer_in} 0 {nace} {ndir}
5 0 1 .{ext} /
'{library}: {zsymam}'/
{mat}
"""
def run(commands, tapein, tapeout, input_filename=None, stdout=False,
njoy_exec='njoy'):
@ -773,3 +780,123 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None,
for temperature in temperatures:
(output_dir / f"ace_{temperature:.1f}").unlink()
(output_dir / f"xsdir_{temperature:.1f}").unlink()
def make_ace_photonuclear(filename, acer=True, xsdir=None,
output_dir=None, pendf=False, error=0.001,
evaluation=None, **kwargs):
"""Generate incident photonuclear ACE file from an ENDF file
File names can be passed to
``[acer, xsdir, pendf]``
to specify the exact output for the given module.
Otherwise, the files will be writen to the current directory
or directory specified by ``output_dir``. Default file
names mirror the variable names, e.g. ``xsdir`` output
will be written to a file named ``xsdir`` unless otherwise
specified.
Parameters
----------
filename : str
Path to ENDF file
acer : bool or str, optional
Flag indicating if acer should be run. If a string is give, write the
resulting ``ace`` file to this location. Path of ACE file to write.
Defaults to ``"ace"``
xsdir : str, optional
Path of xsdir file to write. Defaults to ``"xsdir"`` in the same
directory as ``acer``
output_dir : str, optional
Directory to write output for requested modules. If not provided
and at least one of ``[pendf, broadr, heatr, gaspr, purr, acer]``
is ``True``, then write output files to current directory. If given,
must be a path to a directory.
pendf : str, optional
Path of pendf file to write. If omitted, the pendf file is not saved.
error : float, optional
Fractional error tolerance for NJOY processing
evaluation : openmc.data.endf.Evaluation, optional
If the ENDF file contains multiple material evaluations, this argument
indicates which evaluation should be used.
**kwargs
Keyword arguments passed to :func:`openmc.data.njoy.run`
Raises
------
subprocess.CalledProcessError
If the NJOY process returns with a non-zero status
IOError
If ``output_dir`` does not point to a directory
"""
if output_dir is None:
output_dir = Path()
else:
output_dir = Path(output_dir)
if not output_dir.is_dir():
raise IOError(f"{output_dir} is not a directory")
ev = evaluation if evaluation is not None else endf.Evaluation(filename)
mat = ev.material
zsymam = ev.target['zsymam']
# Determine name of library
library = '{}-{}.{}'.format(*ev.info['library'])
# Create njoy commands by modules
commands = ""
nendf, npendf = 20, 21
tapein = {nendf: filename}
tapeout = {}
if pendf:
tapeout[npendf] = (output_dir / "pendf") if pendf is True else pendf
# reconr
commands += _TEMPLATE_RECONR
nlast = npendf
commands = commands.format(**locals())
# acer
if acer:
nacer_in = nlast
# Extend input with an ACER run
nace = nacer_in + 1
ndir = nace + 1
ext = f'{1:02}'
commands += _PHOTONUCLEAR_TEMPLATE_ACER.format(**locals())
# Indicate tapes to save for each ACER run
tapeout[nace] = output_dir / f"ace_0.0"
tapeout[ndir] = output_dir / f"xsdir_0.0"
commands += 'stop\n'
run(commands, tapein, tapeout, **kwargs)
if acer:
ace = (output_dir / "ace") if acer is True else Path(acer)
xsdir = (ace.parent / "xsdir") if xsdir is None else xsdir
with ace.open('w') as ace_file, xsdir.open('w') as xsdir_file:
# Get contents of ACE file
text = (output_dir / f"ace_0.0").read_text()
# If the target is metastable, make sure that ZAID in the ACE
# file reflects this by adding 400
if ev.target['isomeric_state'] > 0:
mass_first_digit = int(text[3])
if mass_first_digit <= 2:
text = text[:3] + str(mass_first_digit + 4) + text[4:]
# Concatenate into destination ACE file
ace_file.write(text)
# Concatenate into destination xsdir file
xsdir_in = output_dir / f"xsdir_0.0"
xsdir_file.write(xsdir_in.read_text())
# Remove ACE/xsdir files
(output_dir / f"ace_0.0").unlink()
(output_dir / f"xsdir_0.0").unlink()

905
openmc/data/photonuclear.py Normal file
View file

@ -0,0 +1,905 @@
from collections import OrderedDict
from collections.abc import Callable, Iterable
from io import StringIO
from numbers import Integral, Real
from warnings import warn
import os
import tempfile
import h5py
import numpy as np
from openmc.mixin import EqualityMixin
from openmc.stats import Uniform
import openmc.checkvalue as cv
from . import HDF5_VERSION, HDF5_VERSION_MAJOR
from .ace import Table, get_metadata, get_table, Library
from .angle_distribution import AngleDistribution
from .angle_energy import AngleEnergy
from .correlated import CorrelatedAngleEnergy
from .data import ATOMIC_SYMBOL, EV_PER_MEV
from .endf import Evaluation, SUM_RULES, get_head_record, get_tab1_record, get_cont_record
from .fission_energy import FissionEnergyRelease
from .function import Tabulated1D
from .njoy import make_ace_photonuclear
from .reaction import Reaction, REACTION_NAME as _REACTION_NAME, FISSION_MTS, _get_products, _get_fission_products_endf, _get_photon_products_endf, _get_activation_products
from .product import Product
from .energy_distribution import EnergyDistribution, LevelInelastic, \
DiscretePhoton
from .function import Tabulated1D, Polynomial
from .kalbach_mann import KalbachMann
from .laboratory import LaboratoryAngleEnergy
from .nbody import NBodyPhaseSpace
from .product import Product
from .uncorrelated import UncorrelatedAngleEnergy
REACTION_NAME = {50 : '(gamma,n0)'}
REACTION_NAME.update({key:value.replace("(n,","(gamma,") for key,value in _REACTION_NAME.items()})
class PhotonuclearReaction(EqualityMixin):
def __init__(self, mt):
self._center_of_mass = True
self._redundant = False
self._q_value = 0.
self._xs = None
self._products = []
self._derived_products = []
self.mt = mt
def __repr__(self):
if self.mt in _REACTION_NAME:
return f"<PhotonuclearReaction: MT={self.mt} {REACTION_NAME[self.mt]}>"
else:
return f"<PhotonuclearReaction: MT={self.mt}>"
@property
def center_of_mass(self):
return self._center_of_mass
@center_of_mass.setter
def center_of_mass(self, center_of_mass):
cv.check_type('center of mass', center_of_mass, (bool, np.bool_))
self._center_of_mass = center_of_mass
@property
def redundant(self):
return self._redundant
@redundant.setter
def redundant(self, redundant):
cv.check_type('redundant', redundant, (bool, np.bool_))
self._redundant = redundant
@property
def q_value(self):
return self._q_value
@q_value.setter
def q_value(self, q_value):
cv.check_type('Q value', q_value, Real)
self._q_value = q_value
@property
def products(self):
return self._products
@products.setter
def products(self, products):
cv.check_type('reaction products', products, Iterable, Product)
self._products = products
@property
def derived_products(self):
return self._derived_products
@derived_products.setter
def derived_products(self, derived_products):
cv.check_type('reaction derived products', derived_products,
Iterable, Product)
self._derived_products = derived_products
@property
def xs(self):
return self._xs
@xs.setter
def xs(self, xs):
cv.check_type("reaction cross section", xs, Callable)
self._xs = xs
@classmethod
def from_endf(cls, ev, mt):
"""Generate a reaction from an ENDF evaluation
Parameters
----------
ev : openmc.data.endf.Evaluation
ENDF evaluation
mt : int
The MT value of the reaction to get data for
Returns
-------
rx : openmc.data.PhotonuclearReaction
Reaction data
"""
rx = cls(mt)
# Integrated cross section
if (3, mt) in ev.section:
file_obj = StringIO(ev.section[3, mt])
get_head_record(file_obj)
params, rx.xs = get_tab1_record(file_obj)
rx.q_value = params[1]
# Get fission product yields (nu) as well as delayed neutron energy
# distributions
if mt in FISSION_MTS:
rx.products, rx.derived_products = _get_fission_products_endf(ev)
if (6, mt) in ev.section:
# Product angle-energy distribution
for product in _get_products(ev, mt):
if mt in FISSION_MTS and product.particle == 'neutron':
rx.products[0].applicability = product.applicability
rx.products[0].distribution = product.distribution
else:
rx.products.append(product)
elif (4, mt) in ev.section or (5, mt) in ev.section:
# Uncorrelated angle-energy distribution
# Note that the energy distribution for MT=455 is read in
# _get_fission_products_endf rather than here
if (5, mt) in ev.section:
product = Product('photon')
file_obj = StringIO(ev.section[5, mt])
items = get_head_record(file_obj)
nk = items[4]
for i in range(nk):
params, applicability = get_tab1_record(file_obj)
dist = UncorrelatedAngleEnergy()
dist.energy = EnergyDistribution.from_endf(file_obj, params)
product.applicability.append(applicability)
product.distribution.append(dist)
elif mt == 2:
# Elastic scattering -- no energy distribution is given since it
# can be calulcated analytically
product = Product('photon')
dist = UncorrelatedAngleEnergy()
product.distribution.append(dist)
elif mt >= 50 and mt < 91:
# Level inelastic scattering -- no energy distribution is given
# since it can be calculated analytically. Here we determine the
# necessary parameters to create a LevelInelastic object
product = Product('neutron')
dist = UncorrelatedAngleEnergy()
A = ev.target['mass']
threshold = abs(rx.q_value)
mass_ratio = (A-1)/A
dist.energy = LevelInelastic(threshold, mass_ratio)
product.distribution.append(dist)
if (4, mt) in ev.section:
for dist in product.distribution:
dist.angle = AngleDistribution.from_endf(ev, mt)
if mt in FISSION_MTS and (5, mt) in ev.section:
# For fission reactions,
rx.products[0].applicability = product.applicability
rx.products[0].distribution = product.distribution
else:
rx.products.append(product)
if (8, mt) in ev.section:
rx.products += _get_activation_products(ev, rx)
if (12, mt) in ev.section or (13, mt) in ev.section:
rx.products += _get_photon_products_endf(ev, rx)
return rx
def to_hdf5(self, group):
"""Write reaction to an HDF5 group
Parameters
----------
group : h5py.Group
HDF5 group to write to
"""
group.attrs['mt'] = self.mt
if self.mt in REACTION_NAME:
group.attrs['label'] = np.bytes_(REACTION_NAME[self.mt])
else:
group.attrs['label'] = np.bytes_(self.mt)
group.attrs['Q_value'] = self.q_value
group.attrs['center_of_mass'] = 1 if self.center_of_mass else 0
group.attrs['redundant'] = 1 if self.redundant else 0
dset = group.create_dataset('xs', data=self.xs.y)
threshold_idx = getattr(self.xs, '_threshold_idx', 0)
dset.attrs['threshold_idx'] = threshold_idx
for i, p in enumerate(self.products):
pgroup = group.create_group(f'product_{i}')
p.to_hdf5(pgroup)
@classmethod
def from_hdf5(cls, group, energy):
"""Generate reaction from an HDF5 group
Parameters
----------
group : h5py.Group
HDF5 group to read from
energy : array
array of energies at which cross sections are tabulated at.
Returns
-------
openmc.data.Reaction
Reaction data
"""
mt = group.attrs['mt']
rx = cls(mt)
rx.q_value = group.attrs['Q_value']
rx.center_of_mass = bool(group.attrs['center_of_mass'])
rx.redundant = bool(group.attrs.get('redundant', False))
xs = group['xs'][()]
threshold_idx = group['xs'].attrs['threshold_idx']
tabulated_xs = Tabulated1D(energy[threshold_idx:], xs)
tabulated_xs._threshold_idx = threshold_idx
rx.xs = tabulated_xs
# Determine number of products
n_product = 0
for name in group:
if name.startswith('product_'):
n_product += 1
# Read reaction products
for i in range(n_product):
pgroup = group[f'product_{i}']
rx.products.append(Product.from_hdf5(pgroup))
return rx
@classmethod
def from_ace(cls, ace, i_reaction):
# Get nuclide energy grid
n_grid = ace.nxs[3]
grid = ace.xss[ace.jxs[1] : ace.jxs[1] + n_grid] * EV_PER_MEV
if i_reaction > 0:
mt = int(ace.xss[ace.jxs[6] + i_reaction - 1])
rx = cls(mt)
# Get Q-value of reaction
rx.q_value = ace.xss[ace.jxs[7] + i_reaction - 1]*EV_PER_MEV
# ==================================================================
# CROSS SECTION
# Get locator for cross-section data
loc = int(ace.xss[ace.jxs[8] + i_reaction - 1])
# Determine starting index on energy grid
threshold_idx = int(ace.xss[ace.jxs[9] + loc - 1]) - 1
# Determine number of energies in reaction
n_energy = int(ace.xss[ace.jxs[9] + loc])
energy = grid[threshold_idx : threshold_idx + n_energy]
# Read reaction cross section
xs = ace.xss[ace.jxs[9] + loc + 1 : ace.jxs[9] + loc + 1 + n_energy]
# Fix negatives -- known issue for Y89 in JEFF 3.2
if np.any(xs < 0.0):
warn(
"Negative cross sections found for MT={} in {}. Setting "
"to zero.".format(rx.mt, ace.name)
)
xs[xs < 0.0] = 0.0
tabulated_xs = Tabulated1D(energy, xs)
tabulated_xs._threshold_idx = threshold_idx
rx.xs = tabulated_xs
# ==================================================================
# YIELD AND ANGLE-ENERGY DISTRIBUTION
for i_typ in range(ace.nxs[5]):
loc = ace.jxs[10]+i_typ*ace.nxs[7]-1
mts = ace.xss[int(ace.xss[loc+5]):int(ace.xss[loc+5])+int(ace.xss[loc+2])].astype(int)
try:
i_mtr = mts.tolist().index(mt)
except ValueError:
#skip products for other reactions
continue
match int(ace.xss[loc+1]):
case 1:
particle=Product('neutron')
case 2:
particle=Product('photon')
case _:
#TODO: support more product particles
# for now skip particle if it is not a neutron or photon
warn(f"Unsupported secondary particle type in {ace.name} for MT={mt}. "
"This product will be skipped.")
continue
# Determine reference frame
ty = int(ace.xss[int(ace.xss[loc+6])+i_mtr])
rx.center_of_mass = (ty < 0)
# Determine multiplicity
idx = int(ace.xss[loc+8])+int(ace.xss[int(ace.xss[loc+7])+i_mtr])-1
match int(ace.xss[idx]):
case 6 | 16 | 12:
assert int(ace.xss[idx+1]) == mt
yield_ = Tabulated1D.from_ace(ace, idx+2)
case _:
raise NotImplementedError('partial yields not implemented yet')
particle.yield_ = yield_
# Determine locator for energy distribution
idx = int(ace.xss[int(ace.xss[loc+11])+i_mtr])
distribution = AngleEnergy.from_ace(ace, int(ace.xss[loc+12]), idx)
# Determine locator for angular distribution
idx = int(ace.xss[int(ace.xss[loc+9])+i_mtr])
if idx==0:
# No angular distribution data are given for this reaction,
# isotropic scattering is assumed in LAB
energy = np.array([particle.yield_.x[0], particle.yield_.x[-1]])
mu_isotropic = Uniform(-1., 1.)
distribution.angle = AngleDistribution(
energy, [mu_isotropic, mu_isotropic])
elif idx==-1:
pass
else:
assert idx>=0
distribution.angle = AngleDistribution.from_ace(ace, int(ace.xss[loc+10]), idx)
particle.distribution.append(distribution)
rx.products.append(particle)
else:
# elastic cross section
mt = 2
rx = cls(mt)
# Get elastic cross section values
elastic_xs = ace.xss[ace.jxs[4] : ace.jxs[4] + ace.nxs[3]]
# Fix negatives -- known issue for Ti46,49,50 in JEFF 3.2
if np.any(elastic_xs < 0.0):
warn(
"Negative elastic scattering cross section found for {}. "
"Setting to zero.".format(ace.name)
)
elastic_xs[elastic_xs < 0.0] = 0.0
tabulated_xs = Tabulated1D(grid, elastic_xs)
tabulated_xs._threshold_idx = 0
rx.xs = tabulated_xs
return rx
class IncidentPhotonuclear(EqualityMixin):
"""photo-nuclear interaction data.
This class stores data derived from an ENDF-6 format photo-nuclear interaction
sublibrary. Instances of this class are not normally instantiated by the
user but rather created using the factory methods
:meth:`Photonuclear.from_hdf5`, :meth:`Photonuclear.from_ace`, and
:meth:`Photonuclear.from_endf`.
Parameters
----------
name : str
Name of the nuclide using the GND naming convention
atomic_number : int
Number of photo-nuclears in the target nucleus
atomic_number : int
Number of photo-nuclears in the target nucleus
mass_number : int
Number of nucleons in the target nucleus
metastable : int
Metastable state of the target nucleus. A value of zero indicates ground
state.
atomic_weight_ratio : float
Atomic weight ratio of the target nuclide.
Attributes
----------
atomic_number : int
Number of photo-nuclears in the target nucleus
atomic_symbol : str
Atomic symbol of the nuclide, e.g., 'Zr'
atomic_weight_ratio : float
Atomic weight ratio of the target nuclide.
fission_energy : None or openmc.data.FissionEnergyRelease
The energy released by fission, tabulated by component (e.g. prompt
neutrons or beta particles) and dependent on incident neutron energy
mass_number : int
Number of nucleons in the target nucleus
metastable : int
Metastable state of the target nucleus. A value of zero indicates ground
state.
name : str
Name of the nuclide using the GND naming convention
reactions : collections.OrderedDict
Contains the cross sections, secondary angle and energy distributions,
and other associated data for each reaction. The keys are the MT values
and the values are Reaction objects.
"""
def __init__(
self, name, atomic_number, mass_number, metastable, atomic_weight_ratio
):
self.name = name
self.atomic_number = atomic_number
self.mass_number = mass_number
self.reactions = OrderedDict()
self.energy = []
self._fission_energy = None
self.metastable = metastable
self.atomic_weight_ratio = atomic_weight_ratio
def __contains__(self, mt):
return mt in self.reactions
def __getitem__(self, mt):
if mt in self.reactions:
return self.reactions[mt]
else:
# Try to create a redundant cross section
mts = self.get_reaction_components(mt)
if len(mts) > 0:
return self._get_redundant_reaction(mt, mts)
else:
raise KeyError(f'No reaction with MT={mt}.')
def __repr__(self):
return "<IncidentPhotonuclear: {}>".format(self.name)
def __iter__(self):
return iter(self.reactions.values())
@property
def atomic_number(self):
return self._atomic_number
@property
def name(self):
return self._name
@property
def mass_number(self):
return self._mass_number
@property
def metastable(self):
return self._metastable
@property
def atomic_weight_ratio(self):
return self._atomic_weight_ratio
@property
def atomic_symbol(self):
return ATOMIC_SYMBOL[self.atomic_number]
@name.setter
def name(self, name):
cv.check_type("name", name, str)
self._name = name
@atomic_number.setter
def atomic_number(self, atomic_number):
cv.check_type("atomic number", atomic_number, Integral)
cv.check_greater_than("atomic number", atomic_number, 0, True)
self._atomic_number = atomic_number
@mass_number.setter
def mass_number(self, mass_number):
cv.check_type("mass number", mass_number, Integral)
cv.check_greater_than("mass number", mass_number, 0, True)
self._mass_number = mass_number
@metastable.setter
def metastable(self, metastable):
cv.check_type("metastable", metastable, Integral)
cv.check_greater_than("metastable", metastable, 0, True)
self._metastable = metastable
@atomic_weight_ratio.setter
def atomic_weight_ratio(self, atomic_weight_ratio):
cv.check_type("atomic weight ratio", atomic_weight_ratio, Real)
cv.check_greater_than("atomic weight ratio", atomic_weight_ratio, 0.0)
self._atomic_weight_ratio = atomic_weight_ratio
@property
def fission_energy(self):
return self._fission_energy
@fission_energy.setter
def fission_energy(self, fission_energy):
cv.check_type('fission energy release', fission_energy,
FissionEnergyRelease)
self._fission_energy = fission_energy
@classmethod
def from_endf(cls, ev_or_filename):
"""Generate incident photo-nuclear data from an ENDF evaluation
Parameters
----------
ev_or_filename : openmc.data.endf.Evaluation or str
ENDF evaluation to read from. If given as a string, it is assumed to
be the filename for the ENDF file.
Returns
-------
openmc.data.Photonuclear
photo-nuclear interaction data
"""
if isinstance(ev_or_filename, Evaluation):
ev = ev_or_filename
else:
ev = Evaluation(ev_or_filename)
atomic_number = ev.target["atomic_number"]
mass_number = ev.target["mass_number"]
metastable = ev.target["isomeric_state"]
atomic_weight_ratio = ev.target["mass"]
element = ATOMIC_SYMBOL[atomic_number]
if metastable > 0:
name = f'{element}{mass_number}_m{metastable}'
else:
name = f'{element}{mass_number}'
data = cls(name, atomic_number, mass_number, metastable, atomic_weight_ratio)
# Read each reaction
for mf, mt, nc, mod in ev.reaction_list:
if mf == 3:
data.reactions[mt] = PhotonuclearReaction.from_endf(ev, mt)
# Read fission energy release (requires that we already know nu for
# fission)
if (1, 458) in ev.section:
data.fission_energy = FissionEnergyRelease.from_endf(ev, data)
data._evaluation = ev
return data
@classmethod
def from_ace(cls, ace_or_filename, metastable_scheme="nndc"):
"""Generate incident photo-nuclear continuous-energy data from an ACE table
Parameters
----------
ace_or_filename : openmc.data.ace.Table or str
ACE table to read from. If the value is a string, it is assumed to
be the filename for the ACE file.
metastable_scheme : {'nndc', 'mcnp'}
Determine how ZAID identifiers are to be interpreted in the case of
a metastable nuclide. Because the normal ZAID (=1000*Z + A) does not
encode metastable information, different conventions are used among
different libraries. In MCNP libraries, the convention is to add 400
for a metastable nuclide except for Am242m, for which 95242 is
metastable and 95642 (or 1095242 in newer libraries) is the ground
state. For NNDC libraries, ZAID is given as 1000*Z + A + 100*m.
Returns
-------
openmc.data.Photonuclear
Incident photo-nuclear continuous-energy data
"""
# First obtain the data for the first provided ACE table/file
if isinstance(ace_or_filename, Table):
ace = ace_or_filename
else:
ace = get_table(ace_or_filename)
# If mass number hasn't been specified, make an educated guess
zaid, xs = ace.name.split(".")
if not xs.endswith("u"):
raise TypeError(
"{} is not a continuous-energy photo-nuclear ACE table.".format(ace)
)
name, element, Z, mass_number, metastable = get_metadata(
int(zaid), metastable_scheme
)
data = cls(name, Z, mass_number, metastable, ace.atomic_weight_ratio)
# Read energy grid
n_energy = ace.nxs[3]
i = ace.jxs[1]
energy = ace.xss[i : i + n_energy] * EV_PER_MEV
data.energy = energy
total_xs = ace.xss[ace.jxs[2] : ace.jxs[2] + n_energy]
nonelastic_xs = ace.xss[ace.jxs[3] : ace.jxs[3] + n_energy]
elastic_xs = total_xs-nonelastic_xs
heating_number = ace.xss[ace.jxs[5] : ace.jxs[5] + n_energy] * EV_PER_MEV
# Create redundant reaction for total (MT=1)
total = PhotonuclearReaction(1)
total.xs = Tabulated1D(energy, total_xs)
total.redundant = True
data.reactions[1] = total
# Create redundant reaction for nonelastic (MT=3)
nonelastic = PhotonuclearReaction(3)
nonelastic.xs = Tabulated1D(energy, nonelastic_xs)
nonelastic.redundant = True
data.reactions[3] = nonelastic
# Create redundant reaction for heating (MT=301)
heating = PhotonuclearReaction(301)
heating.xs = Tabulated1D(energy, heating_number*total_xs)
heating.redundant = True
data.reactions[301] = heating
# Read each reaction
n_reaction = ace.nxs[4] + 1
for i in range(n_reaction):
if i==0:
if ace.jxs[4]==0:
assert np.allclose(total_xs,nonelastic_xs)
else:
raise NotImplementedError('photonuclear elastic scattering is not supported.')
else:
rx = PhotonuclearReaction.from_ace(ace, i)
data.reactions[rx.mt] = rx
return data
def export_to_hdf5(self, path, mode="a", libver="earliest"):
"""Export incident photo-nuclear data to an HDF5 file.
Parameters
----------
path : str
Path to write HDF5 file to
mode : {'r', 'r+', 'w', 'x', 'a'}
Mode that is used to open the HDF5 file. This is the second argument
to the :class:`h5py.File` constructor.
libver : {'earliest', 'latest'}
Compatibility mode for the HDF5 file. 'latest' will produce files
that are less backwards compatible but have performance benefits.
"""
# Open file and write version
f = h5py.File(str(path), mode, libver=libver)
f.attrs["filetype"] = np.bytes_("data_photonuclear")
if "version" not in f.attrs:
f.attrs["version"] = np.array(HDF5_VERSION)
# If data come from ENDF, don't allow exporting to HDF5
if hasattr(self, '_evaluation'):
raise NotImplementedError('Cannot export incident photonuclear data that '
'originated from an ENDF file.')
# Write basic data
g = f.create_group(self.name)
g.attrs['Z'] = self.atomic_number
g.attrs['A'] = self.mass_number
g.attrs['metastable'] = self.metastable
g.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio
# Determine union energy grid
union_grid = np.array([])
for rx in self:
union_grid = np.union1d(union_grid, rx.xs.x)
for product in rx.products:
union_grid = np.union1d(union_grid, product.yield_.x)
g.create_dataset("energy", data=union_grid)
# Update reactions according to union energy grid
for rx in self:
if tuple(rx.xs.interpolation) != (2,):
raise NotImplementedError('Only linear-linear interpolable reactions are supported.')
threshold_idx, = np.flatnonzero(union_grid == rx.xs.x[0])
xs = rx.xs(union_grid[threshold_idx:])
tab1d = Tabulated1D(union_grid[threshold_idx:], xs)
tab1d._threshold_idx = threshold_idx
rx.xs = tab1d
# Write reaction data
rxs_group = g.create_group('reactions')
for rx in self.reactions.values():
# Skip writing redundant reaction if it doesn't have neutron
# production or is a summed transmutation reaction.
# Also write heating.
if rx.redundant:
neutron_rx = any(p.particle == 'neutron' for p in rx.products)
keep_mts = (301,)
if not (neutron_rx or rx.mt in keep_mts):
continue
rx_group = rxs_group.create_group(f'reaction_{rx.mt:03}')
rx.to_hdf5(rx_group)
# Write fission energy release data
if self.fission_energy is not None:
fer_group = g.create_group('fission_energy_release')
self.fission_energy.to_hdf5(fer_group)
f.close()
@classmethod
def from_hdf5(cls, group_or_filename):
"""Generate photo-nuclear interaction data from HDF5 group
Parameters
----------
group_or_filename : h5py.Group or str
HDF5 group containing interaction data. If given as a string, it is
assumed to be the filename for the HDF5 file, and the first group is
used to read from.
Returns
-------
openmc.data.Photonuclear
photo-nuclear interaction data
"""
if isinstance(group_or_filename, h5py.Group):
group = group_or_filename
else:
h5file = h5py.File(str(group_or_filename), "r")
# Make sure version matches
if "version" in h5file.attrs:
major, minor = h5file.attrs["version"]
# For now all versions of HDF5 data can be read
else:
raise IOError(
"HDF5 data does not indicate a version. Your installation of "
"the OpenMC Python API expects version {}.x data.".format(
HDF5_VERSION_MAJOR
)
)
group = list(h5file.values())[0]
name = group.name[1:]
atomic_number = group.attrs["Z"]
mass_number = group.attrs["A"]
metastable = group.attrs["metastable"]
atomic_weight_ratio = group.attrs["atomic_weight_ratio"]
data = cls(name, atomic_number, mass_number, metastable, atomic_weight_ratio)
# Read energy grid
data.energy = group["energy"][()]
# Read reaction data
rxs_group = group["reactions"]
for name, obj in sorted(rxs_group.items()):
if name.startswith("reaction_"):
rx = PhotonuclearReaction.from_hdf5(obj, data.energy)
data.reactions[rx.mt] = rx
# Read fission energy release data
if 'fission_energy_release' in group:
fer_group = group['fission_energy_release']
data.fission_energy = FissionEnergyRelease.from_hdf5(fer_group)
# Close h5file when done if was opened
if not isinstance(group_or_filename, h5py.Group):
h5file.close()
return data
@classmethod
def from_njoy(cls, filename, evaluation=None, **kwargs):
"""Generate incident photo-nuclear data by running NJOY.
Parameters
----------
filename : str
Path to ENDF file
evaluation : openmc.data.endf.Evaluation, optional
If the ENDF file contains multiple material evaluations, this
argument indicates which evaluation to use.
**kwargs
Keyword arguments passed to :func:`openmc.data.njoy.make_ace_photonuclear`
Returns
-------
data : openmc.data.Photonuclear
Incident photo-nuclear continuous-energy data
"""
with tempfile.TemporaryDirectory() as tmpdir:
# Run NJOY to create an ACE library
kwargs.setdefault("output_dir", tmpdir)
for key in ("acer", "pendf"):
kwargs.setdefault(key, os.path.join(kwargs["output_dir"], key))
kwargs["evaluation"] = evaluation
make_ace_photonuclear(filename, **kwargs)
# Create instance from ACE tables within library
lib = Library(kwargs["acer"])
data = cls.from_ace(lib.tables[0])
# Add fission energy release data
ev = evaluation if evaluation is not None else Evaluation(filename)
if (1, 458) in ev.section:
data.fission_energy = FissionEnergyRelease.from_endf(ev, data)
return data
def get_reaction_components(self, mt):
"""Determine what reactions make up redundant reaction.
Parameters
----------
mt : int
ENDF MT number of the reaction to find components of.
Returns
-------
mts : list of int
ENDF MT numbers of reactions that make up the redundant reaction and
have cross sections provided.
"""
mts = []
if mt in SUM_RULES:
for mt_i in SUM_RULES[mt]:
mts += self.get_reaction_components(mt_i)
if mts:
return mts
else:
return [mt] if mt in self else []
def _get_redundant_reaction(self, mt, mts):
"""Create redundant reaction from its components
Parameters
----------
mt : int
MT value of the desired reaction
mts : iterable of int
MT values of its components
Returns
-------
openmc.Reaction
Redundant reaction
"""
rx = PhotonuclearReaction(mt)
# Get energy grid
energy = self.energy
xss = [self.reactions[mt_i].xs for mt_i in mts]
idx = min([xs._threshold_idx if hasattr(xs, '_threshold_idx')
else 0 for xs in xss])
rx.xs = Tabulated1D(energy[idx:], Sum(xss)(energy[idx:]))
rx.xs._threshold_idx = idx
rx.redundant = True
return rx

View file

@ -99,6 +99,8 @@ def _get_products(ev, mt):
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.
NotImplementedError
When the projectile is not a neutron and not a photon.
Returns
-------
@ -175,11 +177,16 @@ def _get_products(ev, mt):
)
zat = ev.target["atomic_number"] * 1000 + ev.target["mass_number"]
projectile_mass = ev.projectile["mass"]
if ev.projectile['mass'] == 0.0:
projectile_za = 0
elif np.isclose(ev.projectile['mass'], 1.0, atol=1.0e-12, rtol=0.):
projectile_za = 1
else:
raise NotImplementedError('Unknown projectile')
p.distribution = [KalbachMann.from_endf(file_obj,
za,
zat,
projectile_mass)]
projectile_za)]
elif law == 2:
# Discrete two-body scattering

View file

@ -143,6 +143,8 @@ class Element(str):
tree = ET.parse(cross_sections)
root = tree.getroot()
for child in root.findall('library'):
if child.attrib['type'] == 'photonuclear':
continue
nuclide = child.attrib['materials']
if re.match(r'{}\d+'.format(self), nuclide):
library_nuclides.add(nuclide)

View file

@ -53,6 +53,8 @@ Library::Library(pugi::xml_node node, const std::string& directory)
type_ = Type::thermal;
} else if (type == "photon") {
type_ = Type::photon;
} else if (type == "photonuclear") {
type_ = Type::photonuclear;
} else if (type == "wmp") {
type_ = Type::wmp;
} else {

View file

@ -103,16 +103,14 @@ def test_kalbach_slope():
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
)
kalbach_slope(
energy_projectile=energy_projectile,
energy_emitted=energy_emitted,
za_projectile=0,
za_emitted=1,
za_target=6012
)
assert kalbach_slope(
energy_projectile=energy_projectile,

View file

@ -0,0 +1,133 @@
from collections.abc import Mapping, Callable
import os
import numpy as np
import pandas as pd
import pytest
import openmc.data
import openmc
from . import needs_njoy
@pytest.fixture(scope='module')
def pu239():
"""Pu239 HDF5 data."""
directory = os.path.dirname(openmc.config.get('cross_sections'))
filename = os.path.join(directory, 'photonuclear', 'Pu239.h5')
return openmc.data.IncidentPhotonuclear.from_hdf5(filename)
@pytest.fixture(scope='module')
def u235(endf_data):
endf_file = os.path.join(endf_data, 'gammas', 'g-092_U_235.endf')
return openmc.data.IncidentPhotonuclear.from_njoy(endf_file)
@pytest.fixture(scope='module')
def be9(endf_data):
"""Be9 ENDF data (contains laboratory angle-energy distribution)."""
filename = os.path.join(endf_data, 'gammas', 'g-004_Be_009.endf')
return openmc.data.IncidentPhotonuclear.from_endf(filename)
@pytest.fixture(scope='module')
def h2(endf_data):
endf_file = os.path.join(endf_data, 'gammas', 'g-001_H_002.endf')
return openmc.data.IncidentPhotonuclear.from_njoy(endf_file)
def test_attributes(pu239):
assert pu239.name == 'Pu239'
assert pu239.mass_number == 239
assert pu239.metastable == 0
assert pu239.atomic_symbol == 'Pu'
assert pu239.atomic_weight_ratio == pytest.approx(236.9986)
@needs_njoy
def test_fission_energy(u235):
fer = u235.fission_energy
assert isinstance(fer, openmc.data.FissionEnergyRelease)
components = ['betas', 'delayed_neutrons', 'delayed_photons', 'fragments',
'neutrinos', 'prompt_neutrons', 'prompt_photons', 'recoverable',
'total', 'q_prompt', 'q_recoverable', 'q_total']
for c in components:
assert isinstance(getattr(fer, c), Callable)
def test_energy_grid(pu239):
grid = pu239.energy
assert np.all(np.diff(grid) >= 0.0)
def test_reactions(pu239):
assert 18 in pu239.reactions
assert isinstance(pu239.reactions[18], openmc.data.PhotonuclearReaction)
with pytest.raises(KeyError):
pu239.reactions[2]
def test_fission(pu239):
fission = pu239.reactions[18]
assert not fission.center_of_mass
assert fission.q_value == pytest.approx(197380000.0)
assert fission.mt == 18
assert len(fission.products) == 1
prompt = fission.products[0]
assert prompt.particle == 'neutron'
assert prompt.yield_(1.0e-5) == pytest.approx(1.74559)
@needs_njoy
def test_kerma(run_in_tmpdir, h2):
assert 301 in h2
h2.export_to_hdf5("H2.h5")
read_in = openmc.data.IncidentPhotonuclear.from_hdf5("H2.h5")
assert 301 in read_in
assert np.all(read_in[301].xs.y == h2[301].xs.y)
@needs_njoy
def test_get_reaction_components(h2):
assert h2.get_reaction_components(1) == [50]
assert h2.get_reaction_components(51) == []
def test_export_to_hdf5(tmpdir, pu239, be9):
filename = str(tmpdir.join('pu239.h5'))
pu239.export_to_hdf5(filename)
assert os.path.exists(filename)
with pytest.raises(NotImplementedError):
be9.export_to_hdf5('be9.h5')
@needs_njoy
def test_ace_convert(endf_data, run_in_tmpdir):
filename = os.path.join(endf_data, 'gammas', 'g-001_H_002.endf')
ace_ascii = 'ace_ascii'
ace_binary = 'ace_binary'
openmc.data.njoy.make_ace_photonuclear(filename, acer=ace_ascii)
# Convert to binary
openmc.data.ace.ascii_to_binary(ace_ascii, ace_binary)
# Make sure conversion worked
lib_ascii = openmc.data.ace.Library(ace_ascii)
lib_binary = openmc.data.ace.Library(ace_binary)
for tab_a, tab_b in zip(lib_ascii.tables, lib_binary.tables):
assert tab_a.name == tab_b.name
assert tab_a.atomic_weight_ratio == pytest.approx(tab_b.atomic_weight_ratio)
assert tab_a.temperature == pytest.approx(tab_b.temperature)
assert np.all(tab_a.nxs == tab_b.nxs)
assert np.all(tab_a.jxs == tab_b.jxs)
assert tab_a.zaid == tab_b.zaid
assert tab_a.data_type == tab_b.data_type
def test_ace_table_types():
TT = openmc.data.ace.TableType
assert TT.from_suffix('u') == TT.PHOTONUCLEAR
assert TT.from_suffix('80u') == TT.PHOTONUCLEAR

View file

@ -1,13 +1,14 @@
#!/bin/bash
# touching this script invalidate the test data xs cache
set -ex
# Download HDF5 data
if [[ ! -e $HOME/nndc_hdf5/cross_sections.xml ]]; then
wget -q -O - https://anl.box.com/shared/static/teaup95cqv8s9nn56hfn7ku8mmelr95p.xz | tar -C $HOME -xJ
wget -q -O - https://github.com/GuySten/data/releases/download/v1.0.0/nndc_hdf5_test.tar.xz | tar -C $HOME -xJ
fi
# Download ENDF/B-VII.1 distribution
ENDF=$HOME/endf-b-vii.1
if [[ ! -d $ENDF/neutrons || ! -d $ENDF/photoat || ! -d $ENDF/atomic_relax ]]; then
if [[ ! -d $ENDF/neutrons || ! -d $ENDF/photoat || ! -d $ENDF/atomic_relax || ! -d $ENDF/gammas ]]; then
wget -q -O - https://anl.box.com/shared/static/4kd2gxnf4gtk4w1c8eua5fsua22kvgjb.xz | tar -C $HOME -xJ
fi