mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 14:15:42 -04:00
Merge pull request #1546 from paulromano/pullrequestinc-part4
Address review from PullRequest Inc. (Part 4)
This commit is contained in:
Amanda Lund
GitHub
commit
74133e5014
14 changed files with 226 additions and 184 deletions
|
|
@ -1,14 +1,31 @@
|
|||
#include <stdlib.h>
|
||||
|
||||
//! Convert string representation of a floating point number into a double
|
||||
//
|
||||
//! This function handles converting floating point numbers from an ENDF 11
|
||||
//! character field into a double, covering all the corner cases. Floating point
|
||||
//! numbers are allowed to contain whitespace (which is ignored). Also, in
|
||||
//! exponential notation, it allows the 'e' to be omitted. A field containing
|
||||
//! only whitespace is to be interpreted as a zero.
|
||||
//
|
||||
//! \param buffer character input from an ENDF file
|
||||
//! \param n Length of character input
|
||||
//! \return Floating point number
|
||||
|
||||
double cfloat_endf(const char* buffer, int n)
|
||||
{
|
||||
char arr[12]; // 11 characters plus a null terminator
|
||||
int j = 0; // current position in arr
|
||||
int found_significand = 0;
|
||||
int found_exponent = 0;
|
||||
|
||||
// limit n to 11 characters
|
||||
n = n > 11 ? 11 : n;
|
||||
|
||||
for (int i = 0; i < n; ++i) {
|
||||
// Skip whitespace characters
|
||||
char c = buffer[i];
|
||||
|
||||
// Skip whitespace characters
|
||||
if (c == ' ') continue;
|
||||
|
||||
if (found_significand) {
|
||||
|
|
|
|||
|
|
@ -37,12 +37,18 @@ def linearize(x, f, tolerance=0.001):
|
|||
y_stack.insert(0, f(x[i + 1]))
|
||||
|
||||
while True:
|
||||
# Get the bounding points currently on the stack
|
||||
x_high, x_low = x_stack[-2:]
|
||||
y_high, y_low = y_stack[-2:]
|
||||
|
||||
# Evaluate the function at the midpoint
|
||||
x_mid = 0.5*(x_low + x_high)
|
||||
y_mid = f(x_mid)
|
||||
|
||||
# Linearly interpolate between the bounding points
|
||||
y_interp = y_low + (y_high - y_low)/(x_high - x_low)*(x_mid - x_low)
|
||||
|
||||
# Check the error on the interpolated point and compare to tolerance
|
||||
error = abs((y_interp - y_mid)/y_mid)
|
||||
if error > tolerance:
|
||||
x_stack.insert(-1, x_mid)
|
||||
|
|
|
|||
|
|
@ -245,6 +245,7 @@ class KalbachMann(AngleEnergy):
|
|||
eout_i = Mixture([p_discrete, 1. - p_discrete],
|
||||
[eout_discrete, eout_continuous])
|
||||
|
||||
# Precompound factor and slope are on rows 3 and 4, respectively
|
||||
km_r = Tabulated1D(data[0, j:j+n], data[3, j:j+n])
|
||||
km_a = Tabulated1D(data[0, j:j+n], data[4, j:j+n])
|
||||
|
||||
|
|
|
|||
|
|
@ -5,7 +5,8 @@ import h5py
|
|||
import numpy as np
|
||||
|
||||
import openmc.checkvalue as cv
|
||||
from openmc.mixin import EqualityMixin
|
||||
from ..exceptions import DataError
|
||||
from ..mixin import EqualityMixin
|
||||
from . import WMP_VERSION, WMP_VERSION_MAJOR
|
||||
from .data import K_BOLTZMANN
|
||||
|
||||
|
|
@ -332,19 +333,21 @@ class WindowedMultipole(EqualityMixin):
|
|||
|
||||
if isinstance(group_or_filename, h5py.Group):
|
||||
group = group_or_filename
|
||||
need_to_close = False
|
||||
else:
|
||||
h5file = h5py.File(str(group_or_filename), 'r')
|
||||
need_to_close = True
|
||||
|
||||
# Make sure version matches
|
||||
if 'version' in h5file.attrs:
|
||||
major, minor = h5file.attrs['version']
|
||||
if major != WMP_VERSION_MAJOR:
|
||||
raise IOError(
|
||||
raise DataError(
|
||||
'WMP data format uses version {}. {} whereas your '
|
||||
'installation of the OpenMC Python API expects version '
|
||||
'{}.x.'.format(major, minor, WMP_VERSION_MAJOR))
|
||||
else:
|
||||
raise IOError(
|
||||
raise DataError(
|
||||
'WMP data does not indicate a version. Your installation of '
|
||||
'the OpenMC Python API expects version {}.x data.'
|
||||
.format(WMP_VERSION_MAJOR))
|
||||
|
|
@ -382,6 +385,10 @@ class WindowedMultipole(EqualityMixin):
|
|||
raise ValueError("Windowed multipole is only supported for "
|
||||
"curvefits with 3 or more terms.")
|
||||
|
||||
# If HDF5 file was opened here, make sure it gets closed
|
||||
if need_to_close:
|
||||
h5file.close()
|
||||
|
||||
return out
|
||||
|
||||
def _evaluate(self, E, T):
|
||||
|
|
@ -505,7 +512,7 @@ class WindowedMultipole(EqualityMixin):
|
|||
----------
|
||||
path : str
|
||||
Path to write HDF5 file to
|
||||
mode : {'r', r+', 'w', 'x', 'a'}
|
||||
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'}
|
||||
|
|
|
|||
|
|
@ -20,7 +20,7 @@ from .function import Tabulated1D, Sum, ResonancesWithBackground
|
|||
from .grid import linearize, thin
|
||||
from .njoy import make_ace
|
||||
from .product import Product
|
||||
from .reaction import Reaction, _get_photon_products_ace
|
||||
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
|
||||
|
|
@ -345,7 +345,8 @@ class IncidentNeutron(EqualityMixin):
|
|||
# Add grid around each resonance that includes the peak +/- the
|
||||
# width times each value in _RESONANCE_ENERGY_GRID. Values are
|
||||
# constrained so that points around one resonance don't overlap
|
||||
# with points around another. This algorithm is from Fudge.
|
||||
# with points around another. This algorithm is from Fudge
|
||||
# (https://doi.org/10.1063/1.1945057).
|
||||
energies = []
|
||||
for e, g, e_lower, e_upper in zip(e_peak, gamma, e_mid[:-1],
|
||||
e_mid[1:]):
|
||||
|
|
@ -410,7 +411,7 @@ class IncidentNeutron(EqualityMixin):
|
|||
----------
|
||||
path : str
|
||||
Path to write HDF5 file to
|
||||
mode : {'r', 'r+', 'w', 'x', 'a'}
|
||||
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'}
|
||||
|
|
@ -424,64 +425,62 @@ class IncidentNeutron(EqualityMixin):
|
|||
'originated from an ENDF file.')
|
||||
|
||||
# Open file and write version
|
||||
f = h5py.File(str(path), mode, libver=libver)
|
||||
f.attrs['filetype'] = np.string_('data_neutron')
|
||||
f.attrs['version'] = np.array(HDF5_VERSION)
|
||||
with h5py.File(str(path), mode, libver=libver) as f:
|
||||
f.attrs['filetype'] = np.string_('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 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 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 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
|
||||
# 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('reaction_{:03}'.format(rx.mt))
|
||||
rx.to_hdf5(rx_group)
|
||||
rx_group = rxs_group.create_group('reaction_{:03}'.format(rx.mt))
|
||||
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 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 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)
|
||||
|
||||
f.close()
|
||||
# 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):
|
||||
|
|
@ -543,7 +542,7 @@ class IncidentNeutron(EqualityMixin):
|
|||
data.reactions[rx.mt] = rx
|
||||
|
||||
# Read total nu data if available
|
||||
if rx.mt in (18, 19, 20, 21, 38) and 'total_nu' in group:
|
||||
if rx.mt in FISSION_MTS and 'total_nu' in group:
|
||||
tgroup = group['total_nu']
|
||||
rx.derived_products.append(Product.from_hdf5(tgroup))
|
||||
|
||||
|
|
@ -617,18 +616,20 @@ class IncidentNeutron(EqualityMixin):
|
|||
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 reactions (total, absorption, and heating)
|
||||
# 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
|
||||
|
|
|
|||
|
|
@ -676,8 +676,10 @@ class IncidentPhoton(EqualityMixin):
|
|||
"""
|
||||
if isinstance(group_or_filename, h5py.Group):
|
||||
group = group_or_filename
|
||||
need_to_close = False
|
||||
else:
|
||||
h5file = h5py.File(str(group_or_filename), 'r')
|
||||
need_to_close = True
|
||||
|
||||
# Make sure version matches
|
||||
if 'version' in h5file.attrs:
|
||||
|
|
@ -738,6 +740,10 @@ class IncidentPhoton(EqualityMixin):
|
|||
'num_electrons', 'photon_energy'):
|
||||
data.bremsstrahlung[key] = rgroup[key][()]
|
||||
|
||||
# If HDF5 file was opened here, make sure it gets closed
|
||||
if need_to_close:
|
||||
h5file.close()
|
||||
|
||||
return data
|
||||
|
||||
def export_to_hdf5(self, path, mode='a', libver='earliest'):
|
||||
|
|
@ -747,7 +753,7 @@ class IncidentPhoton(EqualityMixin):
|
|||
----------
|
||||
path : str
|
||||
Path to write HDF5 file to
|
||||
mode : {'r', 'r+', 'w', 'x', 'a'}
|
||||
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'}
|
||||
|
|
@ -755,69 +761,69 @@ class IncidentPhoton(EqualityMixin):
|
|||
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.string_('data_photon')
|
||||
if 'version' not in f.attrs:
|
||||
f.attrs['version'] = np.array(HDF5_VERSION)
|
||||
with h5py.File(str(path), mode, libver=libver) as f:
|
||||
# Write filetype and version
|
||||
f.attrs['filetype'] = np.string_('data_photon')
|
||||
if 'version' not in f.attrs:
|
||||
f.attrs['version'] = np.array(HDF5_VERSION)
|
||||
|
||||
group = f.create_group(self.name)
|
||||
group.attrs['Z'] = Z = self.atomic_number
|
||||
group = f.create_group(self.name)
|
||||
group.attrs['Z'] = Z = self.atomic_number
|
||||
|
||||
# Determine union energy grid
|
||||
union_grid = np.array([])
|
||||
for rx in self:
|
||||
union_grid = np.union1d(union_grid, rx.xs.x)
|
||||
group.create_dataset('energy', data=union_grid)
|
||||
# Determine union energy grid
|
||||
union_grid = np.array([])
|
||||
for rx in self:
|
||||
union_grid = np.union1d(union_grid, rx.xs.x)
|
||||
group.create_dataset('energy', data=union_grid)
|
||||
|
||||
# Write cross sections
|
||||
shell_group = group.create_group('subshells')
|
||||
designators = []
|
||||
for mt, rx in self.reactions.items():
|
||||
name, key = _REACTION_NAME[mt]
|
||||
if mt in (502, 504, 515, 517, 522, 525):
|
||||
sub_group = group.create_group(key)
|
||||
elif mt >= 534 and mt <= 572:
|
||||
# Subshell
|
||||
designators.append(key)
|
||||
sub_group = shell_group.create_group(key)
|
||||
# Write cross sections
|
||||
shell_group = group.create_group('subshells')
|
||||
designators = []
|
||||
for mt, rx in self.reactions.items():
|
||||
name, key = _REACTION_NAME[mt]
|
||||
if mt in (502, 504, 515, 517, 522, 525):
|
||||
sub_group = group.create_group(key)
|
||||
elif mt >= 534 and mt <= 572:
|
||||
# Subshell
|
||||
designators.append(key)
|
||||
sub_group = shell_group.create_group(key)
|
||||
|
||||
# Write atomic relaxation
|
||||
if self.atomic_relaxation is not None:
|
||||
if key in self.atomic_relaxation.subshells:
|
||||
self.atomic_relaxation.to_hdf5(sub_group, key)
|
||||
else:
|
||||
continue
|
||||
|
||||
rx.to_hdf5(sub_group, union_grid, Z)
|
||||
|
||||
shell_group.attrs['designators'] = np.array(designators, dtype='S')
|
||||
|
||||
# Write Compton profiles
|
||||
if self.compton_profiles:
|
||||
compton_group = group.create_group('compton_profiles')
|
||||
|
||||
profile = self.compton_profiles
|
||||
compton_group.create_dataset('num_electrons',
|
||||
data=profile['num_electrons'])
|
||||
compton_group.create_dataset('binding_energy',
|
||||
data=profile['binding_energy'])
|
||||
|
||||
# Get electron momentum values
|
||||
compton_group.create_dataset('pz', data=profile['J'][0].x)
|
||||
|
||||
# Create/write 2D array of profiles
|
||||
J = np.array([Jk.y for Jk in profile['J']])
|
||||
compton_group.create_dataset('J', data=J)
|
||||
|
||||
# Write bremsstrahlung
|
||||
if self.bremsstrahlung:
|
||||
brem_group = group.create_group('bremsstrahlung')
|
||||
for key, value in self.bremsstrahlung.items():
|
||||
if key == 'I':
|
||||
brem_group.attrs[key] = value
|
||||
# Write atomic relaxation
|
||||
if self.atomic_relaxation is not None:
|
||||
if key in self.atomic_relaxation.subshells:
|
||||
self.atomic_relaxation.to_hdf5(sub_group, key)
|
||||
else:
|
||||
brem_group.create_dataset(key, data=value)
|
||||
continue
|
||||
|
||||
rx.to_hdf5(sub_group, union_grid, Z)
|
||||
|
||||
shell_group.attrs['designators'] = np.array(designators, dtype='S')
|
||||
|
||||
# Write Compton profiles
|
||||
if self.compton_profiles:
|
||||
compton_group = group.create_group('compton_profiles')
|
||||
|
||||
profile = self.compton_profiles
|
||||
compton_group.create_dataset('num_electrons',
|
||||
data=profile['num_electrons'])
|
||||
compton_group.create_dataset('binding_energy',
|
||||
data=profile['binding_energy'])
|
||||
|
||||
# Get electron momentum values
|
||||
compton_group.create_dataset('pz', data=profile['J'][0].x)
|
||||
|
||||
# Create/write 2D array of profiles
|
||||
J = np.array([Jk.y for Jk in profile['J']])
|
||||
compton_group.create_dataset('J', data=J)
|
||||
|
||||
# Write bremsstrahlung
|
||||
if self.bremsstrahlung:
|
||||
brem_group = group.create_group('bremsstrahlung')
|
||||
for key, value in self.bremsstrahlung.items():
|
||||
if key == 'I':
|
||||
brem_group.attrs[key] = value
|
||||
else:
|
||||
brem_group.create_dataset(key, data=value)
|
||||
|
||||
def _add_bremsstrahlung(self):
|
||||
"""Add the data used in the thick-target bremsstrahlung approximation
|
||||
|
|
@ -837,7 +843,8 @@ class IncidentPhoton(EqualityMixin):
|
|||
}
|
||||
|
||||
filename = os.path.join(os.path.dirname(__file__), 'BREMX.DAT')
|
||||
brem = open(filename, 'r').read().split()
|
||||
with open(filename, 'r') as fh:
|
||||
brem = fh.read().split()
|
||||
|
||||
# Incident electron kinetic energy grid in eV
|
||||
_BREMSSTRAHLUNG['electron_energy'] = np.logspace(3, 9, 200)
|
||||
|
|
|
|||
|
|
@ -31,25 +31,22 @@ class Product(EqualityMixin):
|
|||
delayed neutron precursor). A special value of 'total' is used when the
|
||||
yield represents particles from prompt and delayed sources.
|
||||
particle : str
|
||||
What particle the reaction product is.
|
||||
The particle type of the reaction product
|
||||
yield_ : openmc.data.Function1D
|
||||
Yield of secondary particle in the reaction.
|
||||
|
||||
"""
|
||||
|
||||
def __init__(self, particle='neutron'):
|
||||
self.particle = particle
|
||||
self.decay_rate = 0.0
|
||||
self.emission_mode = 'prompt'
|
||||
self.distribution = []
|
||||
self.applicability = []
|
||||
self.yield_ = Polynomial((1,)) # 0-order polynomial i.e. a constant
|
||||
self.decay_rate = 0.0
|
||||
self.distribution = []
|
||||
self.emission_mode = 'prompt'
|
||||
self.particle = particle
|
||||
self.yield_ = Polynomial((1,)) # 0-order polynomial, i.e., a constant
|
||||
|
||||
def __repr__(self):
|
||||
if isinstance(self.yield_, Real):
|
||||
return "<Product: {}, emission={}, yield={}>".format(
|
||||
self.particle, self.emission_mode, self.yield_)
|
||||
elif isinstance(self.yield_, Tabulated1D):
|
||||
if isinstance(self.yield_, Tabulated1D):
|
||||
if np.all(self.yield_.y == self.yield_.y[0]):
|
||||
return "<Product: {}, emission={}, yield={}>".format(
|
||||
self.particle, self.emission_mode, self.yield_.y[0])
|
||||
|
|
|
|||
|
|
@ -64,6 +64,8 @@ REACTION_NAME.update({i: '(n,3He{})'.format(i - 750) for i in range(750, 799)})
|
|||
REACTION_NAME.update({i: '(n,a{})'.format(i - 800) for i in range(800, 849)})
|
||||
REACTION_NAME.update({i: '(n,2n{})'.format(i - 875) for i in range(875, 891)})
|
||||
|
||||
FISSION_MTS = (18, 19, 20, 21, 38)
|
||||
|
||||
|
||||
def _get_products(ev, mt):
|
||||
"""Generate products from MF=6 in an ENDF evaluation
|
||||
|
|
@ -492,7 +494,7 @@ def _get_activation_products(ev, rx):
|
|||
|
||||
# Determine if file 9/10 are present
|
||||
present = {9: False, 10: False}
|
||||
for i in range(n_states):
|
||||
for _ in range(n_states):
|
||||
if decay_sublib:
|
||||
items = get_cont_record(file_obj)
|
||||
else:
|
||||
|
|
@ -1021,7 +1023,7 @@ class Reaction(EqualityMixin):
|
|||
neutron.yield_ = yield_
|
||||
rx.products.append(neutron)
|
||||
else:
|
||||
assert mt in (18, 19, 20, 21, 38)
|
||||
assert mt in FISSION_MTS
|
||||
rx.products, rx.derived_products = _get_fission_products_ace(ace)
|
||||
|
||||
for p in rx.products:
|
||||
|
|
@ -1126,13 +1128,13 @@ class Reaction(EqualityMixin):
|
|||
|
||||
# Get fission product yields (nu) as well as delayed neutron energy
|
||||
# distributions
|
||||
if mt in (18, 19, 20, 21, 38):
|
||||
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 (18, 19, 20, 21, 38) and product.particle == 'neutron':
|
||||
if mt in FISSION_MTS and product.particle == 'neutron':
|
||||
rx.products[0].applicability = product.applicability
|
||||
rx.products[0].distribution = product.distribution
|
||||
else:
|
||||
|
|
@ -1177,7 +1179,7 @@ class Reaction(EqualityMixin):
|
|||
for dist in neutron.distribution:
|
||||
dist.angle = AngleDistribution.from_endf(ev, mt)
|
||||
|
||||
if mt in (18, 19, 20, 21, 38) and (5, mt) in ev.section:
|
||||
if mt in FISSION_MTS and (5, mt) in ev.section:
|
||||
# For fission reactions,
|
||||
rx.products[0].applicability = neutron.applicability
|
||||
rx.products[0].distribution = neutron.distribution
|
||||
|
|
|
|||
|
|
@ -444,7 +444,7 @@ class ThermalScattering(EqualityMixin):
|
|||
----------
|
||||
path : str
|
||||
Path to write HDF5 file to
|
||||
mode : {'r', r+', 'w', 'x', 'a'}
|
||||
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'}
|
||||
|
|
|
|||
|
|
@ -95,9 +95,11 @@ def replace_missing(product, decay_data):
|
|||
# Iterate until we find an existing nuclide
|
||||
while product not in decay_data:
|
||||
if Z > 98:
|
||||
# Assume alpha decay occurs for Z=99 and above
|
||||
Z -= 2
|
||||
A -= 4
|
||||
else:
|
||||
# Otherwise assume a beta- or beta+
|
||||
if beta_minus:
|
||||
Z += 1
|
||||
else:
|
||||
|
|
@ -314,7 +316,7 @@ class Chain:
|
|||
nuclide.reactions.append(ReactionTuple(
|
||||
name, daughter, q_value, 1.0))
|
||||
|
||||
if any(mt in reactions_available for mt in [18, 19, 20, 21, 38]):
|
||||
if any(mt in reactions_available for mt in openmc.data.FISSION_MTS):
|
||||
if parent in fpy_data:
|
||||
q_value = reactions[parent][18]
|
||||
nuclide.reactions.append(
|
||||
|
|
@ -334,10 +336,10 @@ class Chain:
|
|||
yield_energies = [0.0]
|
||||
|
||||
yield_data = {}
|
||||
for E, table in zip(yield_energies, fpy.independent):
|
||||
for E, yield_table in zip(yield_energies, fpy.independent):
|
||||
yield_replace = 0.0
|
||||
yields = defaultdict(float)
|
||||
for product, y in table.items():
|
||||
for product, y in yield_table.items():
|
||||
# Handle fission products that have no decay data
|
||||
if product not in decay_data:
|
||||
daughter = replace_missing(product, decay_data)
|
||||
|
|
@ -735,26 +737,27 @@ class Chain:
|
|||
rxn_Q = parent.reactions[rxn_index[0]].Q
|
||||
|
||||
# Remove existing reactions
|
||||
|
||||
for ix in reversed(rxn_index):
|
||||
parent.reactions.pop(ix)
|
||||
|
||||
# Add new reactions
|
||||
all_meta = True
|
||||
|
||||
for tgt, br in new_ratios.items():
|
||||
all_meta = all_meta and ("_m" in tgt)
|
||||
for target, br in new_ratios.items():
|
||||
all_meta = all_meta and ("_m" in target)
|
||||
parent.reactions.append(ReactionTuple(
|
||||
reaction, tgt, rxn_Q, br))
|
||||
reaction, target, rxn_Q, br))
|
||||
|
||||
# If branching ratios don't add to unity, add reaction to ground
|
||||
# with remainder of branching ratio
|
||||
if all_meta and sums[parent_name] != 1.0:
|
||||
ground_br = 1.0 - sums[parent_name]
|
||||
ground_tgt = grounds.get(parent_name)
|
||||
if ground_tgt is None:
|
||||
ground_target = grounds.get(parent_name)
|
||||
if ground_target is None:
|
||||
pz, pa, pm = zam(parent_name)
|
||||
ground_tgt = gnd_name(pz, pa + 1, 0)
|
||||
new_ratios[ground_tgt] = ground_br
|
||||
ground_target = gnd_name(pz, pa + 1, 0)
|
||||
new_ratios[ground_target] = ground_br
|
||||
parent.reactions.append(ReactionTuple(
|
||||
reaction, ground_tgt, rxn_Q, ground_br))
|
||||
reaction, ground_target, rxn_Q, ground_br))
|
||||
|
||||
@property
|
||||
def fission_yields(self):
|
||||
|
|
|
|||
|
|
@ -232,9 +232,9 @@ class Nuclide:
|
|||
elem.set('half_life', str(self.half_life))
|
||||
elem.set('decay_modes', str(len(self.decay_modes)))
|
||||
elem.set('decay_energy', str(self.decay_energy))
|
||||
for mode, daughter, br in self.decay_modes:
|
||||
for mode_type, daughter, br in self.decay_modes:
|
||||
mode_elem = ET.SubElement(elem, 'decay')
|
||||
mode_elem.set('type', mode)
|
||||
mode_elem.set('type', mode_type)
|
||||
mode_elem.set('target', daughter or "Nothing")
|
||||
mode_elem.set('branching_ratio', str(br))
|
||||
|
||||
|
|
|
|||
|
|
@ -12,7 +12,7 @@ import numpy as np
|
|||
from . import comm, have_mpi, MPI
|
||||
from .reaction_rates import ReactionRates
|
||||
|
||||
_VERSION_RESULTS = (1, 0)
|
||||
VERSION_RESULTS = (1, 0)
|
||||
|
||||
|
||||
__all__ = ["Results"]
|
||||
|
|
@ -174,8 +174,8 @@ class Results:
|
|||
"""
|
||||
new = Results()
|
||||
new.volume = {lm: self.volume[lm] for lm in local_materials}
|
||||
new.mat_to_ind = dict(zip(
|
||||
local_materials, range(len(local_materials))))
|
||||
new.mat_to_ind = {mat: idx for (idx, mat) in enumerate(local_materials)}
|
||||
|
||||
# Direct transfer
|
||||
direct_attrs = ("time", "k", "power", "nuc_to_ind",
|
||||
"mat_to_hdf5_ind", "proc_time")
|
||||
|
|
@ -228,7 +228,7 @@ class Results:
|
|||
|
||||
# Store concentration mat and nuclide dictionaries (along with volumes)
|
||||
|
||||
handle.attrs['version'] = np.array(_VERSION_RESULTS)
|
||||
handle.attrs['version'] = np.array(VERSION_RESULTS)
|
||||
handle.attrs['filetype'] = np.string_('depletion results')
|
||||
|
||||
mat_list = sorted(self.mat_to_hdf5_ind, key=int)
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
import h5py
|
||||
import numpy as np
|
||||
|
||||
from .results import Results, _VERSION_RESULTS
|
||||
from .results import Results, VERSION_RESULTS
|
||||
from openmc.checkvalue import check_filetype_version, check_value
|
||||
|
||||
|
||||
|
|
@ -30,7 +30,7 @@ class ResultsList(list):
|
|||
New instance of depletion results
|
||||
"""
|
||||
with h5py.File(str(filename), "r") as fh:
|
||||
check_filetype_version(fh, 'depletion results', _VERSION_RESULTS[0])
|
||||
check_filetype_version(fh, 'depletion results', VERSION_RESULTS[0])
|
||||
new = cls()
|
||||
|
||||
# Get number of results stored
|
||||
|
|
@ -68,9 +68,9 @@ class ResultsList(list):
|
|||
|
||||
Returns
|
||||
-------
|
||||
time : numpy.ndarray
|
||||
times : numpy.ndarray
|
||||
Array of times in units of ``time_units``
|
||||
concentration : numpy.ndarray
|
||||
concentrations : numpy.ndarray
|
||||
Concentration of specified nuclide in units of ``nuc_units``
|
||||
|
||||
"""
|
||||
|
|
@ -78,30 +78,30 @@ class ResultsList(list):
|
|||
check_value("nuc_units", nuc_units,
|
||||
{"atoms", "atom/b-cm", "atom/cm3"})
|
||||
|
||||
time = np.empty_like(self, dtype=float)
|
||||
concentration = np.empty_like(self, dtype=float)
|
||||
times = np.empty_like(self, dtype=float)
|
||||
concentrations = np.empty_like(self, dtype=float)
|
||||
|
||||
# Evaluate value in each region
|
||||
for i, result in enumerate(self):
|
||||
time[i] = result.time[0]
|
||||
concentration[i] = result[0, mat, nuc]
|
||||
times[i] = result.time[0]
|
||||
concentrations[i] = result[0, mat, nuc]
|
||||
|
||||
# Unit conversions
|
||||
if time_units == "d":
|
||||
time /= (60 * 60 * 24)
|
||||
times /= (60 * 60 * 24)
|
||||
elif time_units == "h":
|
||||
time /= (60 * 60)
|
||||
times /= (60 * 60)
|
||||
elif time_units == "min":
|
||||
time /= 60
|
||||
times /= 60
|
||||
|
||||
if nuc_units != "atoms":
|
||||
# Divide by volume to get density
|
||||
concentration /= self[0].volume[mat]
|
||||
concentrations /= self[0].volume[mat]
|
||||
if nuc_units == "atom/b-cm":
|
||||
# 1 barn = 1e-24 cm^2
|
||||
concentration *= 1e-24
|
||||
concentrations *= 1e-24
|
||||
|
||||
return time, concentration
|
||||
return times, concentrations
|
||||
|
||||
def get_reaction_rate(self, mat, nuc, rx):
|
||||
"""Get reaction rate in a single material/nuclide over time
|
||||
|
|
@ -125,44 +125,44 @@ class ResultsList(list):
|
|||
|
||||
Returns
|
||||
-------
|
||||
time : numpy.ndarray
|
||||
times : numpy.ndarray
|
||||
Array of times in [s]
|
||||
rate : numpy.ndarray
|
||||
rates : numpy.ndarray
|
||||
Array of reaction rates
|
||||
|
||||
"""
|
||||
time = np.empty_like(self, dtype=float)
|
||||
rate = np.empty_like(self, dtype=float)
|
||||
times = np.empty_like(self, dtype=float)
|
||||
rates = np.empty_like(self, dtype=float)
|
||||
|
||||
# Evaluate value in each region
|
||||
for i, result in enumerate(self):
|
||||
time[i] = result.time[0]
|
||||
rate[i] = result.rates[0].get(mat, nuc, rx) * result[0, mat, nuc]
|
||||
times[i] = result.time[0]
|
||||
rates[i] = result.rates[0].get(mat, nuc, rx) * result[0, mat, nuc]
|
||||
|
||||
return time, rate
|
||||
return times, rates
|
||||
|
||||
def get_eigenvalue(self):
|
||||
"""Evaluates the eigenvalue from a results list.
|
||||
|
||||
Returns
|
||||
-------
|
||||
time : numpy.ndarray
|
||||
times : numpy.ndarray
|
||||
Array of times in [s]
|
||||
eigenvalue : numpy.ndarray
|
||||
eigenvalues : numpy.ndarray
|
||||
k-eigenvalue at each time. Column 0
|
||||
contains the eigenvalue, while column
|
||||
1 contains the associated uncertainty
|
||||
|
||||
"""
|
||||
time = np.empty_like(self, dtype=float)
|
||||
eigenvalue = np.empty((len(self), 2), dtype=float)
|
||||
times = np.empty_like(self, dtype=float)
|
||||
eigenvalues = np.empty((len(self), 2), dtype=float)
|
||||
|
||||
# Get time/eigenvalue at each point
|
||||
for i, result in enumerate(self):
|
||||
time[i] = result.time[0]
|
||||
eigenvalue[i] = result.k[0]
|
||||
times[i] = result.time[0]
|
||||
eigenvalues[i] = result.k[0]
|
||||
|
||||
return time, eigenvalue
|
||||
return times, eigenvalues
|
||||
|
||||
def get_depletion_time(self):
|
||||
"""Return an array of the average time to deplete a material
|
||||
|
|
@ -175,7 +175,7 @@ class ResultsList(list):
|
|||
|
||||
Returns
|
||||
-------
|
||||
times : :class:`numpy.ndarray`
|
||||
times : numpy.ndarray
|
||||
Vector of average time to deplete a single material
|
||||
across all processes and materials.
|
||||
|
||||
|
|
|
|||
|
|
@ -22,6 +22,7 @@ def test_float_endf():
|
|||
assert endf.float_endf('1.+2') == approx(100.0)
|
||||
assert endf.float_endf('-1.+2') == approx(-100.0)
|
||||
assert endf.float_endf(' ') == 0.0
|
||||
assert endf.float_endf('9.876540000000000') == approx(9.87654)
|
||||
|
||||
|
||||
def test_int_endf():
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue