mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-21 14:35:27 -04:00
875 lines
32 KiB
Python
875 lines
32 KiB
Python
from collections.abc import Mapping, MutableMapping
|
|
from io import StringIO
|
|
from math import log10
|
|
from numbers import Integral, Real
|
|
import os
|
|
import tempfile
|
|
from warnings import warn
|
|
|
|
import numpy as np
|
|
import h5py
|
|
|
|
from . import HDF5_VERSION, HDF5_VERSION_MAJOR
|
|
from .ace import Library, Table, get_table, get_metadata
|
|
from .data import ATOMIC_SYMBOL, K_BOLTZMANN, EV_PER_MEV, gnds_name
|
|
from .endf import (
|
|
Evaluation, SUM_RULES, as_evaluation, get_head_record, get_tab1_record,
|
|
get_evaluations)
|
|
from .fission_energy import FissionEnergyRelease
|
|
from .function import Tabulated1D, Sum, ResonancesWithBackground
|
|
from .njoy import make_ace, make_pendf
|
|
from .product import Product
|
|
from .reaction import Reaction, _get_photon_products_ace, FISSION_MTS
|
|
from . import resonance as res
|
|
from . import resonance_covariance as res_cov
|
|
from .urr import ProbabilityTables
|
|
import openmc.checkvalue as cv
|
|
from openmc.mixin import EqualityMixin
|
|
|
|
|
|
# Fractions of resonance widths used for reconstructing resonances
|
|
_RESONANCE_ENERGY_GRID = np.logspace(-3, 3, 61)
|
|
|
|
|
|
class IncidentNeutron(EqualityMixin):
|
|
"""Continuous-energy neutron interaction data.
|
|
|
|
This class stores data derived from an ENDF-6 format neutron interaction
|
|
sublibrary. Instances of this class are not normally instantiated by the
|
|
user but rather created using the factory methods
|
|
:meth:`IncidentNeutron.from_hdf5`, :meth:`IncidentNeutron.from_ace`, and
|
|
:meth:`IncidentNeutron.from_endf`.
|
|
|
|
Parameters
|
|
----------
|
|
name : str
|
|
Name of the nuclide using the GNDS naming convention
|
|
atomic_number : int
|
|
Number of protons 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 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_number : int
|
|
Number of protons 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 GNDS naming convention
|
|
reactions : dict
|
|
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.
|
|
resonances : openmc.data.Resonances or None
|
|
Resonance parameters
|
|
resonance_covariance : openmc.data.ResonanceCovariance or None
|
|
Covariance for resonance parameters
|
|
temperatures : list of str
|
|
List of string representations of 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.
|
|
urr : dict
|
|
Dictionary whose keys are temperatures (e.g., '294K') and values are
|
|
unresolved resonance region probability tables.
|
|
|
|
"""
|
|
|
|
def __init__(self, name, atomic_number, mass_number, metastable,
|
|
atomic_weight_ratio, kTs):
|
|
self.name = name
|
|
self.atomic_number = atomic_number
|
|
self.mass_number = mass_number
|
|
self.metastable = metastable
|
|
self.atomic_weight_ratio = atomic_weight_ratio
|
|
self.kTs = kTs
|
|
self.energy = {}
|
|
self._fission_energy = None
|
|
self.reactions = {}
|
|
self._urr = {}
|
|
self._resonances = None
|
|
|
|
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 f"<IncidentNeutron: {self.name}>"
|
|
|
|
def __iter__(self):
|
|
return iter(self.reactions.values())
|
|
|
|
@property
|
|
def name(self):
|
|
return self._name
|
|
|
|
@name.setter
|
|
def name(self, name):
|
|
cv.check_type('name', name, str)
|
|
self._name = name
|
|
|
|
@property
|
|
def atomic_number(self):
|
|
return self._atomic_number
|
|
|
|
@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
|
|
|
|
@property
|
|
def mass_number(self):
|
|
return self._mass_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
|
|
|
|
@property
|
|
def metastable(self):
|
|
return self._metastable
|
|
|
|
@metastable.setter
|
|
def metastable(self, metastable):
|
|
cv.check_type('metastable', metastable, Integral)
|
|
cv.check_greater_than('metastable', metastable, 0, True)
|
|
self._metastable = metastable
|
|
|
|
@property
|
|
def atomic_weight_ratio(self):
|
|
return self._atomic_weight_ratio
|
|
|
|
@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
|
|
|
|
@property
|
|
def reactions(self):
|
|
return self._reactions
|
|
|
|
@reactions.setter
|
|
def reactions(self, reactions):
|
|
cv.check_type('reactions', reactions, Mapping)
|
|
self._reactions = reactions
|
|
|
|
@property
|
|
def resonances(self):
|
|
return self._resonances
|
|
|
|
@resonances.setter
|
|
def resonances(self, resonances):
|
|
cv.check_type('resonances', resonances, res.Resonances)
|
|
self._resonances = resonances
|
|
|
|
@property
|
|
def resonance_covariance(self):
|
|
return self._resonance_covariance
|
|
|
|
@resonance_covariance.setter
|
|
def resonance_covariance(self, resonance_covariance):
|
|
cv.check_type('resonance covariance', resonance_covariance,
|
|
res_cov.ResonanceCovariances)
|
|
self._resonance_covariance = resonance_covariance
|
|
|
|
@property
|
|
def urr(self):
|
|
return self._urr
|
|
|
|
@urr.setter
|
|
def urr(self, urr):
|
|
cv.check_type('probability table dictionary', urr, MutableMapping)
|
|
for key, value in urr:
|
|
cv.check_type('probability table temperature', key, str)
|
|
cv.check_type('probability tables', value, ProbabilityTables)
|
|
self._urr = urr
|
|
|
|
@property
|
|
def temperatures(self):
|
|
return [f"{int(round(kT / K_BOLTZMANN))}K" for kT in self.kTs]
|
|
|
|
@property
|
|
def atomic_symbol(self):
|
|
return ATOMIC_SYMBOL[self.atomic_number]
|
|
|
|
def add_temperature_from_ace(self, ace_or_filename, metastable_scheme='nndc'):
|
|
"""Append data 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.
|
|
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.
|
|
|
|
"""
|
|
|
|
data = IncidentNeutron.from_ace(ace_or_filename, metastable_scheme)
|
|
|
|
# Check if temprature already exists
|
|
strT = data.temperatures[0]
|
|
if strT in self.temperatures:
|
|
warn(f'Cross sections at T={strT} already exist.')
|
|
return
|
|
|
|
# Check that name matches
|
|
if data.name != self.name:
|
|
raise ValueError('Data provided for an incorrect nuclide.')
|
|
|
|
# Add temperature
|
|
self.kTs += data.kTs
|
|
|
|
# Add energy grid
|
|
self.energy[strT] = data.energy[strT]
|
|
|
|
# Add normal and redundant reactions
|
|
for mt in data.reactions:
|
|
if mt in self:
|
|
self[mt].xs[strT] = data[mt].xs[strT]
|
|
else:
|
|
warn("Tried to add cross sections for MT={} at T={} but this "
|
|
"reaction doesn't exist.".format(mt, strT))
|
|
|
|
# Add probability tables
|
|
if strT in data.urr:
|
|
self.urr[strT] = data.urr[strT]
|
|
|
|
def add_elastic_0K_from_endf(self, filename, overwrite=False, **kwargs):
|
|
"""Append 0K elastic scattering cross section from an ENDF file.
|
|
|
|
Parameters
|
|
----------
|
|
filename : str
|
|
Path to ENDF file
|
|
overwrite : bool
|
|
If existing 0 K data is present, this flag can be used to indicate
|
|
that it should be overwritten. Otherwise, an exception will be
|
|
thrown.
|
|
**kwargs
|
|
Keyword arguments passed to :func:`openmc.data.njoy.make_pendf`
|
|
|
|
Raises
|
|
------
|
|
ValueError
|
|
If 0 K data is already present and the `overwrite` parameter is
|
|
False.
|
|
|
|
"""
|
|
# Check for existing data
|
|
if '0K' in self.energy and not overwrite:
|
|
raise ValueError('0 K data already exists for this nuclide.')
|
|
|
|
with tempfile.TemporaryDirectory() as tmpdir:
|
|
# Set arguments for make_pendf
|
|
pendf_path = os.path.join(tmpdir, 'pendf')
|
|
kwargs.setdefault('output_dir', tmpdir)
|
|
kwargs.setdefault('pendf', pendf_path)
|
|
|
|
# Run NJOY to create a pointwise ENDF file
|
|
make_pendf(filename, **kwargs)
|
|
|
|
# Add 0K elastic scattering cross section
|
|
pendf = Evaluation(pendf_path)
|
|
file_obj = StringIO(pendf.section[3, 2])
|
|
get_head_record(file_obj)
|
|
params, xs = get_tab1_record(file_obj)
|
|
self.energy['0K'] = xs.x
|
|
self[2].xs['0K'] = xs
|
|
|
|
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 export_to_hdf5(self, path, mode='a', libver='earliest'):
|
|
"""Export incident neutron data 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.
|
|
|
|
"""
|
|
# If data come from ENDF, don't allow exporting to HDF5
|
|
if hasattr(self, '_evaluation'):
|
|
raise NotImplementedError('Cannot export incident neutron data that '
|
|
'originated from an ENDF file.')
|
|
|
|
# Open file and write version
|
|
with h5py.File(str(path), mode, libver=libver) as f:
|
|
f.attrs['filetype'] = np.bytes_('data_neutron')
|
|
f.attrs['version'] = np.array(HDF5_VERSION)
|
|
|
|
# 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
|
|
ktg = g.create_group('kTs')
|
|
for i, temperature in enumerate(self.temperatures):
|
|
ktg.create_dataset(temperature, data=self.kTs[i])
|
|
|
|
# Write energy grid
|
|
eg = g.create_group('energy')
|
|
for temperature in self.temperatures:
|
|
eg.create_dataset(temperature, data=self.energy[temperature])
|
|
|
|
# Write 0K energy grid if needed
|
|
if '0K' in self.energy and '0K' not in eg:
|
|
eg.create_dataset('0K', data=self.energy['0K'])
|
|
|
|
# Write reaction data
|
|
rxs_group = g.create_group('reactions')
|
|
for rx in self.reactions.values():
|
|
# Skip writing redundant reaction if it doesn't have photon
|
|
# production or is a summed transmutation reaction. MT=4 is also
|
|
# sometimes needed for probability tables. Also write gas
|
|
# production, heating, and damage energy production.
|
|
if rx.redundant:
|
|
photon_rx = any(p.particle == 'photon' for p in rx.products)
|
|
keep_mts = (4, 16, 103, 104, 105, 106, 107,
|
|
203, 204, 205, 206, 207, 301, 444, 901)
|
|
if not (photon_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 total nu data if available
|
|
if len(rx.derived_products) > 0 and 'total_nu' not in g:
|
|
tgroup = g.create_group('total_nu')
|
|
rx.derived_products[0].to_hdf5(tgroup)
|
|
|
|
# Write unresolved resonance probability tables
|
|
if self.urr:
|
|
urr_group = g.create_group('urr')
|
|
for temperature, urr in self.urr.items():
|
|
tgroup = urr_group.create_group(temperature)
|
|
urr.to_hdf5(tgroup)
|
|
|
|
# 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)
|
|
|
|
@classmethod
|
|
def from_hdf5(cls, group_or_filename):
|
|
"""Generate continuous-energy neutron 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.IncidentNeutron
|
|
Continuous-energy neutron 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']
|
|
kTg = group['kTs']
|
|
kTs = []
|
|
for temp in kTg:
|
|
kTs.append(kTg[temp][()])
|
|
|
|
data = cls(name, atomic_number, mass_number, metastable,
|
|
atomic_weight_ratio, kTs)
|
|
|
|
# Read energy grid
|
|
e_group = group['energy']
|
|
for temperature, dset in e_group.items():
|
|
data.energy[temperature] = dset[()]
|
|
|
|
# Read reaction data
|
|
rxs_group = group['reactions']
|
|
for name, obj in sorted(rxs_group.items()):
|
|
if name.startswith('reaction_'):
|
|
rx = Reaction.from_hdf5(obj, data.energy)
|
|
data.reactions[rx.mt] = rx
|
|
|
|
# Read total nu data if available
|
|
if rx.mt in FISSION_MTS and 'total_nu' in group:
|
|
tgroup = group['total_nu']
|
|
rx.derived_products.append(Product.from_hdf5(tgroup))
|
|
|
|
# Read unresolved resonance probability tables
|
|
if 'urr' in group:
|
|
urr_group = group['urr']
|
|
for temperature, tgroup in urr_group.items():
|
|
data.urr[temperature] = ProbabilityTables.from_hdf5(tgroup)
|
|
|
|
# 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)
|
|
|
|
return data
|
|
|
|
@classmethod
|
|
def from_ace(cls, ace_or_filename, metastable_scheme='nndc'):
|
|
"""Generate incident neutron 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.IncidentNeutron
|
|
Incident neutron 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('c'):
|
|
raise TypeError(
|
|
f"{ace} is not a continuous-energy neutron ACE table.")
|
|
name, element, Z, mass_number, metastable = \
|
|
get_metadata(int(zaid), metastable_scheme)
|
|
|
|
# Assign temperature to the running list
|
|
kTs = [ace.temperature*EV_PER_MEV]
|
|
|
|
data = cls(name, Z, mass_number, metastable,
|
|
ace.atomic_weight_ratio, kTs)
|
|
|
|
# Get string of temperature to use as a dictionary key
|
|
strT = data.temperatures[0]
|
|
|
|
# Read energy grid
|
|
n_energy = ace.nxs[3]
|
|
i = ace.jxs[1]
|
|
energy = ace.xss[i : i + n_energy]*EV_PER_MEV
|
|
data.energy[strT] = energy
|
|
total_xs = ace.xss[i + n_energy : i + 2*n_energy]
|
|
absorption_xs = ace.xss[i + 2*n_energy : i + 3*n_energy]
|
|
heating_number = ace.xss[i + 4*n_energy : i + 5*n_energy]*EV_PER_MEV
|
|
|
|
# Create redundant reaction for total (MT=1)
|
|
total = Reaction(1)
|
|
total.xs[strT] = Tabulated1D(energy, total_xs)
|
|
total.redundant = True
|
|
data.reactions[1] = total
|
|
|
|
# Create redundant reaction for absorption (MT=101)
|
|
if np.count_nonzero(absorption_xs) > 0:
|
|
absorption = Reaction(101)
|
|
absorption.xs[strT] = Tabulated1D(energy, absorption_xs)
|
|
absorption.redundant = True
|
|
data.reactions[101] = absorption
|
|
|
|
# Create redundant reaction for heating (MT=301)
|
|
heating = Reaction(301)
|
|
heating.xs[strT] = 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):
|
|
rx = Reaction.from_ace(ace, i)
|
|
data.reactions[rx.mt] = rx
|
|
|
|
# Some photon production reactions may be assigned to MTs that don't
|
|
# exist, usually MT=4. In this case, we create a new reaction and add
|
|
# them
|
|
n_photon_reactions = ace.nxs[6]
|
|
photon_mts = ace.xss[ace.jxs[13]:ace.jxs[13] +
|
|
n_photon_reactions].astype(int)
|
|
|
|
for mt in np.unique(photon_mts // 1000):
|
|
if mt not in data:
|
|
if mt not in SUM_RULES:
|
|
warn('Photon production is present for MT={} but no '
|
|
'cross section is given.'.format(mt))
|
|
continue
|
|
|
|
# Create redundant reaction with appropriate cross section
|
|
mts = data.get_reaction_components(mt)
|
|
if len(mts) == 0:
|
|
warn('Photon production is present for MT={} but no '
|
|
'reaction components exist.'.format(mt))
|
|
continue
|
|
|
|
# Determine redundant cross section
|
|
rx = data._get_redundant_reaction(mt, mts)
|
|
rx.products += _get_photon_products_ace(ace, rx)
|
|
data.reactions[mt] = rx
|
|
|
|
# For transmutation reactions, sometimes only individual levels are
|
|
# present in an ACE file, e.g. MT=600-649 instead of the summation
|
|
# MT=103. In this case, if a user wants to tally (n,p), OpenMC doesn't
|
|
# know about the total cross section. Here, we explicitly create a
|
|
# redundant reaction for this purpose.
|
|
for mt in (16, 103, 104, 105, 106, 107):
|
|
if mt not in data:
|
|
# Determine if any individual levels are present
|
|
mts = data.get_reaction_components(mt)
|
|
if len(mts) == 0:
|
|
continue
|
|
|
|
# Determine redundant cross section
|
|
rx = data._get_redundant_reaction(mt, mts)
|
|
data.reactions[mt] = rx
|
|
|
|
# Make sure redundant cross sections that are present in an ACE file get
|
|
# marked as such
|
|
for rx in data:
|
|
mts = data.get_reaction_components(rx.mt)
|
|
if mts != [rx.mt]:
|
|
rx.redundant = True
|
|
if rx.mt in (203, 204, 205, 206, 207, 444):
|
|
rx.redundant = True
|
|
|
|
# Read unresolved resonance probability tables
|
|
urr = ProbabilityTables.from_ace(ace)
|
|
if urr is not None:
|
|
data.urr[strT] = urr
|
|
|
|
return data
|
|
|
|
@classmethod
|
|
def from_endf(cls, ev_or_filename, covariance=False):
|
|
"""Generate incident neutron continuous-energy data from an ENDF evaluation
|
|
|
|
Parameters
|
|
----------
|
|
ev_or_filename : openmc.data.endf.Evaluation, endf.Material, or str
|
|
ENDF evaluation to read from. If given as a string, it is assumed to
|
|
be the filename for the ENDF file.
|
|
|
|
covariance : bool
|
|
Flag to indicate whether or not covariance data from File 32 should be
|
|
retrieved
|
|
|
|
Returns
|
|
-------
|
|
openmc.data.IncidentNeutron
|
|
Incident neutron continuous-energy data
|
|
|
|
"""
|
|
ev = as_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']
|
|
temperature = ev.target['temperature']
|
|
|
|
# Determine name
|
|
name = gnds_name(atomic_number, mass_number, metastable)
|
|
|
|
# Instantiate incident neutron data
|
|
data = cls(name, atomic_number, mass_number, metastable,
|
|
atomic_weight_ratio, [temperature])
|
|
|
|
if (2, 151) in ev.section:
|
|
data.resonances = res.Resonances.from_endf(ev)
|
|
|
|
if (32, 151) in ev.section and covariance:
|
|
data.resonance_covariance = (
|
|
res_cov.ResonanceCovariances.from_endf(ev, data.resonances)
|
|
)
|
|
|
|
# Read each reaction
|
|
for mf, mt, nc, mod in ev.reaction_list:
|
|
if mf == 3:
|
|
data.reactions[mt] = Reaction.from_endf(ev, mt)
|
|
|
|
# Replace cross sections for elastic, capture, fission
|
|
try:
|
|
if any(isinstance(r, res._RESOLVED) for r in data.resonances):
|
|
for mt in (2, 102, 18):
|
|
if mt in data.reactions:
|
|
rx = data.reactions[mt]
|
|
rx.xs['0K'] = ResonancesWithBackground(
|
|
data.resonances, rx.xs['0K'], mt)
|
|
except ValueError:
|
|
# Thrown if multiple resolved ranges (e.g. Pu239 in ENDF/B-VII.1)
|
|
pass
|
|
|
|
# If first-chance, second-chance, etc. fission are present, check
|
|
# whether energy distributions were specified in MF=5. If not, copy the
|
|
# energy distribution from MT=18.
|
|
for mt, rx in data.reactions.items():
|
|
if mt in (19, 20, 21, 38):
|
|
if (5, mt) not in ev.section:
|
|
if rx.products:
|
|
neutron = data.reactions[18].products[0]
|
|
rx.products[0].applicability = neutron.applicability
|
|
rx.products[0].distribution = neutron.distribution
|
|
|
|
# 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_njoy(cls, filename, temperatures=None, evaluation=None, **kwargs):
|
|
"""Generate incident neutron data by running NJOY.
|
|
|
|
Parameters
|
|
----------
|
|
filename : str
|
|
Path to ENDF file
|
|
temperatures : iterable of float
|
|
Temperatures in Kelvin to produce data at. If omitted, data is
|
|
produced at room temperature (293.6 K)
|
|
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`
|
|
|
|
Returns
|
|
-------
|
|
data : openmc.data.IncidentNeutron
|
|
Incident neutron 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", "heatr", "broadr", "gaspr", "purr"):
|
|
kwargs.setdefault(key, os.path.join(kwargs["output_dir"], key))
|
|
kwargs['evaluation'] = evaluation
|
|
make_ace(filename, temperatures, **kwargs)
|
|
|
|
# Create instance from ACE tables within library
|
|
lib = Library(kwargs['acer'])
|
|
data = cls.from_ace(lib.tables[0])
|
|
for table in lib.tables[1:]:
|
|
data.add_temperature_from_ace(table)
|
|
|
|
# Use name based on ENDF evaluation. The name assigned by from_ace
|
|
# may be wrong for higher metastable states (e.g., Hf178_m2)
|
|
ev = evaluation if evaluation is not None else Evaluation(filename)
|
|
data.name = ev.gnds_name
|
|
|
|
# Add 0K elastic scattering cross section
|
|
if '0K' not in data.energy:
|
|
pendf = Evaluation(kwargs['pendf'])
|
|
file_obj = StringIO(pendf.section[3, 2])
|
|
get_head_record(file_obj)
|
|
params, xs = get_tab1_record(file_obj)
|
|
data.energy['0K'] = xs.x
|
|
data[2].xs['0K'] = xs
|
|
|
|
# Add fission energy release data
|
|
if (1, 458) in ev.section:
|
|
data.fission_energy = f = FissionEnergyRelease.from_endf(ev, data)
|
|
else:
|
|
f = None
|
|
|
|
# For energy deposition, we want to store two different KERMAs:
|
|
# one calculated assuming outgoing photons deposit their energy
|
|
# locally, and one calculated assuming they carry their energy
|
|
# away. This requires two HEATR runs (which make_ace does by
|
|
# default). Here, we just need to correct for the fact that NJOY
|
|
# uses a fission heating number of h = EFR, whereas we want:
|
|
#
|
|
# 1) h = EFR + EGP + EGD + EB (for local case)
|
|
# 2) h = EFR + EB (for non-local case)
|
|
#
|
|
# The best way to handle this is to subtract off the fission
|
|
# KERMA that NJOY calculates and add back exactly what we want.
|
|
|
|
# If NJOY is not run with HEATR at all, skip everything below
|
|
if not kwargs["heatr"]:
|
|
return data
|
|
|
|
# Helper function to get a cross section from an ENDF file on a
|
|
# given energy grid
|
|
def get_file3_xs(ev, mt, E):
|
|
xs = ev.section_data[3, mt]['sigma']
|
|
return xs(E)
|
|
|
|
heating_local = Reaction(901)
|
|
heating_local.redundant = True
|
|
|
|
heatr_evals = get_evaluations(kwargs["heatr"])
|
|
heatr_local_evals = get_evaluations(kwargs["heatr"] + "_local")
|
|
|
|
for ev, ev_local, temp in zip(heatr_evals, heatr_local_evals, data.temperatures):
|
|
# Get total KERMA (originally from ACE file) and energy grid
|
|
kerma = data.reactions[301].xs[temp]
|
|
E = kerma.x
|
|
|
|
if f is not None:
|
|
# Replace fission KERMA with (EFR + EB)*sigma_f
|
|
fission = data[18].xs[temp]
|
|
kerma_fission = get_file3_xs(ev, 318, E)
|
|
kerma.y = kerma.y - kerma_fission + (
|
|
f.fragments(E) + f.betas(E)) * fission(E)
|
|
|
|
# For local KERMA, we first need to get the values from the
|
|
# HEATR run with photon energy deposited locally and put
|
|
# them on the same energy grid
|
|
kerma_local = get_file3_xs(ev_local, 301, E)
|
|
|
|
if f is not None:
|
|
# When photons deposit their energy locally, we replace the
|
|
# fission KERMA with (EFR + EGP + EGD + EB)*sigma_f
|
|
kerma_fission_local = get_file3_xs(ev_local, 318, E)
|
|
kerma_local = kerma_local - kerma_fission_local + (
|
|
f.fragments(E) + f.prompt_photons(E)
|
|
+ f.delayed_photons(E) + f.betas(E))*fission(E)
|
|
|
|
heating_local.xs[temp] = Tabulated1D(E, kerma_local)
|
|
|
|
data.reactions[901] = heating_local
|
|
|
|
return data
|
|
|
|
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 = Reaction(mt)
|
|
# Get energy grid
|
|
for strT in self.temperatures:
|
|
energy = self.energy[strT]
|
|
xss = [self.reactions[mt_i].xs[strT] for mt_i in mts]
|
|
idx = min([xs._threshold_idx if hasattr(xs, '_threshold_idx')
|
|
else 0 for xs in xss])
|
|
rx.xs[strT] = Tabulated1D(energy[idx:], Sum(xss)(energy[idx:]))
|
|
rx.xs[strT]._threshold_idx = idx
|
|
|
|
rx.redundant = True
|
|
|
|
return rx
|