crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Add FissionProductYields class and Decay class

Browse source
This commit is contained in:
Paul Romano 2016-11-07 11:34:09 -06:00
parent 20a628e537
commit 2395f9a96c
6 changed files with 474 additions and 11 deletions
Download patch file Download diff file Expand all files Collapse all files

2
docs/source/conf.py
View file

@ -25,7 +25,7 @@ except ImportError:
MOCK_MODULES = ['numpy', 'numpy.polynomial', 'numpy.polynomial.polynomial',
'h5py', 'pandas', 'opencg']
'h5py', 'pandas', 'uncertainties', 'opencg']
sys.modules.update((mod_name, MagicMock()) for mod_name in MOCK_MODULES)
import numpy as np

2
docs/source/pythonapi/index.rst
View file

@ -373,6 +373,8 @@ Core Classes
openmc.data.CoherentElastic
openmc.data.FissionEnergyRelease
openmc.data.DataLibrary
openmc.data.Decay
openmc.data.FissionProductYields
Core Functions
--------------

1
openmc/data/__init__.py
View file

@ -6,6 +6,7 @@ HDF5_VERSION = (HDF5_VERSION_MAJOR, HDF5_VERSION_MINOR)
from .data import *
from .neutron import *
from .decay import *
from .reaction import *
from .ace import *
from .angle_distribution import *

458
openmc/data/decay.py Normal file
View file

