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

Merge branch 'develop' into photon-new

Browse source
This commit is contained in:
amandalund 2018-06-28 16:19:54 -05:00
commit c912cdce70
264 changed files with 18402 additions and 17147 deletions
Download patch file Download diff file Expand all files Collapse all files

6
openmc/data/ace.py
View file

@ -24,7 +24,7 @@ import numpy as np
from openmc.mixin import EqualityMixin
import openmc.checkvalue as cv
from .data import ATOMIC_SYMBOL
from .data import ATOMIC_SYMBOL, gnd_name
from .endf import ENDF_FLOAT_RE
@ -88,9 +88,7 @@ def get_metadata(zaid, metastable_scheme='nndc'):
# Determine name
element = ATOMIC_SYMBOL[Z]
name = '{}{}'.format(element, mass_number)
if metastable > 0:
name += '_m{}'.format(metastable)
name = gnd_name(Z, mass_number, metastable)
return (name, element, Z, mass_number, metastable)

34
openmc/data/data.py
View file

@ -1,10 +1,9 @@
import itertools
from math import sqrt
import os
import re
from warnings import warn
from numpy import sqrt
# Isotopic abundances from Meija J, Coplen T B, et al, "Isotopic compositions
# of the elements 2013 (IUPAC Technical Report)", Pure. Appl. Chem. 88 (3),
@ -136,23 +135,24 @@ ATOMIC_NUMBER = {value: key for key, value in ATOMIC_SYMBOL.items()}
_ATOMIC_MASS = {}
_GND_NAME_RE = re.compile(r'([A-Zn][a-z]*)(\d+)((?:_[em]\d+)?)')
def atomic_mass(isotope):
"""Return atomic mass of isotope in atomic mass units.
Atomic mass data comes from the Atomic Mass Evaluation 2012, published in
Chinese Physics C 36 (2012), 1287--1602.
Atomic mass data comes from the `Atomic Mass Evaluation 2012
<https://www-nds.iaea.org/amdc/ame2012/AME2012-1.pdf>`_.
Parameters
----------
isotope : str
Name of isotope, e.g. 'Pu239'
Name of isotope, e.g., 'Pu239'
Returns
-------
float or None
Atomic mass of isotope in atomic mass units. If the isotope listed does
not have a known atomic mass, None is returned.
float
Atomic mass of isotope in [amu]
"""
if not _ATOMIC_MASS:
@ -183,7 +183,7 @@ def atomic_mass(isotope):
if '_' in isotope:
isotope = isotope[:isotope.find('_')]
return _ATOMIC_MASS.get(isotope.lower())
return _ATOMIC_MASS[isotope.lower()]
def atomic_weight(element):
@ -199,16 +199,19 @@ def atomic_weight(element):
Returns
-------
float or None
Atomic weight of element in atomic mass units. If the element listed does
not exist, None is returned.
float
Atomic weight of element in [amu]
"""
weight = 0.
for nuclide, abundance in NATURAL_ABUNDANCE.items():
if re.match(r'{}\d+'.format(element), nuclide):
weight += atomic_mass(nuclide) * abundance
return None if weight == 0. else weight
if weight > 0.:
return weight
else:
raise ValueError("No naturally-occurring isotopes for element '{}'."
.format(element))
def water_density(temperature, pressure=0.1013):
@ -234,7 +237,7 @@ def water_density(temperature, pressure=0.1013):
Returns
-------
float
Water density in units of [g / cm^3]
Water density in units of [g/cm^3]
"""
@ -352,8 +355,7 @@ def zam(name):
"""
try:
symbol, A, state = re.match(r'([A-Zn][a-z]*)(\d+)((?:_[em]\d+)?)',
name).groups()
symbol, A, state = _GND_NAME_RE.match(name).groups()
except AttributeError:
raise ValueError("'{}' does not appear to be a nuclide name in GND "
"format.".format(name))

5
openmc/data/decay.py
View file

@ -7,10 +7,7 @@ import re
from warnings import warn
import numpy as np
try:
from uncertainties import ufloat, unumpy, UFloat
except ImportError:
ufloat = UFloat = namedtuple('UFloat', ['nominal_value', 'std_dev'])
from uncertainties import ufloat, unumpy, UFloat
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin

9
openmc/data/library.py
View file

@ -23,9 +23,14 @@ class DataLibrary(EqualityMixin):
def __init__(self):
self.libraries = []
def get_by_material(self, value):
def get_by_material(self, name):
"""Return the library dictionary containing a given material.
Parameters
----------
name : str
Name of material, e.g. 'Am241'
Returns
-------
library : dict or None
@ -34,7 +39,7 @@ class DataLibrary(EqualityMixin):
"""
for library in self.libraries:
if value in library['materials']:
if name in library['materials']:
return library
return None

170
openmc/data/multipole.py
View file

@ -24,7 +24,7 @@ _RM_RF = 3 # Residue fission
# Multi-level Breit Wigner indices
_MLBW_RT = 1 # Residue total
_MLBW_RX = 2 # Residue compettitive
_MLBW_RX = 2 # Residue competitive
_MLBW_RA = 3 # Residue absorption
_MLBW_RF = 4 # Residue fission
@ -141,6 +141,12 @@ def _broaden_wmp_polynomials(E, dopp, n):
class WindowedMultipole(EqualityMixin):
"""Resonant cross sections represented in the windowed multipole format.
Parameters
----------
formalism : {'MLBW', 'RM'}
The R-matrix formalism used to reconstruct resonances. Either 'MLBW'
for multi-level Breit Wigner or 'RM' for Reich-Moore.
Attributes
----------
num_l : Integral
@ -195,11 +201,9 @@ class WindowedMultipole(EqualityMixin):
a/E + b/sqrt(E) + c + d sqrt(E) + ...
"""
def __init__(self):
self.num_l = None
self.fit_order = None
self.fissionable = None
self.formalism = None
def __init__(self, formalism):
self._num_l = None
self.formalism = formalism
self.spacing = None
self.sqrtAWR = None
self.start_E = None
@ -218,11 +222,15 @@ class WindowedMultipole(EqualityMixin):
@property
def fit_order(self):
return self._fit_order
return self.curvefit.shape[1] - 1
@property
def fissionable(self):
return self._fissionable
if self.formalism == 'RM':
return self.data.shape[1] == 4
else:
# Assume self.formalism == 'MLBW'
return self.data.shape[1] == 5
@property
def formalism(self):
@ -272,35 +280,10 @@ class WindowedMultipole(EqualityMixin):
def curvefit(self):
return self._curvefit
@num_l.setter
def num_l(self, num_l):
if num_l is not None:
cv.check_type('num_l', num_l, Integral)
cv.check_greater_than('num_l', num_l, 1, equality=True)
cv.check_less_than('num_l', num_l, 4, equality=True)
# There is an if block in _evaluate that assumes num_l <= 4.
self._num_l = num_l
@fit_order.setter
def fit_order(self, fit_order):
if fit_order is not None:
cv.check_type('fit_order', fit_order, Integral)
cv.check_greater_than('fit_order', fit_order, 2, equality=True)
# _broaden_wmp_polynomials assumes the curve fit has at least 3
# terms.
self._fit_order = fit_order
@fissionable.setter
def fissionable(self, fissionable):
if fissionable is not None:
cv.check_type('fissionable', fissionable, bool)
self._fissionable = fissionable
@formalism.setter
def formalism(self, formalism):
if formalism is not None:
cv.check_type('formalism', formalism, str)
cv.check_value('formalism', formalism, ('MLBW', 'RM'))
cv.check_type('formalism', formalism, str)
cv.check_value('formalism', formalism, ('MLBW', 'RM'))
self._formalism = formalism
@spacing.setter
@ -337,9 +320,20 @@ class WindowedMultipole(EqualityMixin):
cv.check_type('data', data, np.ndarray)
if len(data.shape) != 2:
raise ValueError('Multipole data arrays must be 2D')
if data.shape[1] not in (3, 4, 5): # 3 or 4 for RM, 4 or 5 for MLBW
raise ValueError('The second dimension of multipole data arrays'
' must have a length of 3, 4 or 5')
if self.formalism == 'RM':
if data.shape[1] not in (3, 4):
raise ValueError('For the Reich-Moore formalism, '
'data.shape[1] must be 3 or 4. One value for the pole.'
' One each for the total and absorption residues. '
'Possibly one more for a fission residue.')
else:
# Assume self.formalism == 'MLBW'
if data.shape[1] not in (4, 5):
raise ValueError('For the Multi-level Breit-Wigner '
'formalism, data.shape[1] must be 4 or 5. One value '
'for the pole. One each for the total, competitive, '
'and absorption residues. Possibly one more for a '
'fission residue.')
if not np.issubdtype(data.dtype, complex):
raise TypeError('Multipole data arrays must be complex dtype')
self._data = data
@ -363,6 +357,12 @@ class WindowedMultipole(EqualityMixin):
if not np.issubdtype(l_value.dtype, int):
raise TypeError('Multipole l_value arrays must be integer'
' dtype')
self._num_l = len(np.unique(l_value))
else:
self._num_l = None
self._l_value = l_value
@w_start.setter
@ -428,6 +428,7 @@ class WindowedMultipole(EqualityMixin):
format.
"""
if isinstance(group_or_filename, h5py.Group):
group = group_or_filename
else:
@ -442,20 +443,12 @@ class WindowedMultipole(EqualityMixin):
'Python API expects version ' + WMP_VERSION)
group = h5file['nuclide']
out = cls()
# Read scalar values. Note that group['max_w'] is ignored.
length = group['length'].value
windows = group['windows'].value
out.num_l = group['num_l'].value
out.fit_order = group['fit_order'].value
out.fissionable = bool(group['fissionable'].value)
# Read scalars.
if group['formalism'].value == _FORM_MLBW:
out.formalism = 'MLBW'
out = cls('MLBW')
elif group['formalism'].value == _FORM_RM:
out.formalism = 'RM'
out = cls('RM')
else:
raise ValueError('Unrecognized/Unsupported R-matrix formalism')
@ -466,39 +459,36 @@ class WindowedMultipole(EqualityMixin):
# Read arrays.
err = "WMP '{}' array shape is not consistent with the '{}' value"
err = "WMP '{}' array shape is not consistent with the '{}' array shape"
out.data = group['data'].value
if out.data.shape[0] != length:
raise ValueError(err.format('data', 'length'))
out.l_value = group['l_value'].value
if out.l_value.shape[0] != out.data.shape[0]:
raise ValueError(err.format('l_value', 'data'))
out.pseudo_k0RS = group['pseudo_K0RS'].value
if out.pseudo_k0RS.shape[0] != out.num_l:
raise ValueError(err.format('pseudo_k0RS', 'num_l'))
out.l_value = group['l_value'].value
if out.l_value.shape[0] != length:
raise ValueError(err.format('l_value', 'length'))
raise ValueError(err.format('pseudo_k0RS', 'l_value'))
out.w_start = group['w_start'].value
if out.w_start.shape[0] != windows:
raise ValueError(err.format('w_start', 'windows'))
out.w_end = group['w_end'].value
if out.w_end.shape[0] != windows:
raise ValueError(err.format('w_end', 'windows'))
if out.w_end.shape[0] != out.w_start.shape[0]:
raise ValueError(err.format('w_end', 'w_start'))
out.broaden_poly = group['broaden_poly'].value.astype(np.bool)
if out.broaden_poly.shape[0] != windows:
raise ValueError(err.format('broaden_poly', 'windows'))
if out.broaden_poly.shape[0] != out.w_start.shape[0]:
raise ValueError(err.format('broaden_poly', 'w_start'))
out.curvefit = group['curvefit'].value
if out.curvefit.shape[0] != windows:
raise ValueError(err.format('curvefit', 'windows'))
if out.curvefit.shape[1] != out.fit_order + 1:
raise ValueError(err.format('curvefit', 'fit_order'))
if out.curvefit.shape[0] != out.w_start.shape[0]:
raise ValueError(err.format('curvefit', 'w_start'))
# Note that all the file 3 data (group['reactions/MT...']) are ignored.
# _broaden_wmp_polynomials assumes the curve fit has at least 3 terms.
if out.fit_order < 2:
raise ValueError("Windowed multipole is only supported for "
"curvefits with 3 or more terms.")
return out
@ -661,3 +651,47 @@ class WindowedMultipole(EqualityMixin):
fun = np.vectorize(lambda x: self._evaluate(x, T))
return fun(E)
def export_to_hdf5(self, path, libver='earliest'):
"""Export windowed multipole data to an HDF5 file.
Parameters
----------
path : str
Path to write HDF5 file to
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(path, 'w', libver=libver) as f:
f.create_dataset('version', (1, ), dtype='S10')
f['version'][:] = WMP_VERSION.encode('ASCII')
# Make a nuclide group.
g = f.create_group('nuclide')
# Write scalars.
if self.formalism == 'MLBW':
g.create_dataset('formalism',
data=np.array(_FORM_MLBW, dtype=np.int32))
else:
# Assume RM.
g.create_dataset('formalism',
data=np.array(_FORM_RM, dtype=np.int32))
g.create_dataset('spacing', data=np.array(self.spacing))
g.create_dataset('sqrtAWR', data=np.array(self.sqrtAWR))
g.create_dataset('start_E', data=np.array(self.start_E))
g.create_dataset('end_E', data=np.array(self.end_E))
# Write arrays.
g.create_dataset('data', data=self.data)
g.create_dataset('l_value', data=self.l_value)
g.create_dataset('pseudo_K0RS', data=self.pseudo_k0RS)
g.create_dataset('w_start', data=self.w_start)
g.create_dataset('w_end', data=self.w_end)
g.create_dataset('broaden_poly',
data=self.broaden_poly.astype(np.int8))
g.create_dataset('curvefit', data=self.curvefit)

2
openmc/data/neutron.py
View file

@ -40,7 +40,7 @@ class IncidentNeutron(EqualityMixin):
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
:math:`IncidentNeutron.from_endf`.
:meth:`IncidentNeutron.from_endf`.
Parameters
----------