diff --git a/docs/source/conf.py b/docs/source/conf.py index 672217c84..4e21ec20b 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -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 diff --git a/docs/source/developers.rst b/docs/source/developers.rst index d82d89744..0bdae2416 100644 --- a/docs/source/developers.rst +++ b/docs/source/developers.rst @@ -13,6 +13,7 @@ Active development of the OpenMC Monte Carlo code is currently led by: * `Jon Walsh `_ * `Sterling Harper `_ * `Will Boyd `_ +* `Samuel Shaner `_ * `Benoit Forget `_ * `Kord Smith `_ -* `Andrew Siegel `_ +* `Andrew Siegel `_ diff --git a/docs/source/devguide/workflow.rst b/docs/source/devguide/workflow.rst index 2f36436ca..a53bd114b 100644 --- a/docs/source/devguide/workflow.rst +++ b/docs/source/devguide/workflow.rst @@ -238,4 +238,4 @@ from your private repository into a public fork. .. _paid plan: https://github.com/plans .. _Bitbucket: https://bitbucket.org .. _ctest: http://www.cmake.org/cmake/help/v2.8.12/ctest.html -.. _NNDC: http://http://www.nndc.bnl.gov/endf/b7.1/acefiles.html +.. _NNDC: http://www.nndc.bnl.gov/endf/b7.1/acefiles.html diff --git a/docs/source/methods/geometry.rst b/docs/source/methods/geometry.rst index c2e7270b6..36c252bb1 100644 --- a/docs/source/methods/geometry.rst +++ b/docs/source/methods/geometry.rst @@ -384,7 +384,7 @@ x - x_0`, :math:`\bar{y} = y - y_0`, and :math:`\bar{z} = z - z_0`. We then have Expanding equation :eq:`dist-xcone-1` and rearranging terms, we obtain .. math:: - :label: dist-xcylinder-2 + :label: dist-xcone-2 (v^2 + w^2 - R^2u^2) d^2 + 2 (\bar{y}v + \bar{z}w - R^2\bar{x}u) d + (\bar{y}^2 + \bar{z}^2 - R^2\bar{x}^2) = 0 @@ -392,7 +392,7 @@ Expanding equation :eq:`dist-xcone-1` and rearranging terms, we obtain Defining the terms .. math:: - :label: dist-quadric-terms + :label: dist-xcone-terms a = v^2 + w^2 - R^2u^2 @@ -896,4 +896,4 @@ Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0`. Thus, the gradient to the surface is .. _surfaces: http://en.wikipedia.org/wiki/Surface .. _MCNP: http://mcnp.lanl.gov .. _Serpent: http://montecarlo.vtt.fi -.. _Monte Carlo Performance benchmark: https://github.com/paulromano/benchmarks/tree/master/mc-performance/openmc +.. _Monte Carlo Performance benchmark: https://github.com/mit-crpg/benchmarks/tree/master/mc-performance/openmc diff --git a/docs/source/publications.rst b/docs/source/publications.rst index c5b29b194..097cd43d0 100644 --- a/docs/source/publications.rst +++ b/docs/source/publications.rst @@ -82,6 +82,10 @@ Coupling and Multi-physics Geometry -------- +- Derek M. Lax, "Memory efficient indexing algorithm for physical properties in + OpenMC," S. M. Thesis, Massachusetts Institute of Technology + (2015). ``_ + - Derek Lax, William Boyd, Nicholas Horelik, Benoit Forget, and Kord Smith, "A memory efficient algorithm for classifying unique regions in constructive solid geometries," *Proc. PHYSOR*, Kyoto, Japan, Sep. 28--Oct. 3 (2014). @@ -106,6 +110,10 @@ Miscellaneous triggers for the OpenMC Monte Carlo code," *Trans. Am. Nucl. Soc.*, **112**, 637-640 (2015). +- Kyungkwan Noh and Deokjung Lee, "Whole Core Analysis using OpenMC Monte Carlo + Code," *Trans. Kor. Nucl. Soc. Autumn Meeting*, Gyeongju, Korea, + Oct. 24-25, 2013. + - Timothy P. Burke, Brian C. Kiedrowski, and William R. Martin, "Flux and Reaction Rate Kernel Density Estimators in OpenMC," *Trans. Am. Nucl. Soc.*, **109**, 683-686 (2013). @@ -139,7 +147,7 @@ Doppler Broadening - Colin Josey, Pablo Ducru, Benoit Forget, and Kord Smith, "Windowed multipole for cross section Doppler broadening," *J. Comput. Phys.*, **307**, 715-727 - (2016). ``_ + (2016). ``_ - Jonathan A. Walsh, Benoit Forget, Kord S. Smith, and Forrest B. Brown, "On-the-fly Doppler Broadening of Unresolved Resonance Region Cross Sections @@ -169,6 +177,11 @@ Nuclear Data - Paul K. Romano and Sterling M. Harper, "Nuclear data processing capabilities in OpenMC", *Proc. Nuclear Data*, Sep. 11-16, 2016. +- Jonathan A. Walsh, Benoit Froget, Kord S. Smith, and Forrest B. Brown, + "Neutron Cross Section Processing Methods for Improved Integral Benchmarking + of Unresolved Resonance Region Evaluations," *Eur. Phys. J. Web Conf.* + **111**, 06001 (2016). ``_ + - Jonathan A. Walsh, Paul K. Romano, Benoit Forget, and Kord S. Smith, "Optimizations of the energy grid search algorithm in continuous-energy Monte Carlo particle transport codes", *Comput. Phys. Commun.*, **196**, 134-142 @@ -213,10 +226,19 @@ Parallelism memory subsystem on Monte Carlo code performance," *Proc. Joint Int. Conf. M&C+SNA+MC*, Nashville, Tennessee, Apr. 19--23 (2015). +- Hajime Fujita, Nan Dun, Aiman Fang, Zachary A. Rubinstein, Ziming Zheng, Kamil + Iskra, Jeff Hammonds, Anshu Dubey, Pavan Balaji, and Andrew A. Chien, "Using + Global View Resilience (GVR) to add Resilience to Exascale Applications," + *Proc. Supercomputing*, New Orleans, Louisiana, Nov. 16--21, 2014. + - Nicholas Horelik, Benoit Forget, Kord Smith, and Andrew Siegel, "Domain decomposition and terabyte tallies with the OpenMC Monte Carlo neutron transport code," *Proc. PHYSOR*, Kyoto Japan, Sep. 28--Oct. 3 (2014). +- John R. Tramm, Andrew R. Siegel, Tanzima Islam, and Martin Schulz, "XSBench -- + the development and verification of a performance abstraction for Monte Carlo + reactor analysis," *Proc. PHYSOR*, Kyoto, Japan, Sep 28--Oct. 3, 2014. + - Nicholas Horelik, Andrew Siegel, Benoit Forget, and Kord Smith, "Monte Carlo domain decomposition for robust nuclear reactor analysis," *Parallel Comput.*, **40**, 646--660 (2014). ``_ @@ -270,6 +292,10 @@ Parallelism Depletion --------- +- Anas Gul, K. S. Chaudri, R. Khan, and M. Azeen, "Development and verification + of LOOP: A Linkage of ORIGEN2.2 and OpenMC," *Ann. Nucl. Energy*, **99**, + 321--327 (2017). ``_ + - Kai Huang, Hongchun Wu, Yunzhao Li, and Liangzhi Cao, "Generalized depletion chain simplification based of significance analysis," *Proc. PHYSOR*, Sun Valley, Idaho, May 1-5, 2016. diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst index e364452b1..fc01b99e7 100644 --- a/docs/source/pythonapi/index.rst +++ b/docs/source/pythonapi/index.rst @@ -373,6 +373,8 @@ Core Classes openmc.data.CoherentElastic openmc.data.FissionEnergyRelease openmc.data.DataLibrary + openmc.data.Decay + openmc.data.FissionProductYields Core Functions -------------- @@ -476,6 +478,7 @@ Functions openmc.data.endf.float_endf openmc.data.endf.get_cont_record + openmc.data.endf.get_evaluations openmc.data.endf.get_head_record openmc.data.endf.get_tab1_record openmc.data.endf.get_tab2_record diff --git a/openmc/clean_xml.py b/openmc/clean_xml.py index 564281a5c..6aaf64c66 100644 --- a/openmc/clean_xml.py +++ b/openmc/clean_xml.py @@ -65,25 +65,25 @@ def sort_xml_elements(tree): tree.extend(sorted_elements) -def clean_xml_indentation(element, level=0): +def clean_xml_indentation(element, level=0, spaces_per_level=4): """ copy and paste from http://effbot.org/zone/elementent-lib.htm#prettyprint it basically walks your tree and adds spaces and newlines so the tree is printed in a nice way """ - i = "\n" + level*" " + i = "\n" + level*spaces_per_level*" " if len(element): if not element.text or not element.text.strip(): - element.text = i + " " + element.text = i + spaces_per_level*" " if not element.tail or not element.tail.strip(): element.tail = i for sub_element in element: - clean_xml_indentation(sub_element, level+1) + clean_xml_indentation(sub_element, level+1, spaces_per_level) if not sub_element.tail or not sub_element.tail.strip(): sub_element.tail = i diff --git a/openmc/data/__init__.py b/openmc/data/__init__.py index 373136538..c361a204d 100644 --- a/openmc/data/__init__.py +++ b/openmc/data/__init__.py @@ -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 * diff --git a/openmc/data/decay.py b/openmc/data/decay.py new file mode 100644 index 000000000..4327b2fc3 --- /dev/null +++ b/openmc/data/decay.py @@ -0,0 +1,485 @@ +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 incident neutron + 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 incident neutron + energy. + nuclide : dict + Properties of the fissioning nuclide. + + Notes + ----- + Neutron fission yields are typically not measured with a monoenergetic + source of neutrons. As such, if the fission yields are given at, e.g., + 0.0253 eV, one should interpret this as meaning that they are derived from a + typical thermal reactor flux spectrum as opposed to a monoenergetic source + at 0.0253 eV. + + """ + 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] + energies = np.zeros(n_energy) + + data = [] + for i in range(n_energy): + # Determine i-th energy and number of products + items, values = get_list_record(file_obj) + 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) + isomeric_state = int(values[4*j + 1]) + name = ATOMIC_SYMBOL[Z] + str(A) + if isomeric_state > 0: + name += '_m{}'.format(isomeric_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, ev_or_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(ev_or_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 (' {}, {}>'.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) + + 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): + if hasattr(self.half_life, 'n'): + return log(2.)/self.half_life + else: + mu, sigma = self.half_life + return ufloat(log(2.)/mu, log(2.)/mu**2*sigma) + + @classmethod + def from_endf(cls, ev_or_filename): + """Generate radioactive decay data from an ENDF evaluation + + Parameters + ---------- + ev_or_filename : str or 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. + + Returns + ------- + openmc.data.Decay + Radioactive decay data + + """ + return cls(ev_or_filename) diff --git a/openmc/data/endf.py b/openmc/data/endf.py index c0049f977..34553ad2d 100644 --- a/openmc/data/endf.py +++ b/openmc/data/endf.py @@ -14,9 +14,11 @@ import os from math import pi from collections import OrderedDict, Iterable +from six import string_types 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 +49,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. @@ -258,14 +250,40 @@ def get_tab2_record(file_obj): return params, Tabulated2D(breakpoints, interpolation) +def get_evaluations(filename): + """Return a list of all evaluations within an ENDF file. + + Parameters + ---------- + filename : str + Path to ENDF-6 formatted file + + Returns + ------- + list + A list of :class:`openmc.data.endf.Evaluation` instances. + + """ + evaluations = [] + with open(filename, 'r') as fh: + while True: + pos = fh.tell() + line = fh.readline() + if line[66:70] == ' -1': + break + fh.seek(pos) + evaluations.append(Evaluation(fh)) + return evaluations + class Evaluation(object): """ENDF material evaluation with multiple files/sections Parameters ---------- - filename : str - Path to ENDF file to read + filename_or_obj : str or file-like + Path to ENDF file to read or an open file positioned at the start of an + ENDF material Attributes ---------- @@ -282,8 +300,11 @@ class Evaluation(object): indicator (MOD). """ - def __init__(self, filename): - fh = open(filename, 'r') + def __init__(self, filename_or_obj): + if isinstance(filename_or_obj, string_types): + fh = open(filename_or_obj, 'r') + else: + fh = filename_or_obj self.section = {} self.info = {} self.target = {} @@ -313,6 +334,7 @@ class Evaluation(object): # If end of material reached, exit loop if MAT == 0: + fh.readline() break section_data = '' @@ -396,6 +418,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. diff --git a/openmc/data/laboratory.py b/openmc/data/laboratory.py index be449b79a..0a8908362 100644 --- a/openmc/data/laboratory.py +++ b/openmc/data/laboratory.py @@ -137,3 +137,6 @@ class LaboratoryAngleEnergy(AngleEnergy): energy_out.append(energy_out_i) return cls(tab2.breakpoints, tab2.interpolation, energy, mu, energy_out) + + def to_hdf5(self, group): + raise NotImplementedError diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index e2f2d98e0..bb976f9fa 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -55,13 +55,14 @@ REACTION_NAME = {1: '(n,total)', 2: '(n,elastic)', 4: '(n,level)', 195: '(n,4n2a)', 196: '(n,4npa)', 197: '(n,3p)', 198: '(n,n3p)', 199: '(n,3n2pa)', 200: '(n,5n2p)', 444: '(n,damage)', 649: '(n,pc)', 699: '(n,dc)', 749: '(n,tc)', 799: '(n,3Hec)', - 849: '(n,ac)'} -REACTION_NAME.update({i: '(n,n{})'.format(i-50) for i in range(50, 91)}) -REACTION_NAME.update({i: '(n,p{})'.format(i-600) for i in range(600, 649)}) -REACTION_NAME.update({i: '(n,d{})'.format(i-650) for i in range(650, 699)}) -REACTION_NAME.update({i: '(n,t{})'.format(i-700) for i in range(700, 749)}) -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)}) + 849: '(n,ac)', 891: '(n,2nc)'} +REACTION_NAME.update({i: '(n,n{})'.format(i - 50) for i in range(50, 91)}) +REACTION_NAME.update({i: '(n,p{})'.format(i - 600) for i in range(600, 649)}) +REACTION_NAME.update({i: '(n,d{})'.format(i - 650) for i in range(650, 699)}) +REACTION_NAME.update({i: '(n,t{})'.format(i - 700) for i in range(700, 749)}) +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)}) def _get_products(ev, mt): diff --git a/openmc/filter.py b/openmc/filter.py index a7a8b934a..77e409c40 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1,8 +1,7 @@ -from abc import ABCMeta, abstractproperty +from abc import ABCMeta from collections import Iterable, OrderedDict import copy from numbers import Real, Integral -import sys from xml.etree import ElementTree as ET from six import add_metaclass @@ -832,7 +831,7 @@ class RealFilter(Filter): def merge(self, other): if not self.can_merge(other): msg = 'Unable to merge "{0}" with "{1}" ' \ - 'filters'.format(self.type, other.type) + 'filters'.format(type(self), type(other)) raise ValueError(msg) # Merge unique filter bins @@ -933,7 +932,7 @@ class EnergyFilter(RealFilter): if bins[index] < bins[index-1]: msg = 'Unable to add bin edges "{0}" to a "{1}" Filter ' \ 'since they are not monotonically ' \ - 'increasing'.format(bins, self.type) + 'increasing'.format(bins, type(self)) raise ValueError(msg) def get_pandas_dataframe(self, data_size, **kwargs): @@ -1328,7 +1327,7 @@ class MuFilter(RealFilter): if bins[index] < bins[index-1]: msg = 'Unable to add bin edges "{0}" to a "{1}" Filter ' \ 'since they are not monotonically ' \ - 'increasing'.format(bins, self.type) + 'increasing'.format(bins, type(self)) raise ValueError(msg) def get_pandas_dataframe(self, data_size, **kwargs): @@ -1431,7 +1430,7 @@ class PolarFilter(RealFilter): if bins[index] < bins[index-1]: msg = 'Unable to add bin edges "{0}" to a "{1}" Filter ' \ 'since they are not monotonically ' \ - 'increasing'.format(bins, self.type) + 'increasing'.format(bins, type(self)) raise ValueError(msg) def get_pandas_dataframe(self, data_size, **kwargs): @@ -1479,7 +1478,8 @@ class PolarFilter(RealFilter): hi_bins = np.tile(hi_bins, tile_factor) # Add the new angle columns to the DataFrame. - df.loc[:, self.type + ' low'] = lo_bins + df.loc[:, 'polar low'] = lo_bins + df.loc[:, 'polar high'] = hi_bins return df @@ -1533,7 +1533,7 @@ class AzimuthalFilter(RealFilter): if bins[index] < bins[index-1]: msg = 'Unable to add bin edges "{0}" to a "{1}" Filter ' \ 'since they are not monotonically ' \ - 'increasing'.format(bins, self.type) + 'increasing'.format(bins, type(self)) raise ValueError(msg) def get_pandas_dataframe(self, data_size, distribcell_paths=True): @@ -1581,7 +1581,8 @@ class AzimuthalFilter(RealFilter): hi_bins = np.tile(hi_bins, tile_factor) # Add the new angle columns to the DataFrame. - df.loc[:, self.type + ' low'] = lo_bins + df.loc[:, 'azimuthal low'] = lo_bins + df.loc[:, 'azimuthal high'] = hi_bins return df diff --git a/setup.py b/setup.py index 0885c28a8..befb9c0d2 100755 --- a/setup.py +++ b/setup.py @@ -43,6 +43,7 @@ if have_setuptools: # Optional dependencies 'extras_require': { + 'decay': ['uncertainties'], 'pandas': ['pandas>=0.17.0'], 'sparse' : ['scipy'], 'vtk': ['vtk', 'silomesh'],