@ -0,0 +1,458 @@
from collections import Iterable, namedtuple
from io import StringIO
from math import log
from numbers import Real
import re
from warnings import warn
from six import string_types
import numpy as np
try:
from uncertainties import ufloat, unumpy, UFloat
except ImportError:
ufloat = UFloat = namedtuple('UFloat', ['nominal_value', 'std_dev'])
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
from .data import ATOMIC_SYMBOL, ATOMIC_NUMBER
from .endf import Evaluation, get_head_record, get_list_record, get_tab1_record
# Gives name and (change in A, change in Z) resulting from decay
_DECAY_MODES = {
0: ('gamma', (0, 0)),
1: ('beta-', (0, 1)),
2: ('ec/beta+', (0, -1)),
3: ('IT', (0, 0)),
4: ('alpha', (-4, -2)),
5: ('n', (-1, 0)),
6: ('sf', None),
7: ('p', (-1, -1)),
8: ('e-', (0, 0)),
9: ('xray', (0, 0)),
10: ('unknown', None)
}
_RADIATION_TYPES = {
0: 'gamma',
1: 'beta-',
2: 'ec/beta+',
4: 'alpha',
5: 'n',
6: 'sf',
7: 'p',
8: 'e-',
9: 'xray',
10: 'anti-neutrino',
11: 'neutrino'
}
def get_decay_modes(value):
"""Return sequence of decay modes given an ENDF RTYP value.
Parameters
----------
value : float
ENDF definition of sequence of decay modes
Returns
-------
list of str
List of successive decays, e.g. ('beta-', 'neutron')
"""
return [_DECAY_MODES[int(x)][0] for x in
str(value).strip('0').replace('.', '')]
class FissionProductYields(EqualityMixin):
"""Independent and cumulative fission product yields.
Parameters
----------
ev_or_filename : str of openmc.data.endf.Evaluation
ENDF fission product yield evaluation to read from. If given as a
string, it is assumed to be the filename for the ENDF file.
Attributes
----------
cumulative : list of dict
Cumulative yields for each tabulated energy. Each item in the list is a
dictionary whose keys are nuclide names and values are cumulative
yields. The i-th dictionary corresponds to the i-th energy.
energies : Iterable of float or None
Energies at which fission product yields are tabulated.
independent : list of dict
Independent yields for each tabulated energy. Each item in the list is a
dictionary whose keys are nuclide names and values are independent
yields. The i-th dictionary corresponds to the i-th energy.
nuclide : dict
Properties of the fissioning nuclide.
"""
def __init__(self, ev_or_filename):
# Define function that can be used to read both independent and
# cumulative yields
def get_yields(file_obj):
# Determine number of energies
n_energy = get_head_record(file_obj)[2]
if n_energy > 1:
energies = np.zeros(n_energy)
else:
energies = None
data = []
for i in range(n_energy):
# Determine i-th energy and number of products
items, values = get_list_record(file_obj)
if n_energy > 1:
energies[i] = items[0]
n_products = items[5]
# Get yields for i-th energy
yields = {}
for j in range(n_products):
Z, A = divmod(int(values[4*j]), 1000)
excited_state = int(values[4*j + 1])
name = ATOMIC_SYMBOL[Z] + str(A)
if excited_state > 0:
name += '_e{}'.format(excited_state)
yield_j = ufloat(values[4*j + 2], values[4*j + 3])
yields[name] = yield_j
data.append(yields)
return energies, data
# Get evaluation if str is passed
if isinstance(ev_or_filename, Evaluation):
ev = ev_or_filename
else:
ev = Evaluation(ev_or_filename)
# Assign basic nuclide properties
self.nuclide = {
'name': ev.gnd_name,
'atomic_number': ev.target['atomic_number'],
'mass_number': ev.target['mass_number'],
'isomeric_state': ev.target['isomeric_state']
}
# Read independent yields
if (8, 454) in ev.section:
file_obj = StringIO(ev.section[8, 454])
self.energies, self.independent = get_yields(file_obj)
# Read cumulative yields
if (8, 459) in ev.section:
file_obj = StringIO(ev.section[8, 459])
energies, self.cumulative = get_yields(file_obj)
assert np.all(energies == self.energies)
@classmethod
def from_endf(cls, filename):
"""Generate fission product yield data from an ENDF evaluation
Parameters
----------
ev_or_filename : str or openmc.data.endf.Evaluation
ENDF fission product yield evaluation to read from. If given as a
string, it is assumed to be the filename for the ENDF file.
Returns
-------
openmc.data.FissionProductYields
Fission product yield data
"""
return cls(filename)
class DecayMode(EqualityMixin):
"""Radioactive decay mode.
Parameters
----------
parent : str
Parent decaying nuclide
modes : list of str
Successive decay modes
daughter_state : int
Metastable state of the daughter nuclide
energy : uncertainties.UFloat
Total decay energy in eV available in the decay process.
branching_ratio : uncertainties.UFloat
Fraction of the decay of the parent nuclide which proceeds by this mode.
Attributes
----------
branching_ratio : uncertainties.UFloat
Fraction of the decay of the parent nuclide which proceeds by this mode.
daughter : str
Name of daughter nuclide produced from decay
energy : uncertainties.UFloat
Total decay energy in eV available in the decay process.
modes : list of str
Successive decay modes
parent : str
Parent decaying nuclide
"""
def __init__(self, parent, modes, daughter_state, energy,
branching_ratio):
self._daughter_state = daughter_state
self.parent = parent
self.modes = modes
self.energy = energy
self.branching_ratio = branching_ratio
def __repr__(self):
return ('<DecayMode: ({}), {} -> {}, {}>'.format(
','.join(self.modes), self.parent, self.daughter,
self.branching_ratio))
@property
def branching_ratio(self):
return self._branching_ratio
@property
def daughter(self):
# Determine atomic number and mass number of parent
symbol, A = re.match(r'([A-Zn][a-z]*)(\d+)', self.parent).groups()
A = int(A)
Z = ATOMIC_NUMBER[symbol]
# Process changes
for mode in self.modes:
for name, changes in _DECAY_MODES.values():
if name == mode:
if changes is not None:
delta_A, delta_Z = changes
A += delta_A
Z += delta_Z
if self._daughter_state > 0:
return '{}{}_m{}'.format(ATOMIC_SYMBOL[Z], A, self._daughter_state)
else:
return '{}{}'.format(ATOMIC_SYMBOL[Z], A)
@property
def energy(self):
return self._energy
@property
def modes(self):
return self._modes
@property
def parent(self):
return self._parent
@branching_ratio.setter
def branching_ratio(self, branching_ratio):
cv.check_type('branching ratio', branching_ratio, UFloat)
cv.check_greater_than('branching ratio',
branching_ratio.nominal_value, 0.0, True)
if branching_ratio.nominal_value == 0.0:
warn('Decay mode {} of parent {} has a zero branching ratio.'
.format(self.modes, self.parent))
cv.check_greater_than('branching ratio uncertainty',
branching_ratio.std_dev, 0.0, True)
self._branching_ratio = branching_ratio
@energy.setter
def energy(self, energy):
cv.check_type('decay energy', energy, UFloat)
cv.check_greater_than('decay energy', energy.nominal_value, 0.0, True)
cv.check_greater_than('decay energy uncertainty',
energy.std_dev, 0.0, True)
self._energy = energy
@modes.setter
def modes(self, modes):
cv.check_type('decay modes', modes, Iterable, string_types)
self._modes = modes
@parent.setter
def parent(self, parent):
cv.check_type('parent nuclide', parent, string_types)
self._parent = parent
class Decay(EqualityMixin):
"""Radioactive decay data.
Parameters
----------
ev_or_filename : str of openmc.data.endf.Evaluation
ENDF radioactive decay data evaluation to read from. If given as a
string, it is assumed to be the filename for the ENDF file.
Attributes
----------
average_energies : dict
Average decay energies in eV of each type of radiation for decay heat
applications.
decay_constant : uncertainties.UFloat
Decay constant in inverse seconds.
half_life : uncertainties.UFloat
Half-life of the decay in seconds.
modes : list
Decay mode information for each mode of decay.
nuclide : dict
Dictionary describing decaying nuclide with keys 'name',
'excited_state', 'mass', 'stable', 'spin', and 'parity'.
spectra : dict
Resulting radiation spectra for each radiation type.
"""
def __init__(self, ev_or_filename):
# Get evaluation if str is passed
if isinstance(ev_or_filename, Evaluation):
ev = ev_or_filename
else:
ev = Evaluation(ev_or_filename)
#assert ev.info['sublibrary'] == 'Radioactive decay data'
file_obj = StringIO(ev.section[8, 457])
self.nuclide = {}
self.modes = []
self.spectra = {}
self.average_energies = {}
# Get head record
items = get_head_record(file_obj)
Z, A = divmod(items[0], 1000)
metastable = items[3]
self.nuclide['atomic_number'] = Z
self.nuclide['mass_number'] = A
self.nuclide['isomeric_state'] = metastable
if metastable > 0:
self.nuclide['name'] = '{}{}_m{}'.format(ATOMIC_SYMBOL[Z], A,
metastable)
else:
self.nuclide['name'] = '{}{}'.format(ATOMIC_SYMBOL[Z], A)
self.nuclide['mass'] = items[1] # AWR
self.nuclide['excited_state'] = items[2] # State of the original nuclide
self.nuclide['stable'] = (items[4] == 1) # Nucleus stability flag
# Determine if radioactive/stable
if not self.nuclide['stable']:
NSP = items[5] # Number of radiation types
# Half-life and decay energies
items, values = get_list_record(file_obj)
self.half_life = ufloat(items[0], items[1])
NC = items[4]//2
pairs = [x for x in zip(values[::2], values[1::2])]
ex = self.average_energies
ex['light'] = ufloat(*pairs[0])
ex['electromagnetic'] = ufloat(*pairs[1])
ex['heavy'] = ufloat(*pairs[2])
if NC == 17:
ex['beta-'] = ufloat(*pairs[3])
ex['beta+'] = ufloat(*pairs[4])
ex['auger'] = ufloat(*pairs[5])
ex['conversion'] = ufloat(*pairs[6])
ex['gamma'] = ufloat(*pairs[7])
ex['xray'] = ufloat(*pairs[8])
ex['Bremsstrahlung'] = ufloat(*pairs[9])
ex['annihilation'] = ufloat(*pairs[10])
ex['alpha'] = ufloat(*pairs[11])
ex['recoil'] = ufloat(*pairs[12])
ex['SF'] = ufloat(*pairs[13])
ex['neutron'] = ufloat(*pairs[14])
ex['proton'] = ufloat(*pairs[15])
ex['neutrino'] = ufloat(*pairs[16])
items, values = get_list_record(file_obj)
spin = items[0]
if spin == -77.777:
self.nuclide['spin'] = None
else:
self.nuclide['spin'] = spin
self.nuclide['parity'] = items[1] # Parity of the nuclide
# Decay mode information
n_modes = items[5] # Number of decay modes
for i in range(n_modes):
decay_type = get_decay_modes(values[6*i])
isomeric_state = int(values[6*i + 1])
energy = ufloat(*values[6*i + 2:6*i + 4])
branching_ratio = ufloat(*values[6*i + 4:6*(i + 1)])
mode = DecayMode(self.nuclide['name'], decay_type, isomeric_state,
energy, branching_ratio)
self.modes.append(mode)
discrete_type = {0.0: None, 1.0: 'allowed', 2.0: 'first-forbidden',
3.0: 'second-forbidden', 4.0: 'third-forbidden',
5.0: 'fourth-forbidden', 6.0: 'fifth-forbidden'}
# Read spectra
for i in range(NSP):
spectrum = {}
items, values = get_list_record(file_obj)
# Decay radiation type
spectrum['type'] = _RADIATION_TYPES[items[1]]
# Continuous spectrum flag
spectrum['continuous_flag'] = {0: 'discrete', 1: 'continuous',
2: 'both'}[items[2]]
spectrum['discrete_normalization'] = ufloat(*values[0:2])
spectrum['energy_average'] = ufloat(*values[2:4])
spectrum['continuous_normalization'] = ufloat(*values[4:6])
NER = items[5] # Number of tabulated discrete energies
if not spectrum['continuous_flag'] == 'continuous':
# Information about discrete spectrum
spectrum['discrete'] = []
for j in range(NER):
items, values = get_list_record(file_obj)
di = {}
di['energy'] = ufloat(*items[0:2])
di['from_mode'] = get_decay_modes(values[0])
di['type'] = discrete_type[values[1]]
di['intensity'] = ufloat(*values[2:4])
if spectrum['type'] == 'ec/beta+':
di['positron_intensity'] = ufloat(*values[4:6])
elif spectrum['type'] == 'gamma':
di['internal_pair'] = ufloat(*values[4:6])
if len(values) >= 8:
di['total_internal_conversion'] = ufloat(*values[6:8])
if len(values) == 12:
di['k_shell_conversion'] = ufloat(*values[8:10])
di['l_shell_conversion'] = ufloat(*values[10:12])
spectrum['discrete'].append(di)
if not spectrum['continuous_flag'] == 'discrete':
# Read continuous spectrum
ci = {}
params, ci['probability'] = get_tab1_record(file_obj)
ci['type'] = get_decay_modes(params[0])
# Read covariance (Ek, Fk) table
LCOV = params[3]
if LCOV != 0:
items, values = get_list_record(file_obj)
ci['covariance_lb'] = items[3]
ci['covariance'] = zip(values[0::2], values[1::2])
spectrum['continuous'] = ci
# Add spectrum to dictionary
self.spectra[spectrum['type']] = spectrum
else:
items, values = get_list_record(file_obj)
items, values = get_list_record(file_obj)
self.nuclide['spin'] = items[0]
self.nuclide['parity'] = items[1]
@property
def decay_constant(self):
return log(2.)/self.half_life

21
openmc/data/endf.py
View file

@ -17,6 +17,7 @@ from collections import OrderedDict, Iterable
import numpy as np
from numpy.polynomial.polynomial import Polynomial
from .data import ATOMIC_SYMBOL
from .function import Tabulated1D, INTERPOLATION_SCHEME
from openmc.stats.univariate import Uniform, Tabular, Legendre
@ -47,16 +48,6 @@ SUM_RULES = {1: [2, 3],
_ENDF_FLOAT_RE = re.compile(r'([\s\-\+]?\d*\.\d+)([\+\-]\d+)')
def radiation_type(value):
p = {0: 'gamma', 1: 'beta-', 2: 'ec/beta+', 3: 'IT',
4: 'alpha', 5: 'neutron', 6: 'sf', 7: 'proton',
8: 'e-', 9: 'xray', 10: 'unknown'}
if value % 1.0 == 0:
return p[int(value)]
else:
return (p[int(value)], p[int(10*value % 10)])
def float_endf(s):
"""Convert string of floating point number in ENDF to float.
@ -396,6 +387,16 @@ class Evaluation(object):
mod = 0
self.reaction_list.append((mf, mt, nc, mod))
@property
def gnd_name(self):
symbol = ATOMIC_SYMBOL[self.target['atomic_number']]
A = self.target['mass_number']
m = self.target['isomeric_state']
if m > 0:
return '{}{}_m{}'.format(symbol, A, m)
else:
return '{}{}'.format(symbol, A)
class Tabulated2D(object):
"""Metadata for a two-dimensional function.

1
setup.py
View file

@ -43,6 +43,7 @@ if have_setuptools:
# Optional dependencies
'extras_require': {
'decay': ['uncertainties'],
'pandas': ['pandas>=0.17.0'],
'sparse' : ['scipy'],
'vtk': ['vtk', 'silomesh'],