from collections.abc import Iterable
from collections import namedtuple
from difflib import get_close_matches
from numbers import Real
from io import StringIO
import itertools
import os
import re
import tempfile
from warnings import warn
import numpy as np
import h5py
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
from openmc.stats import Discrete, Tabular
from . import HDF5_VERSION, HDF5_VERSION_MAJOR, endf
from .data import K_BOLTZMANN, ATOMIC_SYMBOL, EV_PER_MEV, isotopes
from .ace import Table, get_table, Library
from .angle_energy import AngleEnergy
from .function import Tabulated1D, Function1D, Sum
from .njoy import make_ace_thermal
from .thermal_angle_energy import (CoherentElasticAE, IncoherentElasticAE,
IncoherentElasticAEDiscrete,
IncoherentInelasticAEDiscrete,
IncoherentInelasticAE, MixedElasticAE)
_THERMAL_NAMES = {
'c_Al27': ('al', 'al27', 'al-27', '13-al- 27'),
'c_Al_in_Al2O3': ('asap00', 'asap', 'al(al2o3)'),
'c_Be': ('be', 'be-metal', 'be-met', 'be00', 'be-metal', 'be metal', '4-be'),
'c_BeO': ('beo',),
'c_Be_in_BeO': ('bebeo', 'be-beo', 'be-o', 'be/o', 'bbeo00', 'be(beo)'),
'c_Be_in_Be2C': ('bebe2c',),
'c_Be_in_FLiBe': ('beflib', 'be(flibe)'),
'c_C6H6': ('benz', 'c6h6', 'benzine'),
'c_C_in_SiC': ('csic', 'c-sic', 'c(3c-sic)'),
'c_Ca_in_CaH2': ('cah', 'cah00', 'cacah2', 'ca(cah2)'),
'c_D_in_D2O': ('dd2o', 'd-d2o', 'hwtr', 'hw', 'dhw00', 'd(d2o)'),
'c_D_in_D2O_ice': ('dice',),
'c_F_in_FLiBe': ('fflibe', 'f(flibe)'),
'c_Fe56': ('fe', 'fe56', 'fe-56', '26-fe- 56'),
'c_Graphite': ('graph', 'grph', 'gr', 'gr00', 'graphite'),
'c_Graphite_10p': ('grph10', '10p graphit'),
'c_Graphite_30p': ('grph30', '30p graphit'),
'c_H_in_C5O2H8': ('lucite', 'c5o2h8', 'h-luci', 'h(lucite)'),
'c_H_in_CaH2': ('hcah2', 'hca00', 'h(cah2)'),
'c_H_in_CH2': ('hch2', 'poly', 'pol', 'h-poly', 'pol00', 'h(ch2)'),
'c_H_in_CH4_liquid': ('lch4', 'lmeth', 'l-ch4'),
'c_H_in_CH4_solid': ('sch4', 'smeth', 's-ch4'),
'c_H_in_CH4_solid_phase_II': ('sch4p2',),
'c_H_in_H2O': ('hh2o', 'h-h2o', 'lwtr', 'lw', 'lw00', 'h(h2o)'),
'c_H_in_H2O_solid': ('hice', 'h-ice', 'ice00', 'h(ice-ih)', 'h(ice)'),
'c_H_in_HF': ('hhf', 'h(hf)'),
'c_H_in_Mesitylene': ('mesi00', 'mesi', 'mesi-phii'),
'c_H_in_ParaffinicOil': ('hparaf', 'h(paraffin'),
'c_H_in_Toluene': ('tol00', 'tol', 'tolue-phii'),
'c_H_in_UH3': ('huh3', 'h(uh3)'),
'c_H_in_YH2': ('hyh2', 'h-yh2', 'h(yh2)'),
'c_H_in_ZrH': ('hzrh', 'h-zrh', 'h-zr', 'h/zr', 'hzr', 'hzr00', 'h(zrh)'),
'c_H_in_ZrH2': ('hzrh2', 'h(zrh2)'),
'c_H_in_ZrHx': ('hzrhx', 'h(zrhx)'),
'c_Li_in_FLiBe': ('liflib', 'li(flibe)'),
'c_Mg24': ('mg', 'mg24', 'mg00', '24-mg'),
'c_N_in_UN': ('n-un', 'n(un)', 'n(un) l'),
'c_O_in_Al2O3': ('osap00', 'osap', 'o(al2o3)'),
'c_O_in_BeO': ('obeo', 'o-beo', 'o-be', 'o/be', 'obeo00', 'o(beo)'),
'c_O_in_D2O': ('od2o', 'o-d2o', 'ohw00', 'o(d2o)'),
'c_O_in_H2O_solid': ('oice', 'o-ice', 'o(ice-ih)'),
'c_O_in_UO2': ('ouo2', 'o-uo2', 'o2-u', 'o2/u', 'ouo200', 'o(uo2)'),
'c_ortho_D': ('orthod', 'orthoD', 'dortho', 'od200', 'ortod', 'ortho-d'),
'c_ortho_H': ('orthoh', 'orthoH', 'hortho', 'oh200', 'ortoh', 'ortho-h'),
'c_para_D': ('parad', 'paraD', 'dpara', 'pd200', 'para-d'),
'c_para_H': ('parah', 'paraH', 'hpara', 'ph200', 'para-h'),
'c_Si28': ('si00', 'sili', 'si'),
'c_Si_in_SiC': ('sisic', 'si-sic', 'si(3c-sic)'),
'c_SiO2_alpha': ('sio2', 'sio2a', 'sio2alpha'),
'c_SiO2_beta': ('sio2b', 'sio2beta'),
'c_U_in_UN': ('u-un', 'u(un)', 'u(un) l'),
'c_U_in_UO2': ('uuo2', 'u-uo2', 'u-o2', 'u/o2', 'uuo200', 'u(uo2)'),
'c_Y_in_YH2': ('yyh2', 'y-yh2', 'y(yh2)'),
'c_Zr_in_ZrH': ('zrzrh', 'zr-zrh', 'zr-h', 'zr/h', 'zr(zrh)'),
'c_Zr_in_ZrH2': ('zrzrh2', 'zr(zrh2)'),
'c_Zr_in_ZrHx': ('zrzrhx', 'zr(zrhx)'),
}
def _temperature_str(T):
# round() normally returns an int when called with a single argument, but
# numpy floats overload rounding to return another float
return "{}K".format(int(round(T)))
def get_thermal_name(name):
"""Get proper S(a,b) table name, e.g. 'HH2O' -> 'c_H_in_H2O'
Parameters
----------
name : str
Name of an ACE thermal scattering table
Returns
-------
str
GNDS-format thermal scattering name
"""
if name in _THERMAL_NAMES:
return name
else:
for proper_name, names in _THERMAL_NAMES.items():
if name.lower() in names:
return proper_name
# Make an educated guess?? This actually works well for
# JEFF-3.2 which stupidly uses names like lw00.32t,
# lw01.32t, etc. for different temperatures
# First, construct a list of all the values/keys in the names
# dictionary
all_names = itertools.chain(_THERMAL_NAMES.keys(),
*_THERMAL_NAMES.values())
matches = get_close_matches(name, all_names, cutoff=0.5)
if matches:
# Figure out the key for the corresponding match
match = matches[0]
if match not in _THERMAL_NAMES:
for key, value_list in _THERMAL_NAMES.items():
if match in value_list:
match = key
break
warn('Thermal scattering material "{}" is not recognized. '
'Assigning a name of {}.'.format(name, match))
return match
else:
# OK, we give up. Just use the ACE name.
warn('Thermal scattering material "{0}" is not recognized. '
'Assigning a name of c_{0}.'.format(name))
return 'c_' + name
[docs]class CoherentElastic(Function1D):
r"""Coherent elastic scattering data from a crystalline material
The integrated cross section for coherent elastic scattering from a
powdered crystalline material may be represented as:
.. math::
\sigma(E,T) = \frac{1}{E} \sum\limits_{i=1}^{E_i < E} s_i(T)
where :math:`s_i(T)` is proportional the structure factor in [eV-b] at
the moderator temperature :math:`T` in Kelvin.
Parameters
----------
bragg_edges : Iterable of float
Bragg edge energies in eV
factors : Iterable of float
Partial sum of structure factors, :math:`\sum\limits_{i=1}^{E_i<E} s_i`
Attributes
----------
bragg_edges : Iterable of float
Bragg edge energies in eV
factors : Iterable of float
Partial sum of structure factors, :math:`\sum\limits_{i=1}^{E_i<E} s_i`
"""
def __init__(self, bragg_edges, factors):
self.bragg_edges = bragg_edges
self.factors = factors
def __call__(self, E):
idx = np.searchsorted(self.bragg_edges, E) - 1
if isinstance(E, Iterable):
E = np.asarray(E)
nonzero = idx >= 0
xs = np.zeros_like(E)
xs[nonzero] = self.factors[idx[nonzero]] / E[nonzero]
return xs
else:
return self.factors[idx] / E if idx >= 0 else 0.0
def __len__(self):
return len(self.bragg_edges)
@property
def bragg_edges(self):
return self._bragg_edges
@bragg_edges.setter
def bragg_edges(self, bragg_edges):
cv.check_type('Bragg edges', bragg_edges, Iterable, Real)
self._bragg_edges = np.asarray(bragg_edges)
@property
def factors(self):
return self._factors
@factors.setter
def factors(self, factors):
cv.check_type('structure factor cumulative sums', factors,
Iterable, Real)
self._factors = np.asarray(factors)
[docs] def to_hdf5(self, group, name):
"""Write coherent elastic scattering to an HDF5 group
Parameters
----------
group : h5py.Group
HDF5 group to write to
name : str
Name of the dataset to create
"""
dataset = group.create_dataset(name, data=np.vstack(
[self.bragg_edges, self.factors]))
dataset.attrs['type'] = np.bytes_(type(self).__name__)
[docs] @classmethod
def from_hdf5(cls, dataset):
"""Read coherent elastic scattering from an HDF5 dataset
Parameters
----------
dataset : h5py.Dataset
HDF5 dataset to read from
Returns
-------
openmc.data.CoherentElastic
Coherent elastic scattering cross section
"""
bragg_edges = dataset[0, :]
factors = dataset[1, :]
return cls(bragg_edges, factors)
[docs]class IncoherentElastic(Function1D):
r"""Incoherent elastic scattering cross section
Elastic scattering can be treated in the incoherent approximation for
partially ordered systems such as ZrHx and polyethylene. The integrated
cross section can be obtained as:
.. math::
\sigma(E,T) = \frac{\sigma_b}{2} \left ( \frac{1 - e^{-4EW'(T)}}
{2EW'(T)} \right )
where :math:`\sigma_b` is the characteristic bound cross section, and
:math:`W'(T)` is the Debye-Waller integral divided by the atomic mass
in [eV\ :math:`^{-1}`].
Parameters
----------
bound_xs : float
Characteristic bound cross section in [b]
debye_waller : float
Debye-Waller integral in [eV\ :math:`^{-1}`]
Attributes
----------
bound_xs : float
Characteristic bound cross section in [b]
debye_waller : float
Debye-Waller integral in [eV\ :math:`^{-1}`]
"""
def __init__(self, bound_xs, debye_waller):
self.bound_xs = bound_xs
self.debye_waller = debye_waller
def __call__(self, E):
W = self.debye_waller
return self.bound_xs / 2.0 * (1 - np.exp(-4*E*W)) / (2*E*W)
[docs] def to_hdf5(self, group, name):
"""Write incoherent elastic scattering to an HDF5 group
Parameters
----------
group : h5py.Group
HDF5 group to write to
name : str
Name of the dataset to create
"""
data = np.array([self.bound_xs, self.debye_waller])
dataset = group.create_dataset(name, data=data)
dataset.attrs['type'] = np.bytes_(type(self).__name__)
[docs] @classmethod
def from_hdf5(cls, dataset):
"""Read incoherent elastic scattering from an HDF5 dataset
Parameters
----------
dataset : h5py.Dataset
HDF5 dataset to read from
Returns
-------
openmc.data.IncoherentElastic
Incoherent elastic scattering cross section
"""
bound_xs, debye_waller = dataset[()]
return cls(bound_xs, debye_waller)
[docs]class ThermalScatteringReaction(EqualityMixin):
r"""Thermal scattering reaction
This class is used to hold the integral and differential cross sections
for either elastic or inelastic thermal scattering.
Parameters
----------
xs : dict of str to Function1D
Integral cross section at each temperature
distribution : dict of str to AngleEnergy
Secondary angle-energy distribution at each temperature
Attributes
----------
xs : dict of str to Function1D
Integral cross section at each temperature
distribution : dict of str to AngleEnergy
Secondary angle-energy distribution at each temperature
"""
def __init__(self, xs, distribution):
self.xs = xs
self.distribution = distribution
[docs] def to_hdf5(self, group, name):
"""Write thermal scattering reaction to HDF5
Parameters
----------
group : h5py.Group
HDF5 group to write to
name : {'elastic', 'inelastic'}
Name of reaction to write
"""
for T, xs in self.xs.items():
Tgroup = group.require_group(T)
rx_group = Tgroup.create_group(name)
xs.to_hdf5(rx_group, 'xs')
dgroup = rx_group.create_group('distribution')
self.distribution[T].to_hdf5(dgroup)
[docs] @classmethod
def from_hdf5(cls, group, name, temperatures):
"""Generate thermal scattering reaction data from HDF5
Parameters
----------
group : h5py.Group
HDF5 group to read from
name : {'elastic', 'inelastic'}
Name of the reaction to read
temperatures : Iterable of str
Temperatures to read
Returns
-------
openmc.data.ThermalScatteringReaction
Thermal scattering reaction data
"""
xs = {}
distribution = {}
for T in temperatures:
rx_group = group[T][name]
xs[T] = Function1D.from_hdf5(rx_group['xs'])
if isinstance(xs[T], CoherentElastic):
distribution[T] = CoherentElasticAE(xs[T])
else:
distribution[T] = AngleEnergy.from_hdf5(rx_group['distribution'])
return cls(xs, distribution)
[docs]class ThermalScattering(EqualityMixin):
"""A ThermalScattering object contains thermal scattering data as represented by
an S(alpha, beta) table.
Parameters
----------
name : str
Name of the material using GNDS convention, e.g. c_H_in_H2O
atomic_weight_ratio : float
Atomic mass ratio of the target nuclide.
kTs : Iterable of float
List of temperatures of the target nuclide in the data set.
The temperatures have units of eV.
Attributes
----------
atomic_weight_ratio : float
Atomic mass ratio of the target nuclide.
energy_max : float
Maximum energy for thermal scattering data in [eV]
elastic : openmc.data.ThermalScatteringReaction or None
Elastic scattering derived in the coherent or incoherent approximation
inelastic : openmc.data.ThermalScatteringReaction
Inelastic scattering cross section derived in the incoherent
approximation
name : str
Name of the material using GNDS convention, e.g. c_H_in_H2O
temperatures : Iterable of str
List of string representations the temperatures of the target nuclide
in the data set. The temperatures are strings of the temperature,
rounded to the nearest integer; e.g., '294K'
kTs : Iterable of float
List of temperatures of the target nuclide in the data set.
The temperatures have units of eV.
nuclides : Iterable of str
Nuclide names that the thermal scattering data applies to
"""
def __init__(self, name, atomic_weight_ratio, energy_max, kTs):
self.name = name
self.atomic_weight_ratio = atomic_weight_ratio
self.energy_max = energy_max
self.kTs = kTs
self.elastic = None
self.inelastic = None
self.nuclides = []
def __repr__(self):
if hasattr(self, 'name'):
return "<Thermal Scattering Data: {}>".format(self.name)
else:
return "<Thermal Scattering Data>"
@property
def temperatures(self):
return [_temperature_str(kT / K_BOLTZMANN) for kT in self.kTs]
[docs] def export_to_hdf5(self, path, mode='a', libver='earliest'):
"""Export table to an HDF5 file.
Parameters
----------
path : str
Path to write HDF5 file to
mode : {'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
with h5py.File(str(path), mode, libver=libver) as f:
f.attrs['filetype'] = np.bytes_('data_thermal')
f.attrs['version'] = np.array(HDF5_VERSION)
# Write basic data
g = f.create_group(self.name)
g.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio
g.attrs['energy_max'] = self.energy_max
g.attrs['nuclides'] = np.array(self.nuclides, dtype='S')
ktg = g.create_group('kTs')
for i, temperature in enumerate(self.temperatures):
ktg.create_dataset(temperature, data=self.kTs[i])
# Write elastic/inelastic reaction data
if self.elastic is not None:
self.elastic.to_hdf5(g, 'elastic')
self.inelastic.to_hdf5(g, 'inelastic')
[docs] def add_temperature_from_ace(self, ace_or_filename, name=None):
"""Add data to the ThermalScattering object from an ACE file at a
different temperature.
Parameters
----------
ace_or_filename : openmc.data.ace.Table or str
ACE table to read from. If given as a string, it is assumed to be
the filename for the ACE file.
name : str
GNDS-conforming name of the material, e.g. c_H_in_H2O. If none is
passed, the appropriate name is guessed based on the name of the ACE
table.
Returns
-------
openmc.data.ThermalScattering
Thermal scattering data
"""
data = ThermalScattering.from_ace(ace_or_filename, name)
# Check if temprature already exists
strT = data.temperatures[0]
if strT in self.temperatures:
warn('S(a,b) data at T={} already exists.'.format(strT))
return
# Check that name matches
if data.name != self.name:
raise ValueError('Data provided for an incorrect material.')
# Add temperature
self.kTs += data.kTs
# Add inelastic cross section and distributions
if data.inelastic is not None:
self.inelastic.xs.update(data.inelastic.xs)
self.inelastic.distribution.update(data.inelastic.distribution)
# Add elastic cross sectoin and angular distribution
if data.elastic is not None:
self.elastic.xs.update(data.elastic.xs)
self.elastic.distribution.update(data.elastic.distribution)
[docs] @classmethod
def from_hdf5(cls, group_or_filename):
"""Generate thermal scattering 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.ThermalScattering
Neutron thermal scattering 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']
if major != HDF5_VERSION_MAJOR:
raise IOError(
'HDF5 data format uses version {}.{} whereas your '
'installation of the OpenMC Python API expects version '
'{}.x.'.format(major, minor, HDF5_VERSION_MAJOR))
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_weight_ratio = group.attrs['atomic_weight_ratio']
energy_max = group.attrs['energy_max']
kTg = group['kTs']
kTs = [dataset[()] for dataset in kTg.values()]
table = cls(name, atomic_weight_ratio, energy_max, kTs)
table.nuclides = [nuc.decode() for nuc in group.attrs['nuclides']]
# Read thermal elastic scattering
if 'elastic' in group[table.temperatures[0]]:
table.elastic = ThermalScatteringReaction.from_hdf5(
group, 'elastic', table.temperatures
)
# Read thermal inelastic scattering
table.inelastic = ThermalScatteringReaction.from_hdf5(
group, 'inelastic', table.temperatures
)
return table
[docs] @classmethod
def from_ace(cls, ace_or_filename, name=None):
"""Generate thermal scattering data from an ACE table
Parameters
----------
ace_or_filename : openmc.data.ace.Table or str
ACE table to read from. If given as a string, it is assumed to be
the filename for the ACE file.
name : str
GNDS-conforming name of the material, e.g. c_H_in_H2O. If none is
passed, the appropriate name is guessed based on the name of the ACE
table.
Returns
-------
openmc.data.ThermalScattering
Thermal scattering data
"""
if isinstance(ace_or_filename, Table):
ace = ace_or_filename
else:
ace = get_table(ace_or_filename)
# Get new name that is GND-consistent
ace_name, xs = ace.name.split('.')
if not xs.endswith('t'):
raise TypeError("{} is not a thermal scattering ACE table.".format(ace))
if name is None:
name = get_thermal_name(ace_name)
# Assign temperature to the running list
kTs = [ace.temperature*EV_PER_MEV]
# Incoherent inelastic scattering cross section
idx = ace.jxs[1]
n_energy = int(ace.xss[idx])
energy = ace.xss[idx+1 : idx+1+n_energy]*EV_PER_MEV
xs = ace.xss[idx+1+n_energy : idx+1+2*n_energy]
inelastic_xs = Tabulated1D(energy, xs)
energy_max = energy[-1]
# Incoherent inelastic angle-energy distribution
continuous = (ace.nxs[7] == 2)
n_energy_out = ace.nxs[4]
if not continuous:
n_mu = ace.nxs[3]
idx = ace.jxs[3]
energy_out = ace.xss[idx:idx + n_energy * n_energy_out *
(n_mu + 2): n_mu + 2]*EV_PER_MEV
energy_out.shape = (n_energy, n_energy_out)
mu_out = ace.xss[idx:idx + n_energy * n_energy_out * (n_mu + 2)]
mu_out.shape = (n_energy, n_energy_out, n_mu+2)
mu_out = mu_out[:, :, 1:]
skewed = (ace.nxs[7] == 1)
distribution = IncoherentInelasticAEDiscrete(energy_out, mu_out, skewed)
else:
n_mu = ace.nxs[3] - 1
idx = ace.jxs[3]
locc = ace.xss[idx:idx + n_energy].astype(int)
n_energy_out = \
ace.xss[idx + n_energy:idx + 2 * n_energy].astype(int)
energy_out = []
mu_out = []
for i in range(n_energy):
idx = locc[i]
# Outgoing energy distribution for incoming energy i
e = ace.xss[idx + 1:idx + 1 + n_energy_out[i]*(n_mu + 3):
n_mu + 3]*EV_PER_MEV
p = ace.xss[idx + 2:idx + 2 + n_energy_out[i]*(n_mu + 3):
n_mu + 3]/EV_PER_MEV
c = ace.xss[idx + 3:idx + 3 + n_energy_out[i]*(n_mu + 3):
n_mu + 3]
eout_i = Tabular(e, p, 'linear-linear', ignore_negative=True)
eout_i.c = c
# Outgoing angle distribution for each
# (incoming, outgoing) energy pair
mu_i = []
for j in range(n_energy_out[i]):
mu = ace.xss[idx + 4:idx + 4 + n_mu]
# The equiprobable angles produced by NJOY are not always
# sorted. This is problematic when the smearing algorithm
# is applied when sampling the angles. We sort the angles
# here, because they are equiprobable, so the order
# doesn't matter.
mu.sort()
# Older versions of NJOY had a bug, and the discrete
# scattering angles could sometimes be less than -1 or
# greater than 1. We check for this here, and warn users.
if mu[0] < -1. or mu[-1] > 1.:
warn('S(a,b) scattering angle for incident energy index '
f'{i} and exit energy index {j} outside of the '
'interval [-1, 1].')
p_mu = 1. / n_mu * np.ones(n_mu)
mu_ij = Discrete(mu, p_mu)
mu_ij.c = np.cumsum(p_mu)
mu_i.append(mu_ij)
idx += 3 + n_mu
# Check if the CDF for the outgoing energy distribution starts
# at 0. For NJOY and FRENDY evaluations, this is never the case,
# and can very rarely lead to negative energies when sampling
# the outgoing energy. From Eq. 7.6 of the ENDF manual, we can
# add an outgoing energy 0 eV that has a PDF of 0 (and of
# course, a CDF of 0 as well).
if eout_i.c[0] > 0.:
eout_i.x = np.insert(eout_i.x, 0, 0.)
eout_i.p = np.insert(eout_i.p, 0, 0.)
eout_i.c = np.insert(eout_i.c, 0, 0.)
# For this added outgoing energy (of 0 eV) we add a set of
# isotropic discrete angles.
dmu = 2. / n_mu
mu = np.linspace(-1. + 0.5*dmu, 1. - 0.5*dmu, n_mu)
p_mu = 1. / n_mu * np.ones(n_mu)
mu_0 = Discrete(mu, p_mu)
mu_0.c = np.cumsum(p_mu)
mu_i.insert(0, mu_0)
# We don't worry about renormalizing the outgoing energy PDF/CDF
# after this manipulation, because it never seems to be
# normalized to begin with (at least with NJOY).
energy_out.append(eout_i)
mu_out.append(mu_i)
# Create correlated angle-energy distribution
breakpoints = [n_energy]
interpolation = [2]
energy = inelastic_xs.x
distribution = IncoherentInelasticAE(
breakpoints, interpolation, energy, energy_out, mu_out)
table = cls(name, ace.atomic_weight_ratio, energy_max, kTs)
T = table.temperatures[0]
table.inelastic = ThermalScatteringReaction(
{T: inelastic_xs}, {T: distribution}
)
# Incoherent/coherent elastic scattering cross section
idx = ace.jxs[4]
if idx != 0:
if ace.nxs[5] in (4, 5):
# Coherent elastic
n_energy = int(ace.xss[idx])
energy = ace.xss[idx + 1: idx + 1 + n_energy]*EV_PER_MEV
P = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy]
coherent_xs = CoherentElastic(energy, P*EV_PER_MEV)
coherent_dist = CoherentElasticAE(coherent_xs)
# Coherent elastic shouldn't have angular distributions listed
n_mu = ace.nxs[6] + 1
assert n_mu == 0
if ace.nxs[5] in (3, 5):
# Incoherent elastic scattering -- first determine if both
# incoherent and coherent are present (mixed)
mixed = (ace.nxs[5] == 5)
# Get cross section values
idx = ace.jxs[7] if mixed else ace.jxs[4]
n_energy = int(ace.xss[idx])
energy = ace.xss[idx + 1: idx + 1 + n_energy]*EV_PER_MEV
values = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy]
incoherent_xs = Tabulated1D(energy, values)
# Angular distribution
n_mu = (ace.nxs[8] if mixed else ace.nxs[6]) + 1
assert n_mu > 0
idx = ace.jxs[9] if mixed else ace.jxs[6]
mu_out = ace.xss[idx:idx + n_energy * n_mu]
mu_out.shape = (n_energy, n_mu)
incoherent_dist = IncoherentElasticAEDiscrete(mu_out)
if ace.nxs[5] == 3:
xs = incoherent_xs
distribution = incoherent_dist
elif ace.nxs[5] == 4:
xs = coherent_xs
distribution = coherent_dist
else:
# Create mixed cross section -- note that coherent must come
# first due to assumption on C++ side
xs = Sum([coherent_xs, incoherent_xs])
# Create mixed distribution
distribution = MixedElasticAE(coherent_dist, incoherent_dist)
table.elastic = ThermalScatteringReaction({T: xs}, {T: distribution})
# Get relevant nuclides -- NJOY only allows one to specify three
# nuclides that the S(a,b) table applies to. Thus, for all elements
# other than H and Fe, we automatically add all the naturally-occurring
# isotopes.
for zaid, awr in ace.pairs:
if zaid > 0:
Z, A = divmod(zaid, 1000)
element = ATOMIC_SYMBOL[Z]
if element in ['H', 'Fe']:
table.nuclides.append(element + str(A))
else:
if element + '0' not in table.nuclides:
table.nuclides.append(element + '0')
for isotope, _ in isotopes(element):
if isotope not in table.nuclides:
table.nuclides.append(isotope)
return table
[docs] @classmethod
def from_njoy(cls, filename, filename_thermal, temperatures=None,
evaluation=None, evaluation_thermal=None,
use_endf_data=True, **kwargs):
"""Generate thermal scattering data by running NJOY.
Parameters
----------
filename : str
Path to ENDF neutron sublibrary file
filename_thermal : str
Path to ENDF thermal scattering sublibrary file
temperatures : iterable of float
Temperatures in Kelvin to produce data at. If omitted, data is
produced at all temperatures in the ENDF thermal scattering
sublibrary.
evaluation : openmc.data.endf.Evaluation, optional
If the ENDF neutron sublibrary file contains multiple material
evaluations, this argument indicates which evaluation to use.
evaluation_thermal : openmc.data.endf.Evaluation, optional
If the ENDF thermal scattering sublibrary file contains multiple
material evaluations, this argument indicates which evaluation to
use.
use_endf_data : bool
If the material has incoherent elastic scattering, the ENDF data
will be used rather than the ACE data.
**kwargs
Keyword arguments passed to :func:`openmc.data.njoy.make_ace_thermal`
Returns
-------
data : openmc.data.ThermalScattering
Thermal scattering data
"""
with tempfile.TemporaryDirectory() as tmpdir:
# Run NJOY to create an ACE library
kwargs.setdefault('output_dir', tmpdir)
kwargs.setdefault('ace', os.path.join(kwargs['output_dir'], 'ace'))
kwargs['evaluation'] = evaluation
kwargs['evaluation_thermal'] = evaluation_thermal
make_ace_thermal(filename, filename_thermal, temperatures, **kwargs)
# Create instance from ACE tables within library
lib = Library(kwargs['ace'])
name = kwargs.get('table_name')
data = cls.from_ace(lib.tables[0], name=name)
for table in lib.tables[1:]:
data.add_temperature_from_ace(table, name=name)
# Load ENDF data to replace incoherent elastic
if use_endf_data:
data_endf = cls.from_endf(filename_thermal)
if data_endf.elastic is not None:
# Get appropriate temperatures
if temperatures is None:
temperatures = data_endf.temperatures
else:
temperatures = [_temperature_str(t) for t in temperatures]
# Replace ACE data with ENDF data
rx, rx_endf = data.elastic, data_endf.elastic
for t in temperatures:
if isinstance(rx_endf.xs[t], (IncoherentElastic, Sum)):
rx.xs[t] = rx_endf.xs[t]
rx.distribution[t] = rx_endf.distribution[t]
return data
[docs] @classmethod
def from_endf(cls, ev_or_filename):
"""Generate thermal scattering data from an ENDF file
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.ThermalScattering
Thermal scattering data
"""
if isinstance(ev_or_filename, endf.Evaluation):
ev = ev_or_filename
else:
ev = endf.Evaluation(ev_or_filename)
# Read coherent/incoherent elastic data
elastic = None
if (7, 2) in ev.section:
# Define helper functions to avoid duplication
def get_coherent_elastic(file_obj):
# Get structure factor at first temperature
params, S = endf.get_tab1_record(file_obj)
strT = _temperature_str(params[0])
n_temps = params[2]
bragg_edges = S.x
xs = {strT: CoherentElastic(bragg_edges, S.y)}
distribution = {strT: CoherentElasticAE(xs[strT])}
# Get structure factor for subsequent temperatures
for _ in range(n_temps):
params, S = endf.get_list_record(file_obj)
strT = _temperature_str(params[0])
xs[strT] = CoherentElastic(bragg_edges, S)
distribution[strT] = CoherentElasticAE(xs[strT])
return xs, distribution
def get_incoherent_elastic(file_obj):
params, W = endf.get_tab1_record(file_obj)
bound_xs = params[0]
xs = {}
distribution = {}
for T, debye_waller in zip(W.x, W.y):
strT = _temperature_str(T)
xs[strT] = IncoherentElastic(bound_xs, debye_waller)
distribution[strT] = IncoherentElasticAE(debye_waller)
return xs, distribution
file_obj = StringIO(ev.section[7, 2])
lhtr = endf.get_head_record(file_obj)[2]
if lhtr == 1:
# coherent elastic
xs, distribution = get_coherent_elastic(file_obj)
elif lhtr == 2:
# incoherent elastic
xs, distribution = get_incoherent_elastic(file_obj)
elif lhtr == 3:
# mixed coherent / incoherent elastic
xs_c, dist_c = get_coherent_elastic(file_obj)
xs_i, dist_i = get_incoherent_elastic(file_obj)
assert sorted(xs_c) == sorted(xs_i)
xs = {T: Sum([xs_c[T], xs_i[T]]) for T in xs_c}
distribution = {T: MixedElasticAE(dist_c[T], dist_i[T]) for T in dist_c}
elastic = ThermalScatteringReaction(xs, distribution)
# Read incoherent inelastic data
assert (7, 4) in ev.section, 'No MF=7, MT=4 found in thermal scattering'
file_obj = StringIO(ev.section[7, 4])
params = endf.get_head_record(file_obj)
data = {'symmetric': params[4] == 0}
# Get information about principal atom
params, B = endf.get_list_record(file_obj)
data['log'] = bool(params[2])
data['free_atom_xs'] = B[0]
data['epsilon'] = B[1]
data['A0'] = awr = B[2]
data['e_max'] = energy_max = B[3]
data['M0'] = B[5]
# Get information about non-principal atoms
n_non_principal = params[5]
data['non_principal'] = []
NonPrincipal = namedtuple('NonPrincipal', ['func', 'xs', 'A', 'M'])
for i in range(1, n_non_principal + 1):
func = {0.0: 'SCT', 1.0: 'free gas', 2.0: 'diffusive'}[B[6*i]]
xs = B[6*i + 1]
A = B[6*i + 2]
M = B[6*i + 5]
data['non_principal'].append(NonPrincipal(func, xs, A, M))
# Get S(alpha,beta,T)
kTs = []
if data['free_atom_xs'] > 0.0:
params, _ = endf.get_tab2_record(file_obj)
n_beta = params[5]
sab = {'beta': np.empty(n_beta)}
for i in range(n_beta):
params, S = endf.get_tab1_record(file_obj)
T0, beta, lt = params[:3]
if i == 0:
sab['alpha'] = alpha = S.x
sab[T0] = np.empty((alpha.size, n_beta))
kTs.append(K_BOLTZMANN * T0)
sab['beta'][i] = beta
sab[T0][:, i] = S.y
for _ in range(lt):
params, S = endf.get_list_record(file_obj)
T = params[0]
if i == 0:
sab[T] = np.empty((alpha.size, n_beta))
kTs.append(K_BOLTZMANN * T)
sab[T][:, i] = S
data['sab'] = sab
# Get effective temperature for each atom
_, Teff = endf.get_tab1_record(file_obj)
data['effective_temperature'] = [Teff]
for atom in data['non_principal']:
if atom.func == 'SCT':
_, Teff = endf.get_tab1_record(file_obj)
data['effective_temperature'].append(Teff)
name = ev.target['zsymam'].strip()
instance = cls(name, awr, energy_max, kTs)
if elastic is not None:
instance.elastic = elastic
# Currently we don't have a proper cross section or distribution for
# incoherent inelastic, so we just create an empty object and attach
# all the data as a dictionary
instance.inelastic = ThermalScatteringReaction(None, None)
instance.inelastic.data = data
return instance