From e19291118e1398f399e780db5c6f661c265a9fc7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 9 Jun 2016 14:52:58 -0500 Subject: [PATCH] Refactor use of ACE tables to generate neutron/thermal scattering data --- data/convert_hdf5.py | 96 - data/cross_sections.xml | 2076 --------------------- data/cross_sections_ascii.xml | 2066 --------------------- data/cross_sections_nndc.xml | 870 --------- data/cross_sections_serpent.xml | 2707 ---------------------------- data/get_nndc_data.py | 50 +- docs/source/index.rst | 6 +- docs/source/methods/physics.rst | 410 ++--- docs/source/pythonapi/index.rst | 8 +- openmc/data/__init__.py | 3 + openmc/data/ace.py | 1573 +--------------- openmc/data/angle_distribution.py | 68 +- openmc/data/angle_energy.py | 70 + openmc/data/container.py | 39 + openmc/data/correlated.py | 117 +- openmc/data/energy_distribution.py | 227 ++- openmc/data/kalbach_mann.py | 94 +- openmc/data/library.py | 47 + openmc/data/nbody.py | 24 + openmc/data/neutron.py | 405 +++++ openmc/data/reaction.py | 512 ++++++ openmc/data/thermal.py | 305 +++- openmc/data/urr.py | 39 + openmc/mgxs/library.py | 4 +- openmc/mgxs/mgxs.py | 70 +- openmc/mgxs_library.py | 2 +- openmc/nuclide.py | 6 +- openmc/tallies.py | 14 +- scripts/openmc-ace-to-hdf5 | 139 ++ scripts/openmc-ascii-to-binary | 13 - scripts/openmc-update-inputs | 66 + scripts/openmc-xsdata-to-xml | 148 -- scripts/openmc-xsdir-to-xml | 288 --- src/reaction_header.F90 | 2 +- src/relaxng/cross_sections.rnc | 29 +- src/relaxng/cross_sections.rng | 112 +- 36 files changed, 2467 insertions(+), 10238 deletions(-) delete mode 100644 data/convert_hdf5.py delete mode 100644 data/cross_sections.xml delete mode 100644 data/cross_sections_ascii.xml delete mode 100644 data/cross_sections_nndc.xml delete mode 100644 data/cross_sections_serpent.xml create mode 100644 openmc/data/library.py create mode 100644 openmc/data/neutron.py create mode 100644 openmc/data/reaction.py create mode 100755 scripts/openmc-ace-to-hdf5 delete mode 100755 scripts/openmc-ascii-to-binary delete mode 100755 scripts/openmc-xsdata-to-xml delete mode 100755 scripts/openmc-xsdir-to-xml diff --git a/data/convert_hdf5.py b/data/convert_hdf5.py deleted file mode 100644 index 3f5fcfdb7..000000000 --- a/data/convert_hdf5.py +++ /dev/null @@ -1,96 +0,0 @@ -#!/usr/bin/env python - -import glob -import os -from xml.dom.minidom import getDOMImplementation - -import openmc.data.ace - - -if not os.path.isdir('nndc_hdf5'): - os.mkdir('nndc_hdf5') - -nndc_files = glob.glob('nndc/293.6K/*.ace') -nndc_thermal_files = glob.glob('nndc/tsl/*.acer') - -thermal_names = {'al': 'c_Al27', - 'be': 'c_Be', - 'bebeo': 'c_Be_in_BeO', - 'benzine': 'c_Benzine', - 'dd2o': 'c_D_in_D2O', - 'fe': 'c_Fe56', - 'graphite': 'c_Graphite', - 'hch2': 'c_H_in_CH2', - 'hh2o': 'c_H_in_H2O', - 'hzrh': 'c_H_in_ZrH', - 'lch4': 'c_liquid_CH4', - 'obeo': 'c_O_in_BeO', - 'orthod': 'c_ortho_D', - 'orthoh': 'c_ortho_H', - 'ouo2': 'c_O_in_UO2', - 'parad': 'c_para_D', - 'parah': 'c_para_H', - 'sch4': 'c_solid_CH4', - 'uuo2': 'c_U_in_UO2', - 'zrzrh': 'c_Zr_in_ZrH'} - -impl = getDOMImplementation() -doc = impl.createDocument(None, "cross_sections", None) -doc_root = doc.documentElement - -for f in sorted(nndc_files): - print('Converting {}...'.format(f)) - - # Deterine output file name - dirname, basename = os.path.split(f) - root, ext = os.path.splitext(basename) - outfile = os.path.join('nndc_hdf5', root + '.h5') - if os.path.exists(outfile): - os.remove(outfile) - - # Determine elemental symbol, mass number and metastable state - element, mass_number, temp = basename.split('_') - metastable = int(mass_number[-1]) if 'm' in mass_number else 0 - mass_number = int(mass_number[:3]) - - # Parse ACE file, create HDF5 file - t = openmc.data.ace.get_table(f) - t.export_to_hdf5(outfile, element, mass_number, metastable) - xs = t.name.split('.')[1] - if metastable > 0: - name = "{}{}_m{}.{}".format(element, mass_number, metastable, xs) - else: - name = "{}{}.{}".format(element, mass_number, xs) - - # Add entry to XML listing - libraryNode = doc.createElement("library") - libraryNode.setAttribute("path", root + '.h5') - libraryNode.setAttribute("materials", name) - libraryNode.setAttribute("type", "neutron") - doc_root.appendChild(libraryNode) - -for f in sorted(nndc_thermal_files): - print('Converting {}...'.format(f)) - - # Deterine output file name - dirname, basename = os.path.split(f) - root, ext = os.path.splitext(basename) - outfile = os.path.join('nndc_hdf5', root + '.h5') - if os.path.exists(outfile): - os.remove(outfile) - - # Parse ACE file, create HDF5 file - t = openmc.data.ace.get_table(f) - t.export_to_hdf5(outfile, thermal_names[root]) - xs = t.name.split('.')[1] - - # Add entry to XML listing - libraryNode = doc.createElement("library") - libraryNode.setAttribute("path", root + '.h5') - libraryNode.setAttribute("materials", thermal_names[root] + '.' + xs) - libraryNode.setAttribute("type", "thermal") - doc_root.appendChild(libraryNode) - -# Write cross_sections.xml -lines = doc.toprettyxml(indent=' ') -open(os.path.join('nndc_hdf5', 'cross_sections.xml'), 'w').write(lines) diff --git a/data/cross_sections.xml b/data/cross_sections.xml deleted file mode 100644 index 180021064..000000000 --- a/data/cross_sections.xml +++ /dev/null @@ -1,2076 +0,0 @@ - - - - /opt/mcnp/data - - - binary - - - 4096 - - - 512 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/data/cross_sections_ascii.xml b/data/cross_sections_ascii.xml deleted file mode 100644 index b0f083bb2..000000000 --- a/data/cross_sections_ascii.xml +++ /dev/null @@ -1,2066 +0,0 @@ - - - /opt/mcnp/data - ascii - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/data/cross_sections_nndc.xml b/data/cross_sections_nndc.xml deleted file mode 100644 index 3a7a5ca72..000000000 --- a/data/cross_sections_nndc.xml +++ /dev/null @@ -1,870 +0,0 @@ - - - ascii - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/data/cross_sections_serpent.xml b/data/cross_sections_serpent.xml deleted file mode 100644 index 959790309..000000000 --- a/data/cross_sections_serpent.xml +++ /dev/null @@ -1,2707 +0,0 @@ - - - - /opt/serpent/xsdata/endfb7/acedata - - ascii - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/data/get_nndc_data.py b/data/get_nndc_data.py index 844abea81..a0429fba5 100755 --- a/data/get_nndc_data.py +++ b/data/get_nndc_data.py @@ -20,10 +20,6 @@ try: except ImportError: from urllib2 import urlopen -cwd = os.getcwd() -sys.path.insert(0, os.path.join(cwd, '..')) -from openmc.data.ace import ascii_to_binary - baseUrl = 'http://www.nndc.bnl.gov/endf/b7.1/aceFiles/' files = ['ENDF-B-VII.1-neutron-293.6K.tar.gz', 'ENDF-B-VII.1-tsl.tar.gz'] @@ -114,12 +110,6 @@ text = text.replace('6012', '6000', 1) with open(graphite, 'w') as fh: fh.write(text) -# ============================================================================== -# COPY CROSS_SECTIONS.XML - -print('Copying cross_sections_nndc.xml...') -shutil.copyfile('cross_sections_nndc.xml', 'nndc/cross_sections.xml') - # ============================================================================== # PROMPT USER TO DELETE .TAR.GZ FILES @@ -140,44 +130,26 @@ if not response or response.lower().startswith('y'): os.remove(f) # ============================================================================== -# PROMPT USER TO CONVERT ASCII TO BINARY +# PROMPT USER TO GENERATE HDF5 LIBRARY # Ask user to convert if not args.batch: if sys.version_info[0] < 3: - response = raw_input('Convert ACE files to binary? ([y]/n) ') + response = raw_input('Generate HDF5 library? ([y]/n) ') else: - response = input('Convert ACE files to binary? ([y]/n) ') + response = input('Generate HDF5 library? ([y]/n) ') else: response = 'y' # Convert files if requested if not response or response.lower().startswith('y'): + # get a list of all ACE files + ace_files = sorted(glob.glob(os.path.join('nndc', '**', '*.ace*'))) - # get a list of directories - ace_dirs = glob.glob(os.path.join('nndc', '*K')) - ace_dirs += glob.glob(os.path.join('nndc', 'tsl')) + # Ensure 'import openmc.data' works in the openmc-ace-to-xml script + cwd = os.getcwd() + env = os.environ.copy() + env['PYTHONPATH'] = os.path.join(cwd, '..') - # loop around ace directories - for d in ace_dirs: - print('Converting {0}...'.format(d)) - - # get a list of files to convert - ace_files = glob.glob(os.path.join(d, '*.ace*')) - - # convert files - for f in ace_files: - print(' Converting {0}...'.format(os.path.split(f)[1])) - ascii_to_binary(f, f) - - # Change cross_sections.xml file - xs_file = os.path.join('nndc', 'cross_sections.xml') - asc_str = "ascii" - bin_str = " binary \n " - bin_str += " 4096 \n " - bin_str += " 512 " - with open(xs_file) as fh: - text = fh.read() - text = text.replace(asc_str, bin_str) - with open(xs_file, 'w') as fh: - fh.write(text) + subprocess.call(['../scripts/openmc-ace-to-hdf5', '-d', 'nndc_hdf5'] + + ace_files, env=env) diff --git a/docs/source/index.rst b/docs/source/index.rst index 413107c34..8fe680a70 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -5,9 +5,9 @@ The OpenMC Monte Carlo Code OpenMC is a Monte Carlo particle transport simulation code focused on neutron criticality calculations. It is capable of simulating 3D models based on constructive solid geometry with second-order surfaces. OpenMC supports either -continuous-energy or multi-group transport. The continuous-energy -particle interaction data is based on ACE format cross sections, also used -in the MCNP and Serpent Monte Carlo codes. +continuous-energy or multi-group transport. The continuous-energy particle +interaction data is based on a native HDF5 format that can be generated from ACE +files used by the MCNP and Serpent Monte Carlo codes. OpenMC was originally developed by members of the `Computational Reactor Physics Group`_ at the `Massachusetts Institute of Technology`_ starting diff --git a/docs/source/methods/physics.rst b/docs/source/methods/physics.rst index b246584c8..530d4ac41 100644 --- a/docs/source/methods/physics.rst +++ b/docs/source/methods/physics.rst @@ -279,9 +279,9 @@ idiosyncrasies in treating fission. In an eigenvalue calculation, secondary neutrons from fission are only "banked" for use in the next generation rather than being tracked as secondary neutrons from elastic and inelastic scattering would be. On top of this, fission is sometimes broken into first-chance fission, -second-chance fission, etc. An ACE table either lists the partial fission -reactions with secondary energy distributions for each one, or a total fission -reaction with a single secondary energy distribution. +second-chance fission, etc. The nuclear data file either lists the partial +fission reactions with secondary energy distributions for each one, or a total +fission reaction with a single secondary energy distribution. When a fission reaction is sampled in OpenMC (either total fission or, if data exists, first- or second-chance fission), the following algorithm is used to @@ -290,7 +290,7 @@ number of prompt and delayed neutrons must be determined to decide whether the secondary neutrons will be prompt or delayed. This is important because delayed neutrons have a markedly different spectrum from prompt neutrons, one that has a lower average energy of emission. The total number of neutrons emitted -:math:`\nu_t` is given as a function of incident energy in the ACE format. Two +:math:`\nu_t` is given as a function of incident energy in the ENDF format. Two representations exist for :math:`\nu_t`. The first is a polynomial of order :math:`N` with coefficients :math:`c_0,c_1,\dots,c_N`. If :math:`\nu_t` has this format, we can evaluate it at incoming energy :math:`E` by using the equation @@ -347,26 +347,52 @@ provided as group-wise data instead of in a continuous-energy format. In this case, the outgoing energy of the fission neutrons are represented as histograms by way of either the nu-fission matrix or chi vector. ------------------------------------------ -Secondary Angles and Energy Distributions ------------------------------------------ +------------------------------------ +Secondary Angle-Energy Distributions +------------------------------------ Note that this section is specific to continuous-energy mode since the multi-group scattering process has already been described including the secondary energy and angle sampling. -For any reactions with secondary neutrons, it is necessary to sample secondary -angle and energy distributions. This includes elastic and inelastic scattering, -fission, and :math:`(n,xn)` reactions. In some cases, the angle and energy -distributions may be specified separately, and in other cases, they may be -specified as a correlated angle-energy distribution. In the following sections, -we will outline the methods used to sample secondary distributions as well as -how they are used to modify the state of a particle. +For a reaction with secondary products, it is necessary to determine the +outgoing angle and energy of the products. For any reaction other than elastic +and level inelastic scattering, the outgoing energy must be determined based on +tabulated or parameterized data. The `ENDF-6 Format`_ specifies a variety of +ways that the secondary energy distribution can be represented. ENDF File 5 +contains uncorrelated energy distribution whereas ENDF File 6 contains +correlated energy-angle distributions. The ACE format specifies its own +representations based loosely on the formats given in ENDF-6. OpenMC's HDF5 +nuclear data files use a combination of ENDF and ACE distributions; in this +section, we will describe how the outgoing angle and energy of secondary +particles are sampled. + +One of the subtleties in the nuclear data format is the fact that a single +reaction product can have multiple angle-energy distributions. This is mainly +useful for reactions with multiple products of the same type in the exit channel +such as :math:`(n,2n)` or :math:`(n,3n)`. In these types of reactions, each +neutron is emitted corresponding to a different excitation level of the compound +nucleus, and thus in general the neutrons will originate from different energy +distributions. If multiple angle-energy distributions are present, they are +assigned incoming-energy-dependent probabilities that can then be used to +randomly select one. + +Once a distribution has been selected, the procedure for determining the +outgoing angle and energy will depend on the type of the distribution. + +Uncorrelated Angle-Energy Distributions +--------------------------------------- + +The first set of distributions we will look at are uncorrelated angle-energy +distributions, where angle and energy are specified separately. For these +distributions, OpenMC first samples the angular distribution as described +:ref:`sample-angle` and then samples an energy as described in +:ref:`sample-energy`. .. _sample-angle: -Sampling Secondary Angle Distributions --------------------------------------- +Sampling Angular Distributions +++++++++++++++++++++++++++++++ For elastic scattering, it is only necessary to specific a secondary angle distribution since the outgoing energy can be determined analytically. Other @@ -374,15 +400,14 @@ reactions may also have separate secondary angle and secondary energy distributions that are uncorrelated. In these cases, the secondary angle distribution is represented as either -- An Isotropic angular distribution, -- An equiprobable distribution with 32 bins, or +- An isotropic angular distribution, - A tabular distribution. Isotropic Angular Distribution -++++++++++++++++++++++++++++++ +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -In the first case, no data needs to be stored on the ACE table, and the cosine -of the scattering angle is simply calculated as +In the first case, no data is stored in the nuclear data file, and the cosine of +the scattering angle is simply calculated as .. math:: :label: isotropic-angle @@ -392,42 +417,17 @@ of the scattering angle is simply calculated as where :math:`\mu` is the cosine of the scattering angle and :math:`\xi` is a random number sampled uniformly on :math:`[0,1)`. -Equiprobable Angle Bin Distribution -+++++++++++++++++++++++++++++++++++ - -For a 32 equiprobable bin distribution, we select a random number :math:`\xi` to -sample a cosine bin :math:`i` such that - -.. math:: - :label: equiprobable-bin - - i = 1 + \lfloor 32\xi \rfloor. - -The same random number can then also be used to interpolate between neighboring -:math:`\mu` values to get the final scattering cosine: - -.. math:: - :label: equiprobable-cosine - - \mu = \mu_i + (32\xi - i) (\mu_{i+1} - \mu_i) - -where :math:`\mu_i` is the :math:`i`-th scattering cosine. - .. _angle-tabular: Tabular Angular Distribution -++++++++++++++++++++++++++++ +^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -As the `MCNP Manual`_ points out, using an equiprobable bin distribution works -well for high-probability regions of the scattering cosine probability, but for -low-probability regions it is not very accurate. Thus, a more accurate method is -to represent the scattering cosine with a tabular distribution. In this case, we -have a table of cosines and their corresponding values for a probability -distribution function and cumulative distribution function. For each incoming -neutron energy :math:`E_i`, let us call :math:`p_{i,j}` the j-th value in the -probability distribution function and :math:`c_{i,j}` the j-th value in the -cumulative distribution function. We first find the interpolation factor on the -incoming energy grid: +In this case, we have a table of cosines and their corresponding values for a +probability distribution function and cumulative distribution function. For each +incoming neutron energy :math:`E_i`, let us call :math:`p_{i,j}` the j-th value +in the probability distribution function and :math:`c_{i,j}` the j-th value in +the cumulative distribution function. We first find the interpolation factor on +the incoming energy grid: .. math:: :label: interpolation-factor @@ -545,89 +545,11 @@ linear-linear interpolation: .. _sample-energy: -Sampling Secondary Energy and Correlated Angle/Energy Distributions -------------------------------------------------------------------- +Sampling Energy Distributions ++++++++++++++++++++++++++++++ -For a reaction with secondary neutrons, it is necessary to determine the -outgoing energy of the neutrons. For any reaction other than elastic scattering, -the outgoing energy must be determined based on tabulated or parameterized -data. The `ENDF-6 Format`_ specifies a variety of ways that the secondary energy -distribution can be represented. ENDF File 5 contains uncorrelated energy -distribution where ENDF File 6 contains correlated energy-angle -distributions. The ACE format specifies its own representations based loosely on -the formats given in ENDF-6. In this section, we will describe how the outgoing -energy of secondary particles is determined based on each ACE law. - -One of the subtleties in the ACE format is the fact that a single reaction can -have multiple secondary energy distributions. This is mainly useful for -reactions with multiple neutrons in the exit channel such as :math:`(n,2n)` or -:math:`(n,3n)`. In these types of reactions, each neutron is emitted -corresponding to a different excitation level of the compound nucleus, and thus -in general the neutrons will originate from different energy distributions. If -multiple energy distributions are present, they are assigned probabilities that -can then be used to randomly select one. - -Once a secondary energy distribution has been sampled, the procedure for -determining the outgoing energy will depend on which ACE law has been specified -for the data. - -.. _ace-law-1: - -ACE Law 1 - Tabular Equiprobable Energy Bins -++++++++++++++++++++++++++++++++++++++++++++ - -In the tabular equiprobable bin representation, an array of equiprobable -outgoing energy bins is given for a number of incident energies. While the -representation itself is simple, the complexity lies in how one interpolates -between incident as well as outgoing energies on such a table. If one performs -simple interpolation between tables for neighboring incident energies, it is -possible that the resulting energies would violate laws governing the -kinematics, i.e. the outgoing energy may be outside the range of available -energy in the reaction. - -To avoid this situation, the accepted practice is to use a process known as -scaled interpolation [Doyas]_. First, we find the tabulated incident energies -which bound the actual incoming energy of the particle, i.e. find :math:`i` such -that :math:`E_i < E < E_{i+1}` and calculate the interpolation factor :math:`f` -via :eq:`interpolation-factor`. Then, we interpolate between the minimum and -maximum energies of the outgoing energy distributions corresponding to -:math:`E_i` and :math:`E_{i+1}`: - -.. math:: - :label: ace-law-1-minmax - - E_{min} = E_{i,1} + f ( E_{i+1,1} - E_i ) \\ - E_{max} = E_{i,M} + f ( E_{i+1,M} - E_M ) - -where :math:`E_{min}` and :math:`E_{max}` are the minimum and maximum outgoing -energies of a scaled distribution, :math:`E_{i,j}` is the j-th outgoing energy -corresponding to the incoming energy :math:`E_i`, and :math:`M` is the number of -outgoing energy bins. Next, statistical interpolation is performed to choose -between using the outgoing energy distributions corresponding to energy -:math:`E_i` and :math:`E_{i+1}`. Let :math:`\ell` be the chosen table where -:math:`\ell = i` if :math:`\xi_1 > f` and :math:`\ell = i + 1` otherwise, and -:math:`\xi_1` is a random number. Now, we randomly sample an equiprobable -outgoing energy bin :math:`j` and interpolate between successive values on the -outgoing energy distribution: - -.. math:: - :label: ace-law-1-intermediate - - \hat{E} = E_{\ell,j} + \xi_2 (E_{\ell,j+1} - E_{\ell,j}) - -where :math:`\xi_2` is a random number sampled uniformly on :math:`[0,1)`. Since -this outgoing energy may violate reaction kinematics, we then scale it to the -minimum and maximum energies we calculated earlier to get the final outgoing -energy: - -.. math:: - :label: ace-law-1-energy - - E' = E_{min} + \frac{\hat{E} - E_{\ell,1}}{E_{\ell,M} - E_{\ell,1}} - (E_{max} - E_{min}) - -ACE Law 3 - Inelastic Level Scattering -++++++++++++++++++++++++++++++++++++++ +Inelastic Level Scattering +^^^^^^^^^^^^^^^^^^^^^^^^^^ It can be shown (see Foderaro_) that in inelastic level scattering, the outgoing energy of the neutron :math:`E'` can be related to the Q-value of the reaction @@ -640,31 +562,50 @@ and the incoming energy: where :math:`A` is the mass of the target nucleus measured in neutron masses. -.. _ace-law-4: +.. _continuous-tabular: -ACE Law 4 - Continuous Tabular Distribution -+++++++++++++++++++++++++++++++++++++++++++ +Continuous Tabular Distribution +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -This representation is very similar to :ref:`ace-law-1` except that instead of -equiprobable outgoing energy bins, the outgoing energy distribution for each -incoming energy is represented with a probability distribution function. For -each incoming neutron energy :math:`E_i`, let us call :math:`p_{i,j}` the j-th -value in the probability distribution function, :math:`c_{i,j}` the j-th value -in the cumulative distribution function, and :math:`E_{i,j}` the j-th outgoing -energy. +In a continuous tabular distribution, a tabulated energy distribution is +provided for each of a set of incoming energies. While the representation itself +is simple, the complexity lies in how one interpolates between incident as well +as outgoing energies on such a table. If one performs simple interpolation +between tables for neighboring incident energies, it is possible that the +resulting energies would violate laws governing the kinematics, i.e., the +outgoing energy may be outside the range of available energy in the reaction. -We proceed first as we did for ACE Law 1, determining the bounding energies of -the particle's incoming energy such that :math:`E_i < E < E_{i+1}` and -calculating an interpolation factor :math:`f` with equation -:eq:`interpolation-factor`. Next, statistical interpolation is performed to -choose between using the outgoing energy distributions corresponding to energy -:math:`E_i` and :math:`E_{i+1}`. Let :math:`\ell` be the chosen table where -:math:`\ell = i` if :math:`\xi_1 > f` and :math:`\ell = i + 1` otherwise, and -:math:`\xi_1` is a random number. Then, we sample an outgoing energy bin +To avoid this situation, the accepted practice is to use a process known as +scaled interpolation [Doyas]_. First, we find the tabulated incident energies +which bound the actual incoming energy of the particle, i.e., find :math:`i` +such that :math:`E_i < E < E_{i+1}` and calculate the interpolation factor +:math:`f` via :eq:`interpolation-factor`. Then, we interpolate between the +minimum and maximum energies of the outgoing energy distributions corresponding +to :math:`E_i` and :math:`E_{i+1}`: + +.. math:: + :label: continuous-minmax + + E_{min} = E_{i,1} + f ( E_{i+1,1} - E_{i,1} ) \\ + E_{max} = E_{i,M} + f ( E_{i+1,M} - E_{i,M} ) + +where :math:`E_{min}` and :math:`E_{max}` are the minimum and maximum outgoing +energies of a scaled distribution, :math:`E_{i,j}` is the j-th outgoing energy +corresponding to the incoming energy :math:`E_i`, and :math:`M` is the number of +outgoing energy bins. + +Next, statistical interpolation is performed to choose between using the +outgoing energy distributions corresponding to energy :math:`E_i` and +:math:`E_{i+1}`. Let :math:`\ell` be the chosen table where :math:`\ell = i` if +:math:`\xi_1 > f` and :math:`\ell = i + 1` otherwise, and :math:`\xi_1` is a +random number. For each incoming neutron energy :math:`E_i`, let us call +:math:`p_{i,j}` the j-th value in the probability distribution function, +:math:`c_{i,j}` the j-th value in the cumulative distribution function, and +:math:`E_{i,j}` the j-th outgoing energy. We then sample an outgoing energy bin :math:`j` using the cumulative distribution function: .. math:: - :label: ace-law-4-sample-cdf + :label: continuous-sample-cdf c_{\ell,j} < \xi_2 < c_{\ell,j+1} @@ -692,22 +633,22 @@ If linear-linear interpolation is to be used, the outgoing energy on the \right ). Since this outgoing energy may violate reaction kinematics, we then scale it to -minimum and maximum energies interpolated between the neighboring outgoing -energy distributions to get the final outgoing energy: +minimum and maximum energies calculated in equation :eq:`continuous-minmax` to +get the final outgoing energy: .. math:: - :label: ace-law-4-energy + :label: continuous-eout E' = E_{min} + \frac{\hat{E} - E_{\ell,1}}{E_{\ell,M} - E_{\ell,1}} (E_{max} - E_{min}) where :math:`E_{min}` and :math:`E_{max}` are defined the same as in equation -:eq:`ace-law-1-minmax`. +:eq:`continuous-minmax`. .. _maxwell: -ACE Law 7 - Maxwell Fission Spectrum -++++++++++++++++++++++++++++++++++++ +Maxwell Fission Spectrum +^^^^^^^^^^^^^^^^^^^^^^^^ One representation of the secondary energies for neutrons from fission is the so-called Maxwell spectrum. A probability distribution for the Maxwell spectrum @@ -720,7 +661,7 @@ can be written in the form where :math:`E` is the incoming energy of the neutron and :math:`T` is the so-called nuclear temperature, which is a function of the incoming energy of the -neutron. The ACE format contains a list of nuclear temperatures versus incoming +neutron. The ENDF format contains a list of nuclear temperatures versus incoming energies. The nuclear temperature is interpolated between neighboring incoming energies using a specified interpolation law. Once the temperature :math:`T` is determined, we then calculate a candidate outgoing energy based on rule C64 in @@ -740,12 +681,12 @@ interval. The outgoing energy is only accepted if 0 \le E' \le E - U -where :math:`U` is called the restriction energy and is specified on the ACE -table. If the outgoing energy is rejected, it is resampled using equation +where :math:`U` is called the restriction energy and is specified in the ENDF +data. If the outgoing energy is rejected, it is resampled using equation :eq:`maxwell-E-candidate`. -ACE Law 9 - Evaporation Spectrum -++++++++++++++++++++++++++++++++ +Evaporation Spectrum +^^^^^^^^^^^^^^^^^^^^ Evaporation spectra are primarily used in compound nucleus processes where a secondary particle can "evaporate" from the compound nucleus if it has @@ -759,7 +700,7 @@ be written in the form where :math:`E` is the incoming energy of the neutron and :math:`T` is the nuclear temperature, which is a function of the incoming energy of the -neutron. The ACE format contains a list of nuclear temperatures versus incoming +neutron. The ENDF format contains a list of nuclear temperatures versus incoming energies. The nuclear temperature is interpolated between neighboring incoming energies using a specified interpolation law. Once the temperature :math:`T` is determined, we then calculate a candidate outgoing energy based on the algorithm @@ -777,11 +718,11 @@ energy as in equation :eq:`maxwell-restriction`. This algorithm has a much higher rejection efficiency than the standard technique, i.e. rule C45 in the `Monte Carlo Sampler`_. -ACE Law 11 - Energy-Dependent Watt Spectrum -+++++++++++++++++++++++++++++++++++++++++++ +Energy-Dependent Watt Spectrum +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -The probability distribution for a Watt fission spectrum can be written in the -form +The probability distribution for a [Watt]_ fission spectrum can be written in +the form .. math:: :label: watt-spectrum @@ -805,29 +746,37 @@ where :math:`\xi` is a random number sampled on the interval :math:`[0,1)`. The outgoing energy is only accepted according to a specified restriction energy :math:`U` as defined in equation :eq:`maxwell-restriction`. -This algorithm can be found in Forrest Brown's lectures_ on Monte Carlo methods -and is an unpublished sampling scheme based on the original Watt spectrum -derivation [Watt]_. +A derivation of the algorithm described here can be found in a paper by Romano_. -ACE Law 44 - Kalbach-Mann Correlated Scattering -+++++++++++++++++++++++++++++++++++++++++++++++ +Product Angle-Energy Distributions +---------------------------------- -This law is very similar to ACE Law 4 except now the outgoing angle of the -neutron is correlated to the outgoing energy and is not sampled from a separate -distribution. For each incident neutron energy :math:`E_i` tabulated, there is -an array of precompound factors :math:`R_{i,j}` and angular distribution slopes -:math:`A_{i,j}` corresponding to each outgoing energy bin :math:`j` in addition -to the outgoing energies and distribution functions as in ACE Law 4. +If the secondary distribution for a product was given in file 6 in ENDF, the +angle and energy are correlated with one another and cannot be sampled +separately. Several representations exist in ENDF/ACE for correlated +angle-energy distributions. + +Kalbach-Mann Correlated Scattering +++++++++++++++++++++++++++++++++++ + +This law is very similar to the uncorrelated continuous tabular energy +distribution except now the outgoing angle of the neutron is correlated to the +outgoing energy and is not sampled from a separate distribution. For each +incident neutron energy :math:`E_i` tabulated, there is an array of precompound +factors :math:`R_{i,j}` and angular distribution slopes :math:`A_{i,j}` +corresponding to each outgoing energy bin :math:`j` in addition to the outgoing +energies and distribution functions as in :ref:`continuous-tabular`. The calculation of the outgoing energy of the neutron proceeds exactly the same -as in the algorithm described in :ref:`ace-law-4`. In that algorithm, we found -an interpolation factor :math:`f`, statistically sampled an incoming energy bin -:math:`\ell`, and sampled an outgoing energy bin :math:`j` based on the -tabulated cumulative distribution function. Once the outgoing energy has been -determined with equation :eq:`ace-law-4-energy`, we then need to calculate the -outgoing angle based on the tabulated Kalbach-Mann parameters. These parameters -themselves are subject to either histogram or linear-linear interpolation on the -outgoing energy grid. For histogram interpolation, the parameters are +as in the algorithm described in :ref:`continuous-tabular`. In that algorithm, +we found an interpolation factor :math:`f`, statistically sampled an incoming +energy bin :math:`\ell`, and sampled an outgoing energy bin :math:`j` based on +the tabulated cumulative distribution function. Once the outgoing energy has +been determined with equation :eq:`continuous-eout`, we then need to calculate +the outgoing angle based on the tabulated Kalbach-Mann parameters. These +parameters themselves are subject to either histogram or linear-linear +interpolation on the outgoing energy grid. For histogram interpolation, the +parameters are .. math:: :label: KM-parameters-histogram @@ -873,52 +822,55 @@ outgoing angle is \mu = \frac{1}{A} \ln \left ( \xi_4 e^A + (1 - \xi_4) e^{-A} \right ). -.. _ace-law-61: +.. _correlated-energy-angle: -ACE Law 61 - Correlated Energy and Angle Distribution -+++++++++++++++++++++++++++++++++++++++++++++++++++++ +Correlated Energy and Angle Distribution +++++++++++++++++++++++++++++++++++++++++ -This law is very similar to ACE Law 44 in the sense that the outgoing angle of -the neutron is correlated to the outgoing energy and is not sampled from a -separate distribution. In this case though, rather than being determined from an -analytical distribution function, the cosine of the scattering angle is -determined from a tabulated distribution. For each incident energy :math:`i` and -outgoing energy :math:`j`, there is a tabulated angular distribution. +This distribution is very similar to a Kalbach-Mann distribution in the sense +that the outgoing angle of the neutron is correlated to the outgoing energy and +is not sampled from a separate distribution. In this case though, rather than +being determined from an analytical distribution function, the cosine of the +scattering angle is determined from a tabulated distribution. For each incident +energy :math:`i` and outgoing energy :math:`j`, there is a tabulated angular +distribution. The calculation of the outgoing energy of the neutron proceeds exactly the same -as in the algorithm described in :ref:`ace-law-4`. In that algorithm, we found -an interpolation factor :math:`f`, statistically sampled an incoming energy bin -:math:`\ell`, and sampled an outgoing energy bin :math:`j` based on the -tabulated cumulative distribution function. Once the outgoing energy has been -determined with equation :eq:`ace-law-4-energy`, we then need to decide which -angular distribution to use. If histogram interpolation was used on the outgoing -energy bins, then we use the angular distribution corresponding to incoming -energy bin :math:`\ell` and outgoing energy bin :math:`j`. If linear-linear -interpolation was used on the outgoing energy bins, then we use the whichever -angular distribution was closer to the sampled value of the cumulative -distribution function for the outgoing energy. The actual algorithm used to -sample the chosen tabular angular distribution has been previously described in -:ref:`angle-tabular`. +as in the algorithm described in :ref:`continuous-tabular`. In that algorithm, +we found an interpolation factor :math:`f`, statistically sampled an incoming +energy bin :math:`\ell`, and sampled an outgoing energy bin :math:`j` based on +the tabulated cumulative distribution function. Once the outgoing energy has +been determined with equation :eq:`continuous-eout`, we then need to decide +which angular distribution to use. If histogram interpolation was used on the +outgoing energy bins, then we use the angular distribution corresponding to +incoming energy bin :math:`\ell` and outgoing energy bin :math:`j`. If +linear-linear interpolation was used on the outgoing energy bins, then we use +the whichever angular distribution was closer to the sampled value of the +cumulative distribution function for the outgoing energy. The actual algorithm +used to sample the chosen tabular angular distribution has been previously +described in :ref:`angle-tabular`. -ACE Law 66 - N-Body Phase Space Distribution -++++++++++++++++++++++++++++++++++++++++++++ +N-Body Phase Space Distribution ++++++++++++++++++++++++++++++++ Reactions in which there are more than two products of similar masses are sometimes best treated by using what's known as an N-body phase distribution. This distribution has the following probability density function -for outgoing energy of the :math:`i`-th particle in the center-of-mass system: +for outgoing energy and angle of the :math:`i`-th particle in the center-of-mass +system: .. math:: :label: n-body-pdf - p_i(E') dE' = C_n \sqrt{E'} (E_i^{max} - E')^{(3n/2) - 4} dE' + p_i(\mu, E') dE' d\mu = C_n \sqrt{E'} (E_i^{max} - E')^{(3n/2) - 4} dE' d\mu where :math:`n` is the number of outgoing particles, :math:`C_n` is a normalization constant, :math:`E_i^{max}` is the maximum center-of-mass energy -for particle :math:`i`, and :math:`E'` is the outgoing energy. The algorithm for -sampling the outgoing energy is based on algorithms R28, C45, and C64 in the -`Monte Carlo Sampler`_. First we calculate the maximum energy in the -center-of-mass using the following equation: +for particle :math:`i`, and :math:`E'` is the outgoing energy. We see in +equation :eq:`n-body-pdf` that the angle is simply isotropic in the +center-of-mass system. The algorithm for sampling the outgoing energy is based +on algorithms R28, C45, and C64 in the `Monte Carlo Sampler`_. First we +calculate the maximum energy in the center-of-mass using the following equation: .. math:: :label: n-body-emax @@ -961,7 +913,7 @@ distribution. First, the documentation (and code) for MCNP5-1.60 has a mistake in the algorithm for :math:`n = 4`. That being said, there are no existing nuclear data evaluations which use an N-body phase space distribution with :math:`n = 4`, so the error would not affect any calculations. In the -ENDF/B-VII.0 nuclear data evaluation, only one reaction uses an N-body phase +ENDF/B-VII.1 nuclear data evaluation, only one reaction uses an N-body phase space distribution at all, the :math:`(n,2n)` reaction with H-2. .. _transform-coordinates: @@ -1527,16 +1479,16 @@ accordingly. Continuous Outgoing Energies ++++++++++++++++++++++++++++ -If the thermal data was processed with :math:`iwt=2` in NJOY, then the -outgoing energy spectra is represented by a continuous outgoing energy spectra -in tabular form with linear-linear interpolation. The sampling of the outgoing -energy portion of this format is very similar to :ref:`ACE Law 61`, -but the sampling of the correlated angle is performed as it was in the other -two representations discussed in this sub-section. In the Law 61 algorithm, -we found an interpolation factor :math:`f`, statistically sampled an incoming +If the thermal data was processed with :math:`iwt=2` in NJOY, then the outgoing +energy spectra is represented by a continuous outgoing energy spectra in tabular +form with linear-linear interpolation. The sampling of the outgoing energy +portion of this format is very similar to :ref:`correlated-energy-angle`, but +the sampling of the correlated angle is performed as it was in the other two +representations discussed in this sub-section. In the Law 61 algorithm, we +found an interpolation factor :math:`f`, statistically sampled an incoming energy bin :math:`\ell`, and sampled an outgoing energy bin :math:`j` based on the tabulated cumulative distribution function. Once the outgoing energy has -been determined with equation :eq:`ace-law-4-energy`, we then need to decide +been determined with equation :eq:`continuous-eout`, we then need to decide which angular distribution data to use. Like the linear-linear interpolation case in Law 61, the angular distribution closest to the sampled value of the cumulative distribution function for the outgoing energy is utilized. The @@ -1723,6 +1675,8 @@ another. .. _MC21: http://www.osti.gov/bridge/servlets/purl/903083-HT5p1o/903083.pdf +.. _Romano: http://dx.doi.org/10.1016/j.cpc.2014.11.001 + .. _Sutton and Brown: http://www.osti.gov/bridge/product.biblio.jsp?osti_id=307911 .. _lectures: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-05-4983.pdf diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst index a17af8ac6..854c2bc4c 100644 --- a/docs/source/pythonapi/index.rst +++ b/docs/source/pythonapi/index.rst @@ -328,8 +328,11 @@ Core Classes :nosignatures: :template: myclass.rst + openmc.data.IncidentNeutron + openmc.data.Reaction openmc.data.Product openmc.data.Tabulated1D + openmc.data.ThermalScattering openmc.data.CoherentElastic Angle-Energy Distributions @@ -370,11 +373,6 @@ Classes openmc.data.ace.Library openmc.data.ace.Table - openmc.data.ace.NeutronTable - openmc.data.ace.SabTable - openmc.data.ace.PhotoatomicTable - openmc.data.ace.PhotonuclearTable - openmc.data.ace.Reaction Functions +++++++++ diff --git a/openmc/data/__init__.py b/openmc/data/__init__.py index 32de3b598..c6a1baf9d 100644 --- a/openmc/data/__init__.py +++ b/openmc/data/__init__.py @@ -1,4 +1,6 @@ from .data import * +from .neutron import * +from .reaction import * from .ace import * from .angle_distribution import * from .container import * @@ -11,3 +13,4 @@ from .kalbach_mann import * from .nbody import * from .thermal import * from .urr import * +from .library import * diff --git a/openmc/data/ace.py b/openmc/data/ace.py index 358d98789..2a39e208a 100644 --- a/openmc/data/ace.py +++ b/openmc/data/ace.py @@ -21,26 +21,9 @@ from os import SEEK_CUR import struct import sys from warnings import warn -from collections import OrderedDict -from copy import deepcopy import numpy as np -from numpy.polynomial import Polynomial -import h5py -from . import atomic_number, atomic_symbol, reaction_name -from .container import Tabulated1D, interpolation_scheme -from .angle_distribution import AngleDistribution -from .energy_distribution import * -from .product import Product -from .angle_energy import AngleEnergy -from .kalbach_mann import KalbachMann -from .uncorrelated import UncorrelatedAngleEnergy -from .correlated import CorrelatedAngleEnergy -from .nbody import NBodyPhaseSpace -from .thermal import CoherentElastic -from .urr import ProbabilityTables -from openmc.stats import Tabular, Discrete, Uniform, Mixture if sys.version_info[0] >= 3: basestring = str @@ -120,47 +103,6 @@ def ascii_to_binary(ascii_file, binary_file): binary.close() -def _get_tabulated_1d(array, idx=0): - """Create a Tabulated1D object from array. - - Parameters - ---------- - array : numpy.ndarray - Array is formed as a 1 dimensional array as follows: [number of regions, - final pair for each region, interpolation parameters, number of pairs, - x-values, y-values] - idx : int, optional - Offset to read from in array (default of zero) - - Returns - ------- - openmc.data.Tabulated1D - Tabulated data object - - """ - - # Get number of regions and pairs - n_regions = int(array[idx]) - n_pairs = int(array[idx + 1 + 2*n_regions]) - - # Get interpolation information - idx += 1 - if n_regions > 0: - nbt = np.asarray(array[idx:idx + n_regions], dtype=int) - interp = np.asarray(array[idx + n_regions:idx + 2*n_regions], dtype=int) - else: - # NR=0 regions implies linear-linear interpolation by default - nbt = np.array([n_pairs]) - interp = np.array([2]) - - # Get (x,y) pairs - idx += 2*n_regions + 1 - x = array[idx:idx + n_pairs] - y = array[idx + n_pairs:idx + 2*n_pairs] - - return Tabulated1D(x, y, nbt, interp) - - def get_table(filename, name=None): """Read a single table from an ACE file @@ -181,30 +123,14 @@ def get_table(filename, name=None): lib = Library(filename) if name is None: - return list(lib.tables.values())[0] + return lib.tables[0] else: - return lib.tables[name] - - -def get_all_tables(filename): - """Read all tables from an ACE file - - Parameters - ---------- - filename : str - Path of the ACE library to load table from - name : str, optional - Name of table to load, e.g. '92235.71c' - - Returns - ------- - list of openmc.data.ace.Table - ACE tables read from the file - - """ - - lib = Library(filename) - return list(lib.tables.values()) + for table in lib.tables: + if table.name == name: + return table + else: + raise ValueError('Could not find ACE table with name: {}' + .format(name)) class Library(object): @@ -224,10 +150,8 @@ class Library(object): Attributes ---------- - tables : dict - Dictionary whose keys are the names of the ACE tables and whose values - are the instances of subclasses of :class:`Table` - (e.g. :class:`NeutronTable`) + tables : list + List of :class:`Table` instances """ @@ -237,7 +161,7 @@ class Library(object): if table_names is not None: table_names = set(table_names) - self.tables = {} + self.tables = [] # Determine whether file is ASCII or binary try: @@ -291,10 +215,11 @@ class Library(object): # material name, atomic_weight_ratio, temperature, date, comment, mat = \ struct.unpack(str('=10sdd10s70s10s'), fh.read(116)) - name = name.strip() + name = name.decode().strip() # Read ZAID/awr combinations - izaw_pairs = struct.unpack(str('=' + 16*'id'), fh.read(192)) + data = struct.unpack(str('=' + 16*'id'), fh.read(192)) + pairs = list(zip(data[::2], data[1::2])) # Read NXS nxs = list(struct.unpack(str('=16i'), fh.read(64))) @@ -303,30 +228,14 @@ class Library(object): length = nxs[0] n_records = (length + entries - 1)//entries - # name is bytes, make it a string - name = name.decode() # verify that we are supposed to read this table in if (table_names is not None) and (name not in table_names): fh.seek(start_position + recl_length*(n_records + 1)) continue - # ensure we have a valid table type - if len(name) == 0 or name[-1] not in table_types: - warn("Unsupported table: " + name, RuntimeWarning) - fh.seek(start_position + recl_length*(n_records + 1)) - continue - - # get the table - table = table_types[name[-1]](name, atomic_weight_ratio, temperature) - if verbose: temperature_in_K = round(temperature * 1e6 / 8.617342e-5) print("Loading nuclide {0} at {1} K".format(name, temperature_in_K)) - self.tables[name] = table - - # If table is S(a,b), add zaids - zaids = np.array(izaw_pairs[::2]) - table.zaids = zaids[np.nonzero(zaids)] # Read JXS jxs = list(struct.unpack(str('=32i'), fh.read(128))) @@ -336,20 +245,22 @@ class Library(object): xss = list(struct.unpack(str('={0}d'.format(length)), fh.read(length*8))) - # Insert empty object at beginning of NXS, JXS, and XSS arrays so - # that the indexing will be the same as Fortran. This makes it - # easier to follow the ACE format specification. + # Insert zeros at beginning of NXS, JXS, and XSS arrays so that the + # indexing will be the same as Fortran. This makes it easier to + # follow the ACE format specification. nxs.insert(0, 0) - table._nxs = np.array(nxs, dtype=int) + nxs = np.array(nxs, dtype=int) jxs.insert(0, 0) - table._jxs = np.array(jxs, dtype=int) + jxs = np.array(jxs, dtype=int) xss.insert(0, 0.0) - table._xss = np.array(xss) + xss = np.array(xss) - # Read all data blocks - table._read_all() + # Create ACE table with data read in + table = Table(name, atomic_weight_ratio, temperature, pairs, + nxs, jxs, xss) + self.tables.append(table) # Advance to next record fh.seek(start_position + recl_length*(n_records + 1)) @@ -397,7 +308,9 @@ class Library(object): atomic_weight_ratio = float(words[1]) temperature = float(words[2]) - izaw_pairs = (' '.join(lines[2:6])).split() + datastr = ' '.join(lines[2:6]).split() + pairs = list(zip(map(int, datastr[::2]), + map(float, datastr[1::2]))) datastr = '0 ' + ' '.join(lines[6:8]) nxs = np.fromstring(datastr, sep=' ', dtype=int) @@ -417,14 +330,6 @@ class Library(object): lines = [fh.readline() for i in range(13)] continue - # ensure we have a valid table type - if len(name) == 0 or name[-1] not in table_types: - warn("Unsupported table: " + name, RuntimeWarning) - fh.seek(n_bytes, SEEK_CUR) - fh.readline() - lines = [fh.readline() for i in range(13)] - continue - # read and fix over-shoot lines += fh.readlines(n_bytes) if 12 + n_lines < len(lines): @@ -432,41 +337,32 @@ class Library(object): lines = lines[:12+n_lines] fh.seek(-goback, SEEK_CUR) - # get the table - table = table_types[name[-1]](name, atomic_weight_ratio, temperature) - if verbose: temperature_in_K = round(temperature * 1e6 / 8.617342e-5) print("Loading nuclide {0} at {1} K".format(name, temperature_in_K)) - self.tables[name] = table # Read comment - table.comment = lines[1].strip() - - # If table is S(a,b), add zaids - if isinstance(table, SabTable): - zaids = np.fromiter(map(int, izaw_pairs[::2]), int) - table.zaids = zaids[np.nonzero(zaids)] - - # Add NXS, JXS, and XSS arrays to table Insert empty object at - # beginning of NXS, JXS, and XSS arrays so that the indexing will be - # the same as Fortran. This makes it easier to follow the ACE format - # specification. - table._nxs = nxs + comment = lines[1].strip() + # Insert zeros at beginning of NXS, JXS, and XSS arrays so that the + # indexing will be the same as Fortran. This makes it easier to + # follow the ACE format specification. datastr = '0 ' + ' '.join(lines[8:12]) - table._jxs = np.fromstring(datastr, dtype=int, sep=' ') + jxs = np.fromstring(datastr, dtype=int, sep=' ') datastr = '0.0 ' + ''.join(lines[12:12+n_lines]) - table._xss = np.fromstring(datastr, sep=' ') + xss = np.fromstring(datastr, sep=' ') + + table = Table(name, atomic_weight_ratio, temperature, pairs, + nxs, jxs, xss) + self.tables.append(table) # Read all data blocks - table._read_all() lines = [fh.readline() for i in range(13)] class Table(object): - """Abstract superclass of all other classes for cross section tables. + """ACE cross section table Parameters ---------- @@ -475,1387 +371,28 @@ class Table(object): atomic_weight_ratio : float Atomic mass ratio of the target nuclide. temperature : float - Temperature of the target nuclide in eV. - - Attributes - ---------- - name : str - ZAID identifier of the table, e.g. '92235.70c'. - atomic_weight_ratio : float - Atomic mass ratio of the target nuclide. - temperature : float - Temperature of the target nuclide in eV. + Temperature of the target nuclide in MeV. + pairs : list of tuple + 16 pairs of ZAIDs and atomic weight ratios. Used for thermal scattering + tables to indicate what isotopes scattering is applied to. + nxs : numpy.ndarray + Array that defines various lengths with in the table + jxs : numpy.ndarray + Array that gives locations in the ``xss`` array for various blocks of + data + xss : numpy.ndarray + Raw data for the ACE table """ - - def __init__(self, name, atomic_weight_ratio, temperature): + def __init__(self, name, atomic_weight_ratio, temperature, pairs, + nxs, jxs, xss): self.name = name self.atomic_weight_ratio = atomic_weight_ratio self.temperature = temperature - - def _read_all(self): - raise NotImplementedError - - def _get_continuous_tabular(self, idx, ldis): - """Get continuous tabular energy distribution (ACE law 4) starting at specified - index in the XSS array. - - Parameters - ---------- - idx : int - Index in XSS array of the start of the energy distribution data - (LDIS + LOCC - 1) - ldis : int - Index in XSS array of the start of the energy distribution block - (e.g. JXS[11]) - - Returns - ------- - openmc.data.energy_distribution.ContinuousTabular - Continuous tabular energy distribution - - """ - - # Read number of interpolation regions and incoming energies - n_regions = int(self._xss[idx]) - n_energy_in = int(self._xss[idx + 1 + 2*n_regions]) - - # Get interpolation information - idx += 1 - if n_regions > 0: - breakpoints = np.asarray(self._xss[idx:idx + n_regions], dtype=int) - interpolation = np.asarray(self._xss[idx + n_regions: - idx + 2*n_regions], dtype=int) - else: - breakpoints = np.array([n_energy_in]) - interpolation = np.array([2]) - - # Incoming energies at which distributions exist - idx += 2 * n_regions + 1 - energy = self._xss[idx:idx + n_energy_in] - - # Location of distributions - idx += n_energy_in - loc_dist = np.asarray(self._xss[idx:idx + n_energy_in], dtype=int) - - # Initialize variables - energy_out = [] - - # Read each outgoing energy distribution - for i in range(n_energy_in): - idx = ldis + loc_dist[i] - 1 - - # intt = interpolation scheme (1=hist, 2=lin-lin) - INTTp = int(self._xss[idx]) - intt = INTTp % 10 - n_discrete_lines = (INTTp - intt)//10 - if intt not in (1, 2): - warn("Interpolation scheme for continuous tabular distribution " - "is not histogram or linear-linear.") - intt = 2 - - n_energy_out = int(self._xss[idx + 1]) - data = self._xss[idx + 2:idx + 2 + 3*n_energy_out] - data.shape = (3, n_energy_out) - - # Create continuous distribution - eout_continuous = Tabular(data[0][n_discrete_lines:], - data[1][n_discrete_lines:], - interpolation_scheme[intt]) - eout_continuous.c = data[2][n_discrete_lines:] - - # If discrete lines are present, create a mixture distribution - if n_discrete_lines > 0: - eout_discrete = Discrete(data[0][:n_discrete_lines], - data[1][:n_discrete_lines]) - eout_discrete.c = data[2][:n_discrete_lines] - if n_discrete_lines == n_energy_out: - eout_i = eout_discrete - else: - p_discrete = min(sum(eout_discrete.p), 1.0) - eout_i = Mixture([p_discrete, 1. - p_discrete], - [eout_discrete, eout_continuous]) - else: - eout_i = eout_continuous - - energy_out.append(eout_i) - - return ContinuousTabular(breakpoints, interpolation, energy, - energy_out) - - def _get_general_evaporation(self, idx): - # Read nuclear temperature as Tabulated1D - theta = _get_tabulated_1d(array, idx) - - # X-function - nr = int(array[idx]) - ne = int(array[idx + 1 + 2*nr]) - idx += 2 + 2*nr + 2*ne - net = int(array[idx]) - x = array[idx + 1:idx + 1 + net] - - raise NotImplementedError("Where'd you get this ACE file from?") - - def _get_maxwell_energy(self, idx): - # Read nuclear temperature as Tabulated1D - theta = _get_tabulated_1d(self._xss, idx) - - # Restriction energy - nr = int(self._xss[idx]) - ne = int(self._xss[idx + 1 + 2*nr]) - u = self._xss[idx + 2 + 2*nr + 2*ne] - - return MaxwellEnergy(theta, u) - - def _get_evaporation(self, idx): - # Read nuclear temperature as Tabulated1D - theta = _get_tabulated_1d(self._xss, idx) - - # Restriction energy - nr = int(self._xss[idx]) - ne = int(self._xss[idx + 1 + 2*nr]) - u = self._xss[idx + 2 + 2*nr + 2*ne] - - return Evaporation(theta, u) - - def _get_watt_energy(self, idx): - # Energy-dependent a parameter - a = _get_tabulated_1d(self._xss, idx) - - # Advance index - nr = int(self._xss[idx]) - ne = int(self._xss[idx + 1 + 2*nr]) - idx += 2 + 2*nr + 2*ne - - # Energy-dependent b parameter - b = _get_tabulated_1d(self._xss, idx) - - # Advance index - nr = int(self._xss[idx]) - ne = int(self._xss[idx + 1 + 2*nr]) - idx += 2 + 2*nr + 2*ne - - # Restriction energy - u = self._xss[idx] - - return WattEnergy(a, b, u) - - def _get_kalbach_mann(self, idx, ldis): - # Read number of interpolation regions and incoming energies - n_regions = int(self._xss[idx]) - n_energy_in = int(self._xss[idx + 1 + 2*n_regions]) - - # Get interpolation information - idx += 1 - if n_regions > 0: - breakpoints = np.asarray(self._xss[idx:idx + n_regions], dtype=int) - interpolation = np.asarray(self._xss[idx + n_regions: - idx + 2*n_regions], dtype=int) - else: - breakpoints = np.array([n_energy_in]) - interpolation = np.array([2]) - - # Incoming energies at which distributions exist - idx += 2 * n_regions + 1 - energy = self._xss[idx:idx + n_energy_in] - - # Location of distributions - idx += n_energy_in - loc_dist = np.asarray(self._xss[idx:idx + n_energy_in], dtype=int) - - # Initialize variables - energy_out = [] - km_r = [] - km_a = [] - - # Read each outgoing energy distribution - for i in range(n_energy_in): - idx = ldis + loc_dist[i] - 1 - - # intt = interpolation scheme (1=hist, 2=lin-lin) - INTTp = int(self._xss[idx]) - intt = INTTp % 10 - n_discrete_lines = (INTTp - intt)//10 - if intt not in (1, 2): - warn("Interpolation scheme for continuous tabular distribution " - "is not histogram or linear-linear.") - intt = 2 - - n_energy_out = int(self._xss[idx + 1]) - data = self._xss[idx + 2:idx + 2 + 5*n_energy_out] - data.shape = (5, n_energy_out) - - # Create continuous distribution - eout_continuous = Tabular(data[0][n_discrete_lines:], - data[1][n_discrete_lines:], - interpolation_scheme[intt]) - eout_continuous.c = data[2][n_discrete_lines:] - - # If discrete lines are present, create a mixture distribution - if n_discrete_lines > 0: - eout_discrete = Discrete(data[0][:n_discrete_lines], - data[1][:n_discrete_lines]) - eout_discrete.c = data[2][:n_discrete_lines] - if n_discrete_lines == n_energy_out: - eout_i = eout_discrete - else: - p_discrete = min(sum(eout_discrete.p), 1.0) - eout_i = Mixture([p_discrete, 1. - p_discrete], - [eout_discrete, eout_continuous]) - else: - eout_i = eout_continuous - - energy_out.append(eout_i) - km_r.append(Tabulated1D(data[0], data[3])) - km_a.append(Tabulated1D(data[0], data[4])) - - return KalbachMann(breakpoints, interpolation, energy, energy_out, - km_r, km_a) - - def _get_correlated(self, idx, ldis): - # Read number of interpolation regions and incoming energies - n_regions = int(self._xss[idx]) - n_energy_in = int(self._xss[idx + 1 + 2*n_regions]) - - # Get interpolation information - idx += 1 - if n_regions > 0: - breakpoints = np.asarray(self._xss[idx:idx + n_regions], dtype=int) - interpolation = np.asarray(self._xss[idx + n_regions: - idx + 2*n_regions], dtype=int) - else: - breakpoints = np.array([n_energy_in]) - interpolation = np.array([2]) - - # Incoming energies at which distributions exist - idx += 2 * n_regions + 1 - energy = self._xss[idx:idx + n_energy_in] - - # Location of distributions - idx += n_energy_in - loc_dist = np.asarray(self._xss[idx:idx + n_energy_in], dtype=int) - - # Initialize list of distributions - energy_out = [] - mu = [] - - # Read each outgoing energy distribution - for i in range(n_energy_in): - idx = ldis + loc_dist[i] - 1 - - # intt = interpolation scheme (1=hist, 2=lin-lin) - INTTp = int(self._xss[idx]) - intt = INTTp % 10 - n_discrete_lines = (INTTp - intt)//10 - if intt not in (1, 2): - warn("Interpolation scheme for continuous tabular distribution " - "is not histogram or linear-linear.") - intt = 2 - - # Secondary energy distribution - n_energy_out = int(self._xss[idx + 1]) - data = self._xss[idx + 2:idx + 2 + 4*n_energy_out] - data.shape = (4, n_energy_out) - - # Create continuous distribution - eout_continuous = Tabular(data[0][n_discrete_lines:], - data[1][n_discrete_lines:], - interpolation_scheme[intt], - ignore_negative=True) - eout_continuous.c = data[2][n_discrete_lines:] - - # If discrete lines are present, create a mixture distribution - if n_discrete_lines > 0: - eout_discrete = Discrete(data[0][:n_discrete_lines], - data[1][:n_discrete_lines]) - eout_discrete.c = data[2][:n_discrete_lines] - if n_discrete_lines == n_energy_out: - eout_i = eout_discrete - else: - p_discrete = min(sum(eout_discrete.p), 1.0) - eout_i = Mixture([p_discrete, 1. - p_discrete], - [eout_discrete, eout_continuous]) - else: - eout_i = eout_continuous - - energy_out.append(eout_i) - - lc = np.asarray(data[3], dtype=int) - - # Secondary angular distributions - mu_i = [] - for j in range(n_energy_out): - if lc[j] > 0: - idx = ldis + abs(lc[j]) - 1 - - intt = int(self._xss[idx]) - n_cosine = int(self._xss[idx + 1]) - data = self._xss[idx + 2:idx + 2 + 3*n_cosine] - data.shape = (3, n_cosine) - - mu_ij = Tabular(data[0], data[1], interpolation_scheme[intt]) - mu_ij.c = data[2] - else: - # Isotropic distribution - mu_ij = Uniform(-1., 1.) - - mu_i.append(mu_ij) - - # Add cosine distributions for this incoming energy to list - mu.append(mu_i) - - return CorrelatedAngleEnergy(breakpoints, interpolation, energy, - energy_out, mu) - - def _get_energy_distribution(self, location_dist, location_start, rx=None): - """Returns an EnergyDistribution object from data read in starting at - location_start. - - Parameters - ---------- - location_dist : int - Index in the XSS array corresponding to the start of a block, - e.g. JXS(11) for the the DLW block. - location_start : int - Index in the XSS array corresponding to the start of an energy - distribution array - rx : Reaction - Reaction this energy distribution will be associated with - - Returns - ------- - distribution : openmc.data.AngleEnergy - Secondary angle-energy distribution - - """ - - # Set starting index for energy distribution - idx = location_dist + location_start - 1 - - law = int(self._xss[idx + 1]) - location_data = int(self._xss[idx + 2]) - - # Position index for reading law data - idx = location_dist + location_data - 1 - - # Parse energy distribution data - if law == 2: - primary_flag = int(self._xss[idx]) - energy = self._xss[idx + 1] - distribution = UncorrelatedAngleEnergy() - distribution.energy = DiscretePhoton(primary_flag, energy, - self.atomic_weight_ratio) - elif law in (3, 33): - threshold, mass_ratio = self._xss[idx:idx + 2] - distribution = UncorrelatedAngleEnergy() - distribution.energy = LevelInelastic(threshold, mass_ratio) - elif law == 4: - distribution = UncorrelatedAngleEnergy() - distribution.energy = self._get_continuous_tabular(idx, location_dist) - elif law == 5: - distribution = UncorrelatedAngleEnergy() - distribution.energy = self._get_general_evaporation(idx) - elif law == 7: - distribution = UncorrelatedAngleEnergy() - distribution.energy = self._get_maxwell_energy(idx) - elif law == 9: - distribution = UncorrelatedAngleEnergy() - distribution.energy = self._get_evaporation(idx) - elif law == 11: - distribution = UncorrelatedAngleEnergy() - distribution.energy = self._get_watt_energy(idx) - elif law == 44: - distribution = self._get_kalbach_mann(idx, location_dist) - elif law == 61: - distribution = self._get_correlated(idx, location_dist) - elif law == 66: - n_particles = int(self._xss[idx]) - total_mass = self._xss[idx + 1] - distribution = NBodyPhaseSpace(total_mass, n_particles, - self.atomic_weight_ratio, rx.Q_value) - else: - raise IOError("Unsupported ACE secondary energy " - "distribution law {0}".format(law)) - - return distribution - - -class NeutronTable(Table): - """A NeutronTable object contains continuous-energy neutron interaction data - read from an ACE-formatted table. These objects are not normally - instantiated by the user but rather created when reading data using a - Library object and stored within the :attr:`Library.tables` attribute. - - Parameters - ---------- - name : str - ZAID identifier of the table, e.g. '92235.70c'. - atomic_weight_ratio : float - Atomic mass ratio of the target nuclide. - temperature : float - Temperature of the target nuclide in eV. - - Attributes - ---------- - absorption_xs : numpy.ndarray - The microscopic absorption cross section for each value on the energy - grid. - atomic_weight_ratio : float - Atomic weight ratio of the target nuclide. - energy : numpy.ndarray - The energy values (MeV) at which reaction cross-sections are tabulated. - heating_number : numpy.ndarray - The total heating number for each value on the energy grid in MeV-b. - name : str - ZAID identifier of the table, e.g. 92235.70c. - reactions : collections.OrderedDict - 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. - temperature : float - Temperature of the target nuclide in eV. - total_xs : numpy.ndarray - The microscopic total cross section for each value on the energy grid in b. - urr : None or openmc.data.ProbabilityTables - Unresolved resonance region probability tables - - """ - - def __init__(self, name, atomic_weight_ratio, temperature): - super(NeutronTable, self).__init__(name, atomic_weight_ratio, temperature) - self.absorption_xs = None - self.energy = None - self.heating_number = None - self.total_xs = None - self.reactions = OrderedDict() - self.urr = None + self.pairs = pairs + self.nxs = nxs + self.jxs = jxs + self.xss = xss def __repr__(self): - if hasattr(self, 'name'): - return "".format(self.name) - else: - return "" - - def __iter__(self): - return iter(self.reactions.values()) - - def _read_all(self): - self._read_cross_sections() - self._read_nu() - self._read_secondaries() - self._read_photon_production_data() - self._read_unr() - - def _read_cross_sections(self): - """Read reaction cross sections and other data. - - Reads and parses the ESZ, MTR, LQR, TRY, LSIG, and SIG blocks. These - blocks contain the energy grid, all reaction cross sections, the total - cross section, average heating numbers, and a list of reactions with - their Q-values and multiplicites. - """ - - # Determine number of energies on nuclide grid and number of reactions - # excluding elastic scattering - n_energies = self._nxs[3] - n_reactions = self._nxs[4] - - # Read energy grid and total, absorption, elastic scattering, and - # heating cross sections -- note that this appear separate from the rest - # of the reaction cross sections - arr = self._xss[self._jxs[1]:self._jxs[1] + 5*n_energies] - arr.shape = (5, n_energies) - self.energy, self.total_xs, self.absorption_xs, \ - elastic_xs, self.heating_number = arr - - # Create elastic scattering reaction - elastic_scatter = Reaction(2, self) - elastic_scatter.products.append(Product('neutron')) - elastic_scatter.xs = Tabulated1D(self.energy, elastic_xs) - self.reactions[2] = elastic_scatter - - # Create all other reactions with MT values - mts = np.asarray(self._xss[self._jxs[3]:self._jxs[3] + n_reactions], dtype=int) - qvalues = np.asarray(self._xss[self._jxs[4]:self._jxs[4] + - n_reactions], dtype=float) - tys = np.asarray(self._xss[self._jxs[5]:self._jxs[5] + n_reactions], dtype=int) - - # Create all reactions other than elastic scatter - reactions = [(mt, Reaction(mt, self)) for mt in mts] - self.reactions.update(reactions) - - # Loop over all reactions other than elastic scattering - for i, rx in enumerate(list(self.reactions.values())[1:]): - # Copy Q values determine if scattering should be treated in the - # center-of-mass or lab system - rx.Q_value = qvalues[i] - rx.center_of_mass = (tys[i] < 0) - - # For neutron-producing reactions, get yield - if rx.MT < 100: - if tys[i] != 19: - if abs(tys[i]) > 100: - # Energy-dependent neutron yield - idx = self._jxs[11] + abs(tys[i]) - 101 - yield_ = _get_tabulated_1d(self._xss, idx) - else: - yield_ = abs(tys[i]) - - neutron = Product('neutron') - neutron.yield_ = yield_ - rx.products.append(neutron) - - # Get locator for cross-section data - loc = int(self._xss[self._jxs[6] + i]) - - # Determine starting index on energy grid - rx.threshold_idx = int(self._xss[self._jxs[7] + loc - 1]) - 1 - - # Determine number of energies in reaction - n_energies = int(self._xss[self._jxs[7] + loc]) - energy = self.energy[rx.threshold_idx:rx.threshold_idx + n_energies] - - # Read reaction cross section - xs = self._xss[self._jxs[7] + loc + 1:self._jxs[7] + loc + 1 + n_energies] - rx.xs = Tabulated1D(energy, xs) - - def _read_nu(self): - """Read the NU block -- this contains information on the prompt and delayed - neutron precursor yields, decay constants, etc - - """ - # No NU block - if self._jxs[2] == 0: - return - - products = [] - derived_products = [] - - # Either prompt nu or total nu is given - if self._xss[self._jxs[2]] > 0: - whichnu = 'prompt' if self._jxs[24] > 0 else 'total' - - neutron = Product('neutron') - neutron.emission_mode = whichnu - - idx = self._jxs[2] - LNU = int(self._xss[idx]) - if LNU == 1: - # Polynomial function form of nu - NC = int(self._xss[idx+1]) - coefficients = self._xss[idx+2 : idx+2+NC] - neutron.yield_ = Polynomial(coefficients) - elif LNU == 2: - # Tabular data form of nu - neutron.yield_ = _get_tabulated_1d(self._xss, idx + 1) - - products.append(neutron) - - # Both prompt nu and total nu - elif self._xss[self._jxs[2]] < 0: - # Read prompt neutron yield - prompt_neutron = Product('neutron') - prompt_neutron.emission_mode = 'prompt' - - idx = self._jxs[2] + 1 - LNU = int(self._xss[idx]) - if LNU == 1: - # Polynomial function form of nu - NC = int(self._xss[idx+1]) - coefficients = self._xss[idx+2 : idx+2+NC] - prompt_neutron.yield_ = Polynomial(coefficients) - elif LNU == 2: - # Tabular data form of nu - prompt_neutron.yield_ = _get_tabulated_1d(self._xss, idx + 1) - - # Read total neutron yield - total_neutron = Product('neutron') - total_neutron.emission_mode = 'total' - - idx = self._jxs[2] + int(abs(self._xss[self._jxs[2]])) + 1 - LNU = int(self._xss[idx]) - - if LNU == 1: - # Polynomial function form of nu - NC = int(self._xss[idx+1]) - coefficients = self._xss[idx+2 : idx+2+NC] - total_neutron.yield_ = Polynomial(coefficients) - elif LNU == 2: - # Tabular data form of nu - total_neutron.yield_ = _get_tabulated_1d(self._xss, idx + 1) - - products.append(prompt_neutron) - derived_products.append(total_neutron) - - # Check for delayed nu data - if self._jxs[24] > 0: - yield_delayed = _get_tabulated_1d(self._xss, self._jxs[24] + 1) - - # Delayed neutron precursor distribution - idx = self._jxs[25] - n_group = self._nxs[8] - total_group_probability = 0. - for i, group in enumerate(range(n_group)): - delayed_neutron = Product('neutron') - delayed_neutron.emission_mode = 'delayed' - delayed_neutron.decay_rate = self._xss[idx] - - group_probability = _get_tabulated_1d(self._xss, idx + 1) - if np.all(group_probability.y == group_probability.y[0]): - delayed_neutron.yield_ = deepcopy(yield_delayed) - delayed_neutron.yield_.y *= group_probability.y[0] - total_group_probability += group_probability.y[0] - else: - raise NotImplementedError( - 'Delayed neutron with energy-dependent ' - 'group probability') - - # Advance position - nr = int(self._xss[idx + 1]) - ne = int(self._xss[idx + 2 + 2*nr]) - idx += 3 + 2*nr + 2*ne - - # Energy distribution for delayed fission neutrons - location_start = int(self._xss[self._jxs[26] + group]) - delayed_neutron.distribution.append( - self._get_energy_distribution(self._jxs[27], location_start)) - - products.append(delayed_neutron) - - # Renormalize delayed neutron yields to reflect fact that in ACE - # file, the sum of the group probabilities is not exactly one - for product in products[1:]: - product.yield_.y /= total_group_probability - - # Copy fission neutrons to reactions - for MT, rx in self.reactions.items(): - if MT in (18, 19, 20, 21, 38): - rx.products = deepcopy(products) - if derived_products: - rx.derived_products = deepcopy(derived_products) - - def _get_angle_distribution(self, location_dist, location_start): - # Set starting index for angle distribution - idx = location_dist + location_start - 1 - - # Number of energies at which angular distributions are tabulated - n_energies = int(self._xss[idx]) - idx += 1 - - # Incoming energy grid - energy = self._xss[idx:idx + n_energies] - idx += n_energies - - # Read locations for angular distributions - lc = np.asarray(self._xss[idx:idx + n_energies], dtype=int) - idx += n_energies - - mu = [] - for i in range(n_energies): - if lc[i] > 0: - # Equiprobable 32 bin distribution - idx = location_dist + abs(lc[i]) - 1 - cos = self._xss[idx:idx + 33] - pdf = np.zeros(33) - pdf[:32] = 1.0/(32.0*np.diff(cos)) - cdf = np.linspace(0.0, 1.0, 33) - - mu_i = Tabular(cos, pdf, 'histogram', ignore_negative=True) - mu_i.c = cdf - elif lc[i] < 0: - # Tabular angular distribution - idx = location_dist + abs(lc[i]) - 1 - intt = int(self._xss[idx]) - n_points = int(self._xss[idx + 1]) - data = self._xss[idx + 2:idx + 2 + 3*n_points] - data.shape = (3, n_points) - - mu_i = Tabular(data[0], data[1], interpolation_scheme[intt]) - mu_i.c = data[2] - else: - # Isotropic angular distribution - mu_i = Uniform(-1., 1.) - - mu.append(mu_i) - - return AngleDistribution(energy, mu) - - def _read_secondaries(self): - """Read angle/energy distributions for each reaction MT - """ - - # Number of reactions with secondary neutrons (including elastic - # scattering) - n_reactions = self._nxs[5] + 1 - - for i, rx in enumerate(list(self.reactions.values())[:n_reactions]): - if rx.MT == 18: - for p in rx.products: - if p.emission_mode == 'prompt': - neutron = p - break - else: - neutron = rx.products[0] - - if i > 0: - # Determine locator for ith energy distribution - lnw = int(self._xss[self._jxs[10] + i - 1]) - - while lnw > 0: - # Applicability of this distribution - neutron.applicability.append(_get_tabulated_1d( - self._xss, self._jxs[11] + lnw + 2)) - - # Read energy distribution data - neutron.distribution.append(self._get_energy_distribution( - self._jxs[11], lnw, rx)) - - lnw = int(self._xss[self._jxs[11] + lnw - 1]) - else: - # No energy distribution for elastic scattering - neutron.distribution.append(UncorrelatedAngleEnergy()) - - # Check if angular distribution data exist - loc = int(self._xss[self._jxs[8] + i]) - if loc == -1: - # Angular distribution data are given as part of product - # angle-energy distribution - continue - elif loc == 0: - # No angular distribution data are given for this - # reaction, isotropic scattering is asssumed - angle_dist = None - else: - angle_dist = self._get_angle_distribution(self._jxs[9], loc) - - # Apply angular distribution to each uncorrelated angle-energy - # distribution - if angle_dist is not None: - for d in neutron.distribution: - d.angle = angle_dist - - def _read_photon_production_data(self): - """Read cross sections for each photon-production reaction""" - - n_photon_reactions = self._nxs[6] - photon_mts = np.asarray(self._xss[self._jxs[13]:self._jxs[13] + - n_photon_reactions], dtype=int) - - for i, rx in enumerate(photon_mts): - # Determine corresponding reaction - mt = photon_mts[i] // 1000 - reactions = [] - if mt not in self.reactions: - # If the photon is assigned to MT=18 but the file splits fission - # into MT=19,20,21,38, assign the photon product to each of the - # individual reactions - if mt == 18: - for mt_fiss in (19, 20, 21, 38): - if mt_fiss in self.reactions: - reactions.append(self.reactions[mt_fiss]) - if not reactions: - reactions.append(Reaction(mt, self)) - else: - reactions.append(self.reactions[mt]) - - # Create photon product and assign to reactions - photon = Product('photon') - for rx in reactions: - rx.products.append(photon) - - # ================================================================== - # Read photon yield / production cross section - - loca = int(self._xss[self._jxs[14] + i]) - idx = self._jxs[15] + loca - 1 - mftype = int(self._xss[idx]) - idx += 1 - - if mftype in (12, 16): - # Yield data taken from ENDF File 12 or 6 - mtmult = int(self._xss[idx]) - assert mtmult == mt - - # Read photon yield as function of energy - photon.yield_ = _get_tabulated_1d(self._xss, idx + 1) - - elif mftype == 13: - # Cross section data from ENDF File 13 - - # Energy grid index at which data starts - threshold_idx = int(self._xss[idx]) - 1 - - # Get photon production cross section - n_energy = int(self._xss[idx + 1]) - photon._xs = self._xss[idx + 2:idx + 2 + n_energy] - - # Determine yield based on ratio of cross sections - energy = self.energy[threshold_idx:threshold_idx + n_energy] - photon.yield_ = Tabulated1D(energy, photon._xs) - - else: - raise ValueError("MFTYPE must be 12, 13, 16. Got {0}".format( - mftype)) - - # ================================================================== - # Read photon energy distribution - - location_start = int(self._xss[self._jxs[18] + i]) - - # Read energy distribution data - distribution = self._get_energy_distribution( - self._jxs[19], location_start) - assert isinstance(distribution, UncorrelatedAngleEnergy) - - # ================================================================== - # Read photon angular distribution - loc = int(self._xss[self._jxs[16] + i]) - - if loc == 0: - # No angular distribution data are given for this reaction, - # isotropic scattering is asssumed in LAB - energy = np.array([photon.yield_.x[0], photon.yield_.x[-1]]) - mu_isotropic = Uniform(-1., 1.) - distribution.angle = AngleDistribution( - energy, [mu_isotropic, mu_isotropic]) - else: - distribution.angle = self._get_angle_distribution(self._jxs[17], loc) - - # Add to list of distributions - photon.distribution.append(distribution) - - def _read_unr(self): - """Read the unresolved resonance range probability tables if present. - """ - - # Check if URR probability tables are present - idx = self._jxs[23] - if idx == 0: - return - - N = int(self._xss[idx]) # Number of incident energies - M = int(self._xss[idx+1]) # Length of probability table - interpolation = int(self._xss[idx+2]) - inelastic_flag = int(self._xss[idx+3]) - absorption_flag = int(self._xss[idx+4]) - multiply_smooth = (int(self._xss[idx+5]) == 1) - idx += 6 - - # Get energies at which tables exist - energy = self._xss[idx : idx+N] - idx += N - - # Get probability tables - table = self._xss[idx:idx+N*6*M] - table.shape = (N, 6, M) - - # Create object - self.urr = ProbabilityTables(energy, table, interpolation, inelastic_flag, - absorption_flag, multiply_smooth) - - def export_to_hdf5(self, path, element=None, mass_number=None, metastable=0): - """Export table to an HDF5 file. - - Parameters - ---------- - path : str - Path to write HDF5 file to - element : str or None - Elemental symbol, e.g. Zr. If not specified, the atomic - number/symbol are inferred from the name of the table. - mass_number : int or None - Mass number of the nuclide. For natural elements, a value of zero - should be given. If not specified, the mass number is inferred from - the name of the table. - metastable : int - Metastable level of the nuclide. Defaults to 0. - - """ - - f = h5py.File(path, 'a') - - # If element and/or mass number haven't been specified, make an educated - # guess - zaid, xs = self.name.split('.') - if element is None: - Z = int(zaid) // 1000 - element = atomic_symbol[Z] - else: - Z = atomic_number[element] - if mass_number is None: - mass_number = int(zaid) % 1000 - - # Write basic data - if metastable > 0: - name = '{}{}_m{}.{}'.format(element, mass_number, metastable, xs) - else: - name = '{}{}.{}'.format(element, mass_number, xs) - g = f.create_group(name) - g.attrs['Z'] = Z - g.attrs['A'] = mass_number - g.attrs['metastable'] = metastable - g.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio - g.attrs['temperature'] = self.temperature - g.attrs['n_reaction'] = len(self.reactions) - - # Write energy grid - g.create_dataset('energy', data=self.energy) - - # Write reaction data - for i, rx in enumerate(self.reactions.values()): - rx_group = g.create_group('reaction_{}'.format(i)) - rx.to_hdf5(rx_group) - - # Write total nu data if available - if hasattr(rx, 'derived_products') 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 is not None: - urr_group = g.create_group('urr') - self.urr.to_hdf5(urr_group) - - f.close() - - @classmethod - def from_hdf5(self, group): - """Generate continuous-energy neutron interaction data from HDF5 group - - Parameters - ---------- - group : h5py.Group - HDF5 group containing interaction data - - Returns - ------- - openmc.data.ace.NeutronTable - Continuous-energy neutron interaction data - - """ - name = group.name[1:] - atomic_weight_ratio = group.attrs['atomic_weight_ratio'] - temperature = group.attrs['temperature'] - table = NeutronTable(name, atomic_weight_ratio, temperature) - - # Read energy grid - table.energy = group['energy'].value - - # Read reaction data - n_reaction = group.attrs['n_reaction'] - - # Write reaction data - for i in range(n_reaction): - rx_group = group['reaction_{}'.format(i)] - rx = Reaction.from_hdf5(rx_group, table) - table.reactions[rx.MT] = rx - - # Read total nu data if available - if 'total_nu' in rx_group: - tgroup = rx_group['total_nu'] - rx.derived_products = [Product.from_hdf5(tgroup)] - - # Read unresolved resonance probability tables - if 'urr' in group: - urr_group = group['urr'] - table.urr = ProbabilityTables.from_hdf5(urr_group) - - return table - - -class SabTable(Table): - """A SabTable object contains thermal scattering data as represented by - an S(alpha, beta) table. - - Parameters - ---------- - name : str - ZAID identifier of the table, e.g. lwtr.10t. - atomic_weight_ratio : float - Atomic mass ratio of the target nuclide. - temperature : float - Temperature of the target nuclide in eV. - - Attributes - ---------- - atomic_weight_ratio : float - Atomic mass ratio of the target nuclide. - elastic_xs : openmc.data.Tabulated1D or openmc.data.CoherentElastic - Elastic scattering cross section derived in the coherent or incoherent - approximation - inelastic_xs : openmc.data.Tabulated1D - Inelastic scattering cross section derived in the incoherent - approximation - name : str - ZAID identifier of the table, e.g. 92235.70c. - temperature : float - Temperature of the target nuclide in eV. - - """ - - def __init__(self, name, atomic_weight_ratio, temperature): - super(SabTable, self).__init__(name, atomic_weight_ratio, temperature) - self.elastic_xs = None - self.elastic_mu_out = None - - self.inelastic_xs = None - self.inelastic_e_out = None - self.inelastic_mu_out = None - self.secondary_mode = None - - def __repr__(self): - if hasattr(self, 'name'): - return "".format(self.name) - else: - return "" - - def _read_all(self): - self._read_itie() - self._read_itce() - self._read_itxe() - self._read_itca() - - def _read_itie(self): - """Read energy-dependent inelastic scattering cross sections. - """ - idx = self._jxs[1] - n_energies = int(self._xss[idx]) - energy = self._xss[idx+1 : idx+1+n_energies] - xs = self._xss[idx+1+n_energies : idx+1+2*n_energies] - self.inelastic_xs = Tabulated1D(energy, xs) - - def _read_itce(self): - """Read energy-dependent elastic scattering cross sections. - """ - # Determine if ITCE block exists - idx = self._jxs[4] - if idx == 0: - return - - # Read values - n_energies = int(self._xss[idx]) - energy = self._xss[idx+1 : idx+1+n_energies] - P = self._xss[idx+1+n_energies : idx+1+2*n_energies] - - if self._nxs[5] == 4: - self.elastic_xs = CoherentElastic(energy, P) - else: - self.elastic_xs = Tabulated1D(energy, P) - - def _read_itxe(self): - """Read coupled energy/angle distributions for inelastic scattering. - """ - # Determine number of energies and angles - NE_in = len(self.inelastic_xs) - NE_out = self._nxs[4] - - if self._nxs[7] == 0: - self.secondary_mode = 'equal' - elif self._nxs[7] == 1: - self.secondary_mode = 'skewed' - elif self._nxs[7] == 2: - self.secondary_mode = 'continuous' - - if self.secondary_mode in ('equal', 'skewed'): - NMU = self._nxs[3] - idx = self._jxs[3] - self.inelastic_e_out = self._xss[idx:idx+NE_in*NE_out*(NMU+2):NMU+2] - self.inelastic_e_out.shape = (NE_in, NE_out) - - self.inelastic_mu_out = self._xss[idx:idx+NE_in*NE_out*(NMU+2)] - self.inelastic_mu_out.shape = (NE_in, NE_out, NMU+2) - self.inelastic_mu_out = self.inelastic_mu_out[:, :, 1:] - else: - NMU = self._nxs[3] - 1 - idx = self._jxs[3] - locc = self._xss[idx:idx + NE_in].astype(int) - NE_out = self._xss[idx + NE_in:idx + 2*NE_in].astype(int) - energy_out = [] - mu_out = [] - for i in range(NE_in): - idx = locc[i] - - # Outgoing energy distribution for incoming energy i - e = self._xss[idx + 1:idx + 1 + NE_out[i]*(NMU + 3):NMU + 3] - p = self._xss[idx + 2:idx + 2 + NE_out[i]*(NMU + 3):NMU + 3] - c = self._xss[idx + 3:idx + 3 + NE_out[i]*(NMU + 3):NMU + 3] - eout_i = Tabular(e, p, 'linear-linear', ignore_negative=True) - eout_i.c = c - - # Outgoing angle distribution for each (incoming, outgoing) energy pair - mu_i = [] - for j in range(NE_out[i]): - mu = self._xss[idx + 4:idx + 4 + NMU] - p_mu = 1./NMU*np.ones(NMU) - mu_ij = Discrete(mu, p_mu) - mu_ij.c = np.cumsum(p_mu) - mu_i.append(mu_ij) - idx += 3 + NMU - - energy_out.append(eout_i) - mu_out.append(mu_i) - - # Create correlated angle-energy distribution - breakpoints = [NE_in] - interpolation = [2] - energy = self.inelastic_xs.x - self.inelastic_dist = CorrelatedAngleEnergy( - breakpoints, interpolation, energy, energy_out, mu_out) - - def _read_itca(self): - """Read angular distributions for elastic scattering. - """ - NMU = self._nxs[6] - if self._jxs[4] == 0 or NMU == -1: - return - idx = self._jxs[6] - - NE = len(self.elastic_xs) - self.elastic_mu_out = self._xss[idx:idx+NE*NMU] - self.elastic_mu_out.shape = (NE, NMU) - - def export_to_hdf5(self, path, name): - """Export table to an HDF5 file. - - Parameters - ---------- - path : str - Path to write HDF5 file to - name : str - Name of compound (used as first group in HDF5 file) - - """ - - f = h5py.File(path, 'a') - - # Write basic data - g = f.create_group('{}.{}'.format(name, self.name.split('.')[1])) - g.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio - g.attrs['temperature'] = self.temperature - g.attrs['zaids'] = self.zaids - - # Write thermal elastic scattering - if self.elastic_xs is not None: - elastic_group = g.create_group('elastic') - self.elastic_xs.to_hdf5(elastic_group, 'xs') - if self.elastic_mu_out is not None: - elastic_group.create_dataset('mu_out', data=self.elastic_mu_out) - - # Write thermal inelastic scattering - if self.inelastic_xs is not None: - inelastic_group = g.create_group('inelastic') - self.inelastic_xs.to_hdf5(inelastic_group, 'xs') - inelastic_group.attrs['secondary_mode'] = np.string_(self.secondary_mode) - if self.secondary_mode in ('equal', 'skewed'): - inelastic_group.create_dataset('energy_out', data=self.inelastic_e_out) - inelastic_group.create_dataset('mu_out', data=self.inelastic_mu_out) - elif self.secondary_mode == 'continuous': - self.inelastic_dist.to_hdf5(inelastic_group) - - @classmethod - def from_hdf5(self, group): - """Generate thermal scattering data from HDF5 group - - Parameters - ---------- - group : h5py.Group - HDF5 group to read from - - Returns - ------- - openmc.data.SabTable - Neutron thermal scattering data - - """ - name = group.name[1:] - atomic_weight_ratio = group.attrs['atomic_weight_ratio'] - temperature = group.attrs['temperature'] - table = SabTable(name, atomic_weight_ratio, temperature) - table.zaids = group.attrs['zaids'] - - # Read thermal elastic scattering - if 'elastic' in group: - elastic_group = group['elastic'] - - # Cross section - elastic_xs_type = elastic_group['xs'].attrs['type'].decode() - if elastic_xs_type == 'tab1': - table.elastic_xs = Tabulated1D.from_hdf5(elastic_group['xs']) - elif elastic_xs_type == 'bragg': - table.elastic_xs = CoherentElastic.from_hdf5(elastic_group['xs']) - - # Angular distribution - if 'mu_out' in elastic_group: - table.elastic_mu_out = elastic_group['mu_out'].value - - # Read thermal inelastic scattering - if 'inelastic' in group: - inelastic_group = group['inelastic'] - table.secondary_mode = inelastic_group.attrs['secondary_mode'].decode() - table.inelastic_xs = Tabulated1D.from_hdf5(inelastic_group['xs']) - if table.secondary_mode in ('equal', 'skewed'): - table.inelastic_e_out = inelastic_group['energy_out'] - table.inelastic_mu_out = inelastic_group['mu_out'] - elif table.secondary_mode == 'continuous': - table.inelastic_dist = AngleEnergy.from_hdf5(inelastic_group) - - return table - - -class Reaction(object): - """Reaction(MT, table=None) - - A Reaction object represents a single reaction channel for a nuclide with - an associated cross section and, if present, a secondary angle and energy - distribution. These objects are stored within the ``reactions`` attribute on - subclasses of Table, e.g. NeutronTable. - - Parameters - ---------- - MT : int - The ENDF MT number for this reaction. On occasion, MCNP uses MT numbers - that don't correspond exactly to the ENDF specification. - table : openmc.data.ace.Table - The ACE table which contains this reaction. This is useful if data on - the parent nuclide is needed (for instance, the energy grid at which - cross sections are tabulated) - - Attributes - ---------- - center_of_mass : bool - Indicates whether scattering kinematics should be performed in the - center-of-mass or laboratory reference frame. - grid above the threshold value in barns. - MT : int - The ENDF MT number for this reaction. - Q_value : float - The Q-value of this reaction in MeV. - table : openmc.data.ace.Table - The ACE table which contains this reaction. - threshold : float - Threshold of the reaction in MeV - threshold_idx : int - The index on the energy grid corresponding to the threshold of this - reaction. - xs : openmc.data.Tabulated1D - Microscopic cross section for this reaction as a function of incident - energy - products : Iterable of openmc.data.Product - Reaction products - - """ - - def __init__(self, MT, table=None): - self.center_of_mass = True - self.table = table - self.MT = MT - self.Q_value = 0. - self.threshold_idx = 0 - self._xs = None - self.products = [] - - def __repr__(self): - if self.MT in reaction_name: - return "".format(self.MT, reaction_name[self.MT]) - else: - return "".format(self.MT) - - @property - def center_of_mass(self): - return self._center_of_mass - - @property - def products(self): - return self._products - - @property - def threshold(self): - return self.xs.x[0] - - @property - def xs(self): - return self._xs - - @center_of_mass.setter - def center_of_mass(self, center_of_mass): - cv.check_type('center of mass', center_of_mass, (bool, np.bool_)) - self._center_of_mass = center_of_mass - - @products.setter - def products(self, products): - cv.check_type('reaction products', products, Iterable, Product) - self._products = products - - @xs.setter - def xs(self, xs): - cv.check_type('reaction cross section', xs, Tabulated1D) - for y in xs.y: - cv.check_greater_than('reaction cross section', y, 0.0, True) - self._xs = xs - - def to_hdf5(self, group): - """Write reaction to an HDF5 group - - Parameters - ---------- - group : h5py.Group - HDF5 group to write to - - """ - - group.attrs['MT'] = self.MT - if self.MT in reaction_name: - group.attrs['label'] = np.string_(reaction_name[self.MT]) - else: - group.attrs['label'] = np.string_(self.MT) - group.attrs['Q_value'] = self.Q_value - group.attrs['threshold_idx'] = self.threshold_idx + 1 - group.attrs['center_of_mass'] = 1 if self.center_of_mass else 0 - group.attrs['n_product'] = len(self.products) - if self.xs is not None: - group.create_dataset('xs', data=self.xs.y) - for i, p in enumerate(self.products): - pgroup = group.create_group('product_{}'.format(i)) - p.to_hdf5(pgroup) - - @classmethod - def from_hdf5(cls, group, table): - """Generate reaction from an HDF5 group - - Parameters - ---------- - group : h5py.Group - HDF5 group to write to - - Returns - ------- - openmc.data.ace.Reaction - Reaction data - - """ - MT = group.attrs['MT'] - rxn = cls(MT) - rxn.table = table - rxn.Q_value = group.attrs['Q_value'] - rxn.threshold_idx = group.attrs['threshold_idx'] - 1 - rxn.center_of_mass = bool(group.attrs['center_of_mass']) - - # Read cross section - if 'xs' in group: - xs = group['xs'].value - rxn.xs = Tabulated1D(table.energy, xs) - - # Read reaction products - n_product = group.attrs['n_product'] - products = [] - for i in range(n_product): - pgroup = group['product_{}'.format(i)] - products.append(Product.from_hdf5(pgroup)) - rxn.products = products - - return rxn - - -table_types = { - "c": NeutronTable, - "t": SabTable} + return "".format(self.name) diff --git a/openmc/data/angle_distribution.py b/openmc/data/angle_distribution.py index 5f19450bc..3f086a99c 100644 --- a/openmc/data/angle_distribution.py +++ b/openmc/data/angle_distribution.py @@ -4,7 +4,7 @@ from numbers import Real import numpy as np import openmc.checkvalue as cv -from openmc.stats import Univariate, Tabular +from openmc.stats import Univariate, Tabular, Uniform from .container import interpolation_scheme @@ -132,3 +132,69 @@ class AngleDistribution(object): mu.append(mu_i) return cls(energy, mu) + + @classmethod + def from_ace(cls, ace, location_dist, location_start): + """Generate an angular distribution from ACE data + + Parameters + ---------- + ace : openmc.data.ace.Table + ACE table to read from + location_dist : int + Index in the XSS array corresponding to the start of a block, + e.g. JXS(9). + location_start : int + Index in the XSS array corresponding to the start of an angle + distribution array + + Returns + ------- + openmc.data.AngleDistribution + Angular distribution + + """ + # Set starting index for angle distribution + idx = location_dist + location_start - 1 + + # Number of energies at which angular distributions are tabulated + n_energies = int(ace.xss[idx]) + idx += 1 + + # Incoming energy grid + energy = ace.xss[idx:idx + n_energies] + idx += n_energies + + # Read locations for angular distributions + lc = ace.xss[idx:idx + n_energies].astype(int) + idx += n_energies + + mu = [] + for i in range(n_energies): + if lc[i] > 0: + # Equiprobable 32 bin distribution + idx = location_dist + abs(lc[i]) - 1 + cos = ace.xss[idx:idx + 33] + pdf = np.zeros(33) + pdf[:32] = 1.0/(32.0*np.diff(cos)) + cdf = np.linspace(0.0, 1.0, 33) + + mu_i = Tabular(cos, pdf, 'histogram', ignore_negative=True) + mu_i.c = cdf + elif lc[i] < 0: + # Tabular angular distribution + idx = location_dist + abs(lc[i]) - 1 + intt = int(ace.xss[idx]) + n_points = int(ace.xss[idx + 1]) + data = ace.xss[idx + 2:idx + 2 + 3*n_points] + data.shape = (3, n_points) + + mu_i = Tabular(data[0], data[1], interpolation_scheme[intt]) + mu_i.c = data[2] + else: + # Isotropic angular distribution + mu_i = Uniform(-1., 1.) + + mu.append(mu_i) + + return cls(energy, mu) diff --git a/openmc/data/angle_energy.py b/openmc/data/angle_energy.py index fc17d7bec..ff9f41a44 100644 --- a/openmc/data/angle_energy.py +++ b/openmc/data/angle_energy.py @@ -38,3 +38,73 @@ class AngleEnergy(object): return openmc.data.KalbachMann.from_hdf5(group) elif dist_type == 'nbody': return openmc.data.NBodyPhaseSpace.from_hdf5(group) + + @staticmethod + def from_ace(ace, location_dist, location_start, rx=None): + """Generate an AngleEnergy object from ACE data + + Parameters + ---------- + ace : openmc.data.ace.Table + ACE table to read from + location_dist : int + Index in the XSS array corresponding to the start of a block, + e.g. JXS(11) for the the DLW block. + location_start : int + Index in the XSS array corresponding to the start of an energy + distribution array + rx : Reaction + Reaction this energy distribution will be associated with + + Returns + ------- + distribution : openmc.data.AngleEnergy + Secondary angle-energy distribution + + """ + # Set starting index for energy distribution + idx = location_dist + location_start - 1 + + law = int(ace.xss[idx + 1]) + location_data = int(ace.xss[idx + 2]) + + # Position index for reading law data + idx = location_dist + location_data - 1 + + # Parse energy distribution data + if law == 2: + distribution = openmc.data.UncorrelatedAngleEnergy() + distribution.energy = openmc.data.DiscretePhoton.from_ace(ace, idx) + elif law in (3, 33): + distribution = openmc.data.UncorrelatedAngleEnergy() + distribution.energy = openmc.data.LevelInelastic.from_ace(ace, idx) + elif law == 4: + distribution = openmc.data.UncorrelatedAngleEnergy() + distribution.energy = openmc.data.ContinuousTabular.from_ace( + ace, idx, location_dist) + elif law == 5: + distribution = openmc.data.UncorrelatedAngleEnergy() + distribution.energy = openmc.data.GeneralEvaporation.from_ace(ace, idx) + elif law == 7: + distribution = openmc.data.UncorrelatedAngleEnergy() + distribution.energy = openmc.data.MaxwellEnergy.from_ace(ace, idx) + elif law == 9: + distribution = openmc.data.UncorrelatedAngleEnergy() + distribution.energy = openmc.data.Evaporation.from_ace(ace, idx) + elif law == 11: + distribution = openmc.data.UncorrelatedAngleEnergy() + distribution.energy = openmc.data.WattEnergy.from_ace(ace, idx) + elif law == 44: + distribution = openmc.data.KalbachMann.from_ace( + ace, idx, location_dist) + elif law == 61: + distribution = openmc.data.CorrelatedAngleEnergy.from_ace( + ace, idx, location_dist) + elif law == 66: + distribution = openmc.data.NBodyPhaseSpace.from_ace( + ace, idx, rx.q_value) + else: + raise IOError("Unsupported ACE secondary energy " + "distribution law {0}".format(law)) + + return distribution diff --git a/openmc/data/container.py b/openmc/data/container.py index 4bd293c36..b2b323ddd 100644 --- a/openmc/data/container.py +++ b/openmc/data/container.py @@ -267,3 +267,42 @@ class Tabulated1D(object): breakpoints = dataset.attrs['breakpoints'] interpolation = dataset.attrs['interpolation'] return cls(x, y, breakpoints, interpolation) + + @classmethod + def from_ace(cls, ace, idx=0): + """Create a Tabulated1D object from an ACE table. + + Parameters + ---------- + ace : openmc.data.ace.Table + An ACE table + idx : int + Offset to read from in XSS array (default of zero) + + Returns + ------- + openmc.data.Tabulated1D + Tabulated data object + + """ + + # Get number of regions and pairs + n_regions = int(ace.xss[idx]) + n_pairs = int(ace.xss[idx + 1 + 2*n_regions]) + + # Get interpolation information + idx += 1 + if n_regions > 0: + breakpoints = ace.xss[idx:idx + n_regions].astype(int) + interpolation = ace.xss[idx + n_regions:idx + 2*n_regions].astype(int) + else: + # 0 regions implies linear-linear interpolation by default + breakpoints = np.array([n_pairs]) + interpolation = np.array([2]) + + # Get (x,y) pairs + idx += 2*n_regions + 1 + x = ace.xss[idx:idx + n_pairs] + y = ace.xss[idx + n_pairs:idx + 2*n_pairs] + + return Tabulated1D(x, y, breakpoints, interpolation) diff --git a/openmc/data/correlated.py b/openmc/data/correlated.py index 58a90da9d..f82a28b75 100644 --- a/openmc/data/correlated.py +++ b/openmc/data/correlated.py @@ -203,8 +203,8 @@ class CorrelatedAngleEnergy(AngleEnergy): """ interp_data = group['energy'].attrs['interpolation'] - energy_breakpoints = interp_data[0,:] - energy_interpolation = interp_data[1,:] + energy_breakpoints = interp_data[0, :] + energy_interpolation = interp_data[1, :] energy = group['energy'].value offsets = group['energy_out'].attrs['offsets'] @@ -289,3 +289,116 @@ class CorrelatedAngleEnergy(AngleEnergy): return cls(energy_breakpoints, energy_interpolation, energy, energy_out, mu) + + @classmethod + def from_ace(cls, ace, idx, ldis): + """Generate correlated angle-energy distribution from ACE data + + Parameters + ---------- + ace : openmc.data.ace.Table + ACE table to read from + idx : int + Index in XSS array of the start of the energy distribution data + (LDIS + LOCC - 1) + ldis : int + Index in XSS array of the start of the energy distribution block + (e.g. JXS[11]) + + Returns + ------- + openmc.data.CorrelatedAngleEnergy + Correlated angle-energy distribution + + """ + # Read number of interpolation regions and incoming energies + n_regions = int(ace.xss[idx]) + n_energy_in = int(ace.xss[idx + 1 + 2*n_regions]) + + # Get interpolation information + idx += 1 + if n_regions > 0: + breakpoints = ace.xss[idx:idx + n_regions].astype(int) + interpolation = ace.xss[idx + n_regions:idx + 2*n_regions].astype(int) + else: + breakpoints = np.array([n_energy_in]) + interpolation = np.array([2]) + + # Incoming energies at which distributions exist + idx += 2*n_regions + 1 + energy = ace.xss[idx:idx + n_energy_in] + + # Location of distributions + idx += n_energy_in + loc_dist = ace.xss[idx:idx + n_energy_in].astype(int) + + # Initialize list of distributions + energy_out = [] + mu = [] + + # Read each outgoing energy distribution + for i in range(n_energy_in): + idx = ldis + loc_dist[i] - 1 + + # intt = interpolation scheme (1=hist, 2=lin-lin) + INTTp = int(ace.xss[idx]) + intt = INTTp % 10 + n_discrete_lines = (INTTp - intt)//10 + if intt not in (1, 2): + warn("Interpolation scheme for continuous tabular distribution " + "is not histogram or linear-linear.") + intt = 2 + + # Secondary energy distribution + n_energy_out = int(ace.xss[idx + 1]) + data = ace.xss[idx + 2:idx + 2 + 4*n_energy_out] + data.shape = (4, n_energy_out) + + # Create continuous distribution + eout_continuous = Tabular(data[0][n_discrete_lines:], + data[1][n_discrete_lines:], + interpolation_scheme[intt], + ignore_negative=True) + eout_continuous.c = data[2][n_discrete_lines:] + + # If discrete lines are present, create a mixture distribution + if n_discrete_lines > 0: + eout_discrete = Discrete(data[0][:n_discrete_lines], + data[1][:n_discrete_lines]) + eout_discrete.c = data[2][:n_discrete_lines] + if n_discrete_lines == n_energy_out: + eout_i = eout_discrete + else: + p_discrete = min(sum(eout_discrete.p), 1.0) + eout_i = Mixture([p_discrete, 1. - p_discrete], + [eout_discrete, eout_continuous]) + else: + eout_i = eout_continuous + + energy_out.append(eout_i) + + lc = data[3].astype(int) + + # Secondary angular distributions + mu_i = [] + for j in range(n_energy_out): + if lc[j] > 0: + idx = ldis + abs(lc[j]) - 1 + + intt = int(ace.xss[idx]) + n_cosine = int(ace.xss[idx + 1]) + data = ace.xss[idx + 2:idx + 2 + 3*n_cosine] + data.shape = (3, n_cosine) + + mu_ij = Tabular(data[0], data[1], interpolation_scheme[intt]) + mu_ij.c = data[2] + else: + # Isotropic distribution + mu_ij = Uniform(-1., 1.) + + mu_i.append(mu_ij) + + # Add cosine distributions for this incoming energy to list + mu.append(mu_i) + + return cls(breakpoints, interpolation, energy, energy_out, mu) diff --git a/openmc/data/energy_distribution.py b/openmc/data/energy_distribution.py index afbc39014..09b0c6ebf 100644 --- a/openmc/data/energy_distribution.py +++ b/openmc/data/energy_distribution.py @@ -1,6 +1,7 @@ from abc import ABCMeta, abstractmethod from collections import Iterable from numbers import Integral, Real +from warnings import warn import h5py import numpy as np @@ -200,6 +201,33 @@ class MaxwellEnergy(EnergyDistribution): u = group.attrs['u'] return cls(theta, u) + @classmethod + def from_ace(cls, ace, idx=0): + """Create a Maxwell distribution from an ACE table + + Parameters + ---------- + ace : openmc.data.ace.Table + An ACE table + idx : int + Offset to read from in XSS array (default of zero) + + Returns + ------- + openmc.data.MaxwellEnergy + Maxwell distribution + + """ + # Read nuclear temperature + theta = Tabulated1D.from_ace(ace, idx) + + # Restriction energy + nr = int(ace.xss[idx]) + ne = int(ace.xss[idx + 1 + 2*nr]) + u = ace.xss[idx + 2 + 2*nr + 2*ne] + + return cls(theta, u) + class Evaporation(EnergyDistribution): r"""Evaporation spectrum represented as @@ -273,13 +301,40 @@ class Evaporation(EnergyDistribution): Returns ------- openmc.data.Evaporation - Evaporation spectrum distribution + Evaporation spectrum """ theta = Tabulated1D.from_hdf5(group['theta']) u = group.attrs['u'] return cls(theta, u) + @classmethod + def from_ace(cls, ace, idx=0): + """Create an evaporation spectrum from an ACE table + + Parameters + ---------- + ace : openmc.data.ace.Table + An ACE table + idx : int + Offset to read from in XSS array (default of zero) + + Returns + ------- + openmc.data.Evaporation + Evaporation spectrum + + """ + # Read nuclear temperature + theta = Tabulated1D.from_ace(ace, idx) + + # Restriction energy + nr = int(ace.xss[idx]) + ne = int(ace.xss[idx + 1 + 2*nr]) + u = ace.xss[idx + 2 + 2*nr + 2*ne] + + return cls(theta, u) + class WattEnergy(EnergyDistribution): r"""Energy-dependent Watt spectrum represented as @@ -375,6 +430,45 @@ class WattEnergy(EnergyDistribution): u = group.attrs['u'] return cls(a, b, u) + @classmethod + def from_ace(cls, ace, idx): + """Create a Watt fission spectrum from an ACE table + + Parameters + ---------- + ace : openmc.data.ace.Table + An ACE table + idx : int + Offset to read from in XSS array (default of zero) + + Returns + ------- + openmc.data.WattEnergy + Watt fission spectrum + + """ + # Energy-dependent a parameter + a = Tabulated1D.from_ace(ace, idx) + + # Advance index + nr = int(ace.xss[idx]) + ne = int(ace.xss[idx + 1 + 2*nr]) + idx += 2 + 2*nr + 2*ne + + # Energy-dependent b parameter + b = Tabulated1D.from_ace(ace, idx) + + # Advance index + nr = int(ace.xss[idx]) + ne = int(ace.xss[idx + 1 + 2*nr]) + idx += 2 + 2*nr + 2*ne + + # Restriction energy + u = ace.xss[idx] + + return cls(a, b, u) + + class MadlandNix(EnergyDistribution): r"""Energy-dependent fission neutron spectrum (Madland and Nix) given in ENDF MF=5, LF=12 represented as @@ -574,6 +668,27 @@ class DiscretePhoton(EnergyDistribution): awr = group.attrs['atomic_weight_ratio'] return cls(primary_flag, energy, awr) + @classmethod + def from_ace(cls, ace, idx): + """Generate discrete photon energy distribution from an ACE table + + Parameters + ---------- + ace : openmc.data.ace.Table + An ACE table + idx : int + Offset to read from in XSS array (default of zero) + + Returns + ------- + openmc.data.DiscretePhoton + Discrete photon energy distribution + + """ + primary_flag = int(ace.xss[idx]) + energy = ace.xss[idx + 1] + return cls(primary_flag, energy, ace.atomic_weight_ratio) + class LevelInelastic(EnergyDistribution): r"""Level inelastic scattering @@ -650,6 +765,26 @@ class LevelInelastic(EnergyDistribution): mass_ratio = group.attrs['mass_ratio'] return cls(threshold, mass_ratio) + @classmethod + def from_ace(cls, ace, idx): + """Generate level inelastic distribution from an ACE table + + Parameters + ---------- + ace : openmc.data.ace.Table + An ACE table + idx : int + Offset to read from in XSS array (default of zero) + + Returns + ------- + openmc.data.LevelInelastic + Level inelastic scattering distribution + + """ + threshold, mass_ratio = ace.xss[idx:idx + 2] + return cls(threshold, mass_ratio) + class ContinuousTabular(EnergyDistribution): """Continuous tabular distribution @@ -802,8 +937,8 @@ class ContinuousTabular(EnergyDistribution): """ interp_data = group['energy'].attrs['interpolation'] - energy_breakpoints = interp_data[0,:] - energy_interpolation = interp_data[1,:] + energy_breakpoints = interp_data[0, :] + energy_interpolation = interp_data[1, :] energy = group['energy'].value data = group['distribution'] @@ -848,3 +983,89 @@ class ContinuousTabular(EnergyDistribution): return cls(energy_breakpoints, energy_interpolation, energy, energy_out) + + @classmethod + def from_ace(cls, ace, idx, ldis): + """Generate continuous tabular energy distribution from ACE data + + Parameters + ---------- + ace : openmc.data.ace.Table + ACE table to read from + idx : int + Index in XSS array of the start of the energy distribution data + (LDIS + LOCC - 1) + ldis : int + Index in XSS array of the start of the energy distribution block + (e.g. JXS[11]) + + Returns + ------- + openmc.data.ContinuousTabular + Continuous tabular energy distribution + + """ + # Read number of interpolation regions and incoming energies + n_regions = int(ace.xss[idx]) + n_energy_in = int(ace.xss[idx + 1 + 2*n_regions]) + + # Get interpolation information + idx += 1 + if n_regions > 0: + breakpoints = ace.xss[idx:idx + n_regions].astype(int) + interpolation = ace.xss[idx + n_regions:idx + 2*n_regions].astype(int) + else: + breakpoints = np.array([n_energy_in]) + interpolation = np.array([2]) + + # Incoming energies at which distributions exist + idx += 2*n_regions + 1 + energy = ace.xss[idx:idx + n_energy_in] + + # Location of distributions + idx += n_energy_in + loc_dist = ace.xss[idx:idx + n_energy_in].astype(int) + + # Initialize variables + energy_out = [] + + # Read each outgoing energy distribution + for i in range(n_energy_in): + idx = ldis + loc_dist[i] - 1 + + # intt = interpolation scheme (1=hist, 2=lin-lin) + INTTp = int(ace.xss[idx]) + intt = INTTp % 10 + n_discrete_lines = (INTTp - intt)//10 + if intt not in (1, 2): + warn("Interpolation scheme for continuous tabular distribution " + "is not histogram or linear-linear.") + intt = 2 + + n_energy_out = int(ace.xss[idx + 1]) + data = ace.xss[idx + 2:idx + 2 + 3*n_energy_out] + data.shape = (3, n_energy_out) + + # Create continuous distribution + eout_continuous = Tabular(data[0][n_discrete_lines:], + data[1][n_discrete_lines:], + interpolation_scheme[intt]) + eout_continuous.c = data[2][n_discrete_lines:] + + # If discrete lines are present, create a mixture distribution + if n_discrete_lines > 0: + eout_discrete = Discrete(data[0][:n_discrete_lines], + data[1][:n_discrete_lines]) + eout_discrete.c = data[2][:n_discrete_lines] + if n_discrete_lines == n_energy_out: + eout_i = eout_discrete + else: + p_discrete = min(sum(eout_discrete.p), 1.0) + eout_i = Mixture([p_discrete, 1. - p_discrete], + [eout_discrete, eout_continuous]) + else: + eout_i = eout_continuous + + energy_out.append(eout_i) + + return cls(breakpoints, interpolation, energy, energy_out) diff --git a/openmc/data/kalbach_mann.py b/openmc/data/kalbach_mann.py index a8fc8d4c6..7899ae2a6 100644 --- a/openmc/data/kalbach_mann.py +++ b/openmc/data/kalbach_mann.py @@ -197,8 +197,8 @@ class KalbachMann(AngleEnergy): """ interp_data = group['energy'].attrs['interpolation'] - energy_breakpoints = interp_data[0,:] - energy_interpolation = interp_data[1,:] + energy_breakpoints = interp_data[0, :] + energy_interpolation = interp_data[1, :] energy = group['energy'].value data = group['distribution'] @@ -251,3 +251,93 @@ class KalbachMann(AngleEnergy): return cls(energy_breakpoints, energy_interpolation, energy, energy_out, precompound, slope) + + @classmethod + def from_ace(cls, ace, idx, ldis): + """Generate Kalbach-Mann energy-angle distribution from ACE data + + Parameters + ---------- + ace : openmc.data.ace.Table + ACE table to read from + idx : int + Index in XSS array of the start of the energy distribution data + (LDIS + LOCC - 1) + ldis : int + Index in XSS array of the start of the energy distribution block + (e.g. JXS[11]) + + Returns + ------- + openmc.data.KalbachMann + Kalbach-Mann energy-angle distribution + + """ + # Read number of interpolation regions and incoming energies + n_regions = int(ace.xss[idx]) + n_energy_in = int(ace.xss[idx + 1 + 2*n_regions]) + + # Get interpolation information + idx += 1 + if n_regions > 0: + breakpoints = ace.xss[idx:idx + n_regions].astype(int) + interpolation = ace.xss[idx + n_regions:idx + 2*n_regions].astype(int) + else: + breakpoints = np.array([n_energy_in]) + interpolation = np.array([2]) + + # Incoming energies at which distributions exist + idx += 2*n_regions + 1 + energy = ace.xss[idx:idx + n_energy_in] + + # Location of distributions + idx += n_energy_in + loc_dist = ace.xss[idx:idx + n_energy_in].astype(int) + + # Initialize variables + energy_out = [] + km_r = [] + km_a = [] + + # Read each outgoing energy distribution + for i in range(n_energy_in): + idx = ldis + loc_dist[i] - 1 + + # intt = interpolation scheme (1=hist, 2=lin-lin) + INTTp = int(ace.xss[idx]) + intt = INTTp % 10 + n_discrete_lines = (INTTp - intt)//10 + if intt not in (1, 2): + warn("Interpolation scheme for continuous tabular distribution " + "is not histogram or linear-linear.") + intt = 2 + + n_energy_out = int(ace.xss[idx + 1]) + data = ace.xss[idx + 2:idx + 2 + 5*n_energy_out] + data.shape = (5, n_energy_out) + + # Create continuous distribution + eout_continuous = Tabular(data[0][n_discrete_lines:], + data[1][n_discrete_lines:], + interpolation_scheme[intt]) + eout_continuous.c = data[2][n_discrete_lines:] + + # If discrete lines are present, create a mixture distribution + if n_discrete_lines > 0: + eout_discrete = Discrete(data[0][:n_discrete_lines], + data[1][:n_discrete_lines]) + eout_discrete.c = data[2][:n_discrete_lines] + if n_discrete_lines == n_energy_out: + eout_i = eout_discrete + else: + p_discrete = min(sum(eout_discrete.p), 1.0) + eout_i = Mixture([p_discrete, 1. - p_discrete], + [eout_discrete, eout_continuous]) + else: + eout_i = eout_continuous + + energy_out.append(eout_i) + km_r.append(Tabulated1D(data[0], data[3])) + km_a.append(Tabulated1D(data[0], data[4])) + + return cls(breakpoints, interpolation, energy, energy_out, km_r, km_a) diff --git a/openmc/data/library.py b/openmc/data/library.py new file mode 100644 index 000000000..58f32609e --- /dev/null +++ b/openmc/data/library.py @@ -0,0 +1,47 @@ +import os +import xml.etree.ElementTree as ET + +import h5py + +from openmc.clean_xml import clean_xml_indentation + +class DataLibrary(object): + def __init__(self): + self.libraries = [] + + def register_file(self, filename, filetype='neutron'): + h5file = h5py.File(filename, 'r') + + materials = [] + for name, group in h5file.items(): + materials.append(name) + + library = {'path': filename, 'type': filetype, 'materials': materials} + self.libraries.append(library) + + def export_to_xml(self, path='cross_sections.xml'): + root = ET.Element('cross_sections') + + # Determine common directory for library paths + common_dir = os.path.commonpath([lib['path'] for lib in self.libraries]) + if common_dir == '': + common_dir = '.' + + directory = os.path.relpath(common_dir, os.path.dirname(path)) + if directory != '.': + dir_element = ET.SubElement(root, "directory") + dir_element.text = directory + + for library in self.libraries: + lib_element = ET.SubElement(root, "library") + lib_element.set('materials', ' '.join(library['materials'])) + lib_element.set('path', os.path.relpath(library['path'], common_dir)) + lib_element.set('type', library['type']) + + # Clean the indentation to be user-readable + clean_xml_indentation(root) + + # Write XML file + tree = ET.ElementTree(root) + tree.write(path, xml_declaration=True, encoding='utf-8', + method='xml') diff --git a/openmc/data/nbody.py b/openmc/data/nbody.py index 4fa3d05b2..a34380dc9 100644 --- a/openmc/data/nbody.py +++ b/openmc/data/nbody.py @@ -116,3 +116,27 @@ class NBodyPhaseSpace(AngleEnergy): awr = group.attrs['atomic_weight_ratio'] q_value = group.attrs['q_value'] return cls(total_mass, n_particles, awr, q_value) + + @classmethod + def from_ace(cls, ace, idx, q_value): + """Generate N-body phase space distribution from ACE data + + Parameters + ---------- + ace : openmc.data.ace.Table + ACE table to read from + idx : int + Index in XSS array of the start of the energy distribution data + (LDIS + LOCC - 1) + q_value : float + Q-value for reaction in MeV + + Returns + ------- + openmc.data.NBodyPhaseSpace + N-body phase space distribution + + """ + n_particles = int(ace.xss[idx]) + total_mass = ace.xss[idx + 1] + return cls(total_mass, n_particles, ace.atomic_weight_ratio, q_value) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py new file mode 100644 index 000000000..ad5e370f1 --- /dev/null +++ b/openmc/data/neutron.py @@ -0,0 +1,405 @@ +from __future__ import division, unicode_literals +import io +import sys +from warnings import warn +from collections import OrderedDict, Iterable, Mapping +from copy import deepcopy +from numbers import Integral, Real +import sys + +import numpy as np +from numpy.polynomial import Polynomial +import h5py + +from . import atomic_number, atomic_symbol +from .ace import Table, get_table +from .container import Tabulated1D +from .energy_distribution import * +from .product import Product +from .reaction import Reaction, _get_photon_products +from .thermal import CoherentElastic +from .urr import ProbabilityTables +from openmc.stats import Tabular, Discrete, Uniform, Mixture +import openmc.checkvalue as cv + +if sys.version_info[0] >= 3: + basestring = str + + +class IncidentNeutron(object): + """Continuous-energy neutron interaction data. + + Instances of this class are not normally instantiated by the user but rather + created using the factory methods :meth:`IncidentNeutron.from_hdf5` and + :meth:`IncidentNeutron.from_ace`. + + Parameters + ---------- + name : str + Name of the table + atomic_number : int + Number of protons in the nucleus + mass_number : int + Number of nucleons in the nucleus + metastable : int + Metastable state of the nucleus. A value of zero indicates ground state. + atomic_weight_ratio : float + Atomic mass ratio of the target nuclide. + temperature : float + Temperature of the target nuclide in eV. + + Attributes + ---------- + atomic_number : int + Number of protons in the nucleus + atomic_symbol : str + Atomic symbol of the nuclide, e.g., 'Zr' + atomic_weight_ratio : float + Atomic weight ratio of the target nuclide. + energy : numpy.ndarray + The energy values (MeV) at which reaction cross-sections are tabulated. + mass_number : int + Number of nucleons in the nucleus + metastable : int + Metastable state of the nucleus. A value of zero indicates ground state. + name : str + ZAID identifier of the table, e.g. 92235.70c. + reactions : collections.OrderedDict + 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. + summed_reactions : collections.OrderedDict + Contains summed cross sections, e.g., the total cross section. The keys + are the MT values and the values are Reaction objects. + temperature : float + Temperature of the target nuclide in eV. + urr : None or openmc.data.ProbabilityTables + Unresolved resonance region probability tables + + """ + + def __init__(self, name, atomic_number, mass_number, metastable, + atomic_weight_ratio, temperature): + self.name = name + self.atomic_number = atomic_number + self.mass_number = mass_number + self.metastable = metastable + self.atomic_weight_ratio = atomic_weight_ratio + self.temperature = temperature + + self._energy = None + self.reactions = OrderedDict() + self.summed_reactions = OrderedDict() + self.urr = None + + def __repr__(self): + return "".format(self.name) + + def __iter__(self): + return iter(self.reactions.values()) + + @property + def name(self): + return self._name + + @property + def atomic_number(self): + return self._atomic_number + + @property + def mass_number(self): + return self._mass_number + + @property + def metastable(self): + return self._metastable + + @property + def atomic_weight_ratio(self): + return self._atomic_weight_ratio + + @property + def energy(self): + return self._energy + + @property + def temperature(self): + return self._temperature + + @property + def reactions(self): + return self._reactions + + @property + def summed_reactions(self): + return self._summed_reactions + + @property + def urr(self): + return self._urr + + @name.setter + def name(self, name): + cv.check_type('name', name, basestring) + self._name = name + + @property + def atomic_symbol(self): + return atomic_symbol[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) + self._atomic_number = atomic_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 + + @metastable.setter + def metastable(self, metastable): + cv.check_type('metastable', metastable, Integral) + cv.check_greater_than('metastable', metastable, 0, True) + self._metastable = metastable + + @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 + + @temperature.setter + def temperature(self, temperature): + cv.check_type('temperature', temperature, Real) + cv.check_greater_than('temperature', temperature, 0.0) + self._temperature = temperature + + @energy.setter + def energy(self, energy): + cv.check_type('energy grid', energy, Iterable, Real) + self._energy = energy + + @reactions.setter + def reactions(self, reactions): + cv.check_type('reactions', reactions, Mapping) + self._reactions = reactions + + @summed_reactions.setter + def summed_reactions(self, summed_reactions): + cv.check_type('summed reactions', summed_reactions, Mapping) + self._summed_reactions = summed_reactions + + @urr.setter + def urr(self, urr): + cv.check_type('probability tables', urr, + (ProbabilityTables, type(None))) + self._urr = urr + + def export_to_hdf5(self, path, mode='a'): + """Export table to an HDF5 file. + + Parameters + ---------- + path : str + Path to write HDF5 file to + mode : {'r', r+', 'w', 'x', 'a'} + Mode that is used to open the HDF5 file. This is the second argument + to the :class:`h5py.File` constructor. + + """ + + f = h5py.File(path, mode) + + # 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 + g.attrs['temperature'] = self.temperature + g.attrs['n_reaction'] = len(self.reactions) + + # Write energy grid + g.create_dataset('energy', data=self.energy) + + # Write reaction data + for i, rx in enumerate(self.reactions.values()): + rx_group = g.create_group('reaction_{}'.format(i)) + 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 is not None: + urr_group = g.create_group('urr') + self.urr.to_hdf5(urr_group) + + f.close() + + @classmethod + def from_hdf5(self, 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.ace.IncidentNeutron + Continuous-energy neutron interaction data + + """ + if isinstance(group_or_filename, h5py.Group): + group = group_or_filename + else: + h5file = h5py.File(group_or_filename, 'r') + 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'] + temperature = group.attrs['temperature'] + + data = IncidentNeutron(name, atomic_number, mass_number, metastable, + atomic_weight_ratio, temperature) + + # Read energy grid + data.energy = group['energy'].value + + # Read reaction data + n_reaction = group.attrs['n_reaction'] + + # Write reaction data + for i in range(n_reaction): + rx_group = group['reaction_{}'.format(i)] + rx = Reaction.from_hdf5(rx_group, data.energy) + data.reactions[rx.mt] = rx + + # Read total nu data if available + if 'total_nu' in rx_group: + tgroup = rx_group['total_nu'] + rx.derived_products = [Product.from_hdf5(tgroup)] + + # Read unresolved resonance probability tables + if 'urr' in group: + urr_group = group['urr'] + data.urr = ProbabilityTables.from_hdf5(urr_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 : 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. + + Returns + ------- + openmc.data.IncidentNeutron + Incident neutron continuous-energy data + + """ + 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('.') + zaid = int(zaid) + Z = zaid // 1000 + mass_number = zaid % 1000 + + if metastable_scheme == 'mcnp': + if zaid > 1000000: + # New SZA format + Z = Z % 1000 + if zaid == 1095242: + metastable = 0 + else: + metastable = zaid // 1000000 + else: + if zaid == 95242: + metastable = 1 + elif zaid == 95642: + metastable = 0 + else: + metastable = 1 if mass_number > 300 else 0 + elif metastable_scheme == 'nndc': + metastable = 1 if mass_number > 300 else 0 + + while mass_number > 3*Z: + mass_number -= 100 + + # Determine name for group + element = atomic_symbol[Z] + if metastable > 0: + name = '{}{}_m{}.{}'.format(element, mass_number, metastable, xs) + else: + name = '{}{}.{}'.format(element, mass_number, xs) + + data = IncidentNeutron(name, Z, mass_number, metastable, + ace.atomic_weight_ratio, ace.temperature) + + # Read energy grid + n_energy = ace.nxs[3] + energy = ace.xss[ace.jxs[1]:ace.jxs[1] + n_energy] + data.energy = energy + total_xs = ace.xss[ace.jxs[1] + n_energy:ace.jxs[1] + 2*n_energy] + absorption_xs = ace.xss[ace.jxs[1] + 2*n_energy:ace.jxs[1] + 3*n_energy] + + # Create summed reactions (total and absorption) + total = Reaction(1) + total.xs = Tabulated1D(energy, total_xs) + data.summed_reactions[1] = total + absorption = Reaction(27) + absorption.xs = Tabulated1D(energy, absorption_xs) + data.summed_reactions[27] = absorption + + # 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.reactions: + rx = Reaction(mt) + rx.products += _get_photon_products(ace, mt) + data.summed_reactions[mt] = rx + + # Read unresolved resonance probability tables + data.urr = ProbabilityTables.from_ace(ace) + + return data diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py new file mode 100644 index 000000000..40e0910ba --- /dev/null +++ b/openmc/data/reaction.py @@ -0,0 +1,512 @@ +from __future__ import division, unicode_literals +from collections import Iterable +from copy import deepcopy +from numbers import Real + +import numpy as np +from numpy.polynomial import Polynomial + +import openmc.checkvalue as cv +from openmc.stats import Uniform +from .angle_distribution import AngleDistribution +from .angle_energy import AngleEnergy +from .container import Tabulated1D +from .data import reaction_name +from .product import Product +from .uncorrelated import UncorrelatedAngleEnergy + + +def _get_fission_products(ace): + """Generate fission products from an ACE table + + Parameters + ---------- + ace : openmc.data.ace.Table + ACE table to read from + + Returns + ------- + products : list of openmc.data.Product + Prompt and delayed fission neutrons + derived_products : list of openmc.data.Product + "Total" fission neutron + + """ + # No NU block + if ace.jxs[2] == 0: + return None, None + + products = [] + derived_products = [] + + # Either prompt nu or total nu is given + if ace.xss[ace.jxs[2]] > 0: + whichnu = 'prompt' if ace.jxs[24] > 0 else 'total' + + neutron = Product('neutron') + neutron.emission_mode = whichnu + + idx = ace.jxs[2] + LNU = int(ace.xss[idx]) + if LNU == 1: + # Polynomial function form of nu + NC = int(ace.xss[idx+1]) + coefficients = ace.xss[idx+2 : idx+2+NC] + neutron.yield_ = Polynomial(coefficients) + elif LNU == 2: + # Tabular data form of nu + neutron.yield_ = Tabulated1D.from_ace(ace, idx + 1) + + products.append(neutron) + + # Both prompt nu and total nu + elif ace.xss[ace.jxs[2]] < 0: + # Read prompt neutron yield + prompt_neutron = Product('neutron') + prompt_neutron.emission_mode = 'prompt' + + idx = ace.jxs[2] + 1 + LNU = int(ace.xss[idx]) + if LNU == 1: + # Polynomial function form of nu + NC = int(ace.xss[idx+1]) + coefficients = ace.xss[idx+2 : idx+2+NC] + prompt_neutron.yield_ = Polynomial(coefficients) + elif LNU == 2: + # Tabular data form of nu + prompt_neutron.yield_ = Tabulated1D.from_ace(ace, idx + 1) + + # Read total neutron yield + total_neutron = Product('neutron') + total_neutron.emission_mode = 'total' + + idx = ace.jxs[2] + int(abs(ace.xss[ace.jxs[2]])) + 1 + LNU = int(ace.xss[idx]) + + if LNU == 1: + # Polynomial function form of nu + NC = int(ace.xss[idx+1]) + coefficients = ace.xss[idx+2 : idx+2+NC] + total_neutron.yield_ = Polynomial(coefficients) + elif LNU == 2: + # Tabular data form of nu + total_neutron.yield_ = Tabulated1D.from_ace(ace, idx + 1) + + products.append(prompt_neutron) + derived_products.append(total_neutron) + + # Check for delayed nu data + if ace.jxs[24] > 0: + yield_delayed = Tabulated1D.from_ace(ace, ace.jxs[24] + 1) + + # Delayed neutron precursor distribution + idx = ace.jxs[25] + n_group = ace.nxs[8] + total_group_probability = 0. + for i, group in enumerate(range(n_group)): + delayed_neutron = Product('neutron') + delayed_neutron.emission_mode = 'delayed' + delayed_neutron.decay_rate = ace.xss[idx] + + group_probability = Tabulated1D.from_ace(ace, idx + 1) + if np.all(group_probability.y == group_probability.y[0]): + delayed_neutron.yield_ = deepcopy(yield_delayed) + delayed_neutron.yield_.y *= group_probability.y[0] + total_group_probability += group_probability.y[0] + else: + raise NotImplementedError( + 'Delayed neutron with energy-dependent group probability') + + # Advance position + nr = int(ace.xss[idx + 1]) + ne = int(ace.xss[idx + 2 + 2*nr]) + idx += 3 + 2*nr + 2*ne + + # Energy distribution for delayed fission neutrons + location_start = int(ace.xss[ace.jxs[26] + group]) + delayed_neutron.distribution.append( + AngleEnergy.from_ace(ace, ace.jxs[27], location_start)) + + products.append(delayed_neutron) + + # Renormalize delayed neutron yields to reflect fact that in ACE + # file, the sum of the group probabilities is not exactly one + for product in products[1:]: + product.yield_.y /= total_group_probability + + return products, derived_products + + +def _get_photon_products(ace, mt): + """Generate photon products from an ACE table + + Parameters + ---------- + ace : openmc.data.ace.Table + ACE table to read from + mt : int + MT number for the desired reaction + + Returns + ------- + photons : list of openmc.Products + Photons produced from reaction with given MT + + """ + n_photon_reactions = ace.nxs[6] + photon_mts = ace.xss[ace.jxs[13]:ace.jxs[13] + + n_photon_reactions].astype(int) + + photons = [] + for i in range(n_photon_reactions): + # Determine corresponding reaction + neutron_mt = photon_mts[i] // 1000 + + # Restrict to photons that match the requested MT. Note that if the + # photon is assigned to MT=18 but the file splits fission into + # MT=19,20,21,38, we assign the photon product to each of the individual + # reactions + if neutron_mt == 18: + if mt not in (18, 19, 20, 21, 38): + continue + elif neutron_mt != mt: + continue + + # Create photon product and assign to reactions + photon = Product('photon') + + # ================================================================== + # Photon yield / production cross section + + loca = int(ace.xss[ace.jxs[14] + i]) + idx = ace.jxs[15] + loca - 1 + mftype = int(ace.xss[idx]) + idx += 1 + + if mftype in (12, 16): + # Yield data taken from ENDF File 12 or 6 + mtmult = int(ace.xss[idx]) + assert mtmult == neutron_mt + + # Read photon yield as function of energy + photon.yield_ = Tabulated1D.from_ace(ace, idx + 1) + + elif mftype == 13: + # Cross section data from ENDF File 13 + + # Energy grid index at which data starts + threshold_idx = int(ace.xss[idx]) - 1 + + # Get photon production cross section + n_energy = int(ace.xss[idx + 1]) + photon._xs = ace.xss[idx + 2:idx + 2 + n_energy] + + # Determine yield based on ratio of cross sections + energy = ace.xss[ace.jxs[1] + threshold_idx: + ace.jxs[1] + threshold_idx + n_energy] + photon.yield_ = Tabulated1D(energy, photon._xs) + + else: + raise ValueError("MFTYPE must be 12, 13, 16. Got {0}".format( + mftype)) + + # ================================================================== + # Photon energy distribution + + location_start = int(ace.xss[ace.jxs[18] + i]) + distribution = AngleEnergy.from_ace(ace, ace.jxs[19], location_start) + assert isinstance(distribution, UncorrelatedAngleEnergy) + + # ================================================================== + # Photon angular distribution + loc = int(ace.xss[ace.jxs[16] + i]) + + if loc == 0: + # No angular distribution data are given for this reaction, + # isotropic scattering is asssumed in LAB + energy = np.array([photon.yield_.x[0], photon.yield_.x[-1]]) + mu_isotropic = Uniform(-1., 1.) + distribution.angle = AngleDistribution( + energy, [mu_isotropic, mu_isotropic]) + else: + distribution.angle = AngleDistribution.from_ace(ace, ace.jxs[17], loc) + + # Add to list of distributions + photon.distribution.append(distribution) + photons.append(photon) + + return photons + + +class Reaction(object): + """A nuclear reaction + + A Reaction object represents a single reaction channel for a nuclide with + an associated cross section and, if present, a secondary angle and energy + distribution. + + Parameters + ---------- + mt : int + The ENDF MT number for this reaction. On occasion, MCNP uses MT numbers + that don't correspond exactly to the ENDF specification. + + Attributes + ---------- + center_of_mass : bool + Indicates whether scattering kinematics should be performed in the + center-of-mass or laboratory reference frame. + grid above the threshold value in barns. + mt : int + The ENDF MT number for this reaction. + q_value : float + The Q-value of this reaction in MeV. + table : openmc.data.ace.Table + The ACE table which contains this reaction. + threshold : float + Threshold of the reaction in MeV + threshold_idx : int + The index on the energy grid corresponding to the threshold of this + reaction. + xs : openmc.data.Tabulated1D + Microscopic cross section for this reaction as a function of incident + energy + products : Iterable of openmc.data.Product + Reaction products + derived_products : Iterable of openmc.data.Product + Derived reaction products. Used for 'total' fission neutron data when + prompt/delayed data also exists. + + """ + + def __init__(self, mt): + self.center_of_mass = True + self.mt = mt + self.q_value = 0. + self.threshold_idx = 0 + self._xs = None + self.products = [] + self.derived_products = [] + + def __repr__(self): + if self.mt in reaction_name: + return "".format(self.mt, reaction_name[self.mt]) + else: + return "".format(self.mt) + + @property + def center_of_mass(self): + return self._center_of_mass + + @property + def q_value(self): + return self._q_value + + @property + def products(self): + return self._products + + @property + def threshold(self): + return self.xs.x[0] + + @property + def xs(self): + return self._xs + + @center_of_mass.setter + def center_of_mass(self, center_of_mass): + cv.check_type('center of mass', center_of_mass, (bool, np.bool_)) + self._center_of_mass = center_of_mass + + @q_value.setter + def q_value(self, q_value): + cv.check_type('Q value', q_value, Real) + self._q_value = q_value + + @products.setter + def products(self, products): + cv.check_type('reaction products', products, Iterable, Product) + self._products = products + + @xs.setter + def xs(self, xs): + cv.check_type('reaction cross section', xs, Tabulated1D) + for y in xs.y: + cv.check_greater_than('reaction cross section', y, 0.0, True) + self._xs = xs + + def to_hdf5(self, group): + """Write reaction to an HDF5 group + + Parameters + ---------- + group : h5py.Group + HDF5 group to write to + + """ + + group.attrs['mt'] = self.mt + if self.mt in reaction_name: + group.attrs['label'] = np.string_(reaction_name[self.mt]) + else: + group.attrs['label'] = np.string_(self.mt) + group.attrs['Q_value'] = self.q_value + group.attrs['threshold_idx'] = self.threshold_idx + 1 + group.attrs['center_of_mass'] = 1 if self.center_of_mass else 0 + group.attrs['n_product'] = len(self.products) + if self.xs is not None: + group.create_dataset('xs', data=self.xs.y) + for i, p in enumerate(self.products): + pgroup = group.create_group('product_{}'.format(i)) + p.to_hdf5(pgroup) + + @classmethod + def from_hdf5(cls, group, energy): + """Generate reaction from an HDF5 group + + Parameters + ---------- + group : h5py.Group + HDF5 group to write to + energy : Iterable of float + Array of energies at which cross sections are tabulated at + + Returns + ------- + openmc.data.ace.Reaction + Reaction data + + """ + mt = group.attrs['mt'] + rx = cls(mt) + rx.q_value = group.attrs['Q_value'] + rx.threshold_idx = group.attrs['threshold_idx'] - 1 + rx.center_of_mass = bool(group.attrs['center_of_mass']) + + # Read cross section + if 'xs' in group: + xs = group['xs'].value + rx.xs = Tabulated1D(energy, xs) + + # Read reaction products + n_product = group.attrs['n_product'] + products = [] + for i in range(n_product): + pgroup = group['product_{}'.format(i)] + products.append(Product.from_hdf5(pgroup)) + rx.products = products + + return rx + + @classmethod + def from_ace(cls, ace, i_reaction): + # Get nuclide energy grid + n_grid = ace.nxs[3] + grid = ace.xss[ace.jxs[1]:ace.jxs[1] + n_grid] + + if i_reaction > 0: + mt = int(ace.xss[ace.jxs[3] + i_reaction - 1]) + rx = cls(mt) + + # Get Q-value of reaction + rx.q_value = ace.xss[ace.jxs[4] + i_reaction - 1] + + # ================================================================== + # CROSS SECTION + + # Get locator for cross-section data + loc = int(ace.xss[ace.jxs[6] + i_reaction - 1]) + + # Determine starting index on energy grid + rx.threshold_idx = int(ace.xss[ace.jxs[7] + loc - 1]) - 1 + + # Determine number of energies in reaction + n_energy = int(ace.xss[ace.jxs[7] + loc]) + energy = grid[rx.threshold_idx:rx.threshold_idx + n_energy] + + # Read reaction cross section + xs = ace.xss[ace.jxs[7] + loc + 1:ace.jxs[7] + loc + 1 + n_energy] + rx.xs = Tabulated1D(energy, xs) + + # ================================================================== + # YIELD AND ANGLE-ENERGY DISTRIBUTION + + # Determine multiplicity + ty = ace.xss[ace.jxs[5] + i_reaction - 1] + rx.center_of_mass = (ty < 0) + if i_reaction < ace.nxs[5] + 1: + if ty != 19: + if abs(ty) > 100: + # Energy-dependent neutron yield + idx = ace.jxs[11] + abs(ty) - 101 + yield_ = Tabulated1D.from_ace(ace, idx) + else: + yield_ = abs(ty) + + neutron = Product('neutron') + neutron.yield_ = yield_ + rx.products.append(neutron) + else: + assert mt in (18, 19, 20, 21, 38) + rx.products, rx.derived_products = _get_fission_products(ace) + + for p in rx.products: + if p.emission_mode in ('prompt', 'total'): + neutron = p + break + else: + raise Exception("Couldn't find prompt/total fission neutron") + + # Determine locator for ith energy distribution + lnw = int(ace.xss[ace.jxs[10] + i_reaction - 1]) + while lnw > 0: + # Applicability of this distribution + neutron.applicability.append(Tabulated1D.from_ace( + ace, ace.jxs[11] + lnw + 2)) + + # Read energy distribution data + neutron.distribution.append(AngleEnergy.from_ace( + ace, ace.jxs[11], lnw, rx)) + + lnw = int(ace.xss[ace.jxs[11] + lnw - 1]) + + else: + # Elastic scattering + mt = 2 + rx = cls(mt) + + elastic_xs = ace.xss[ace.jxs[1] + 3*n_grid:ace.jxs[1] + 4*n_grid] + rx.xs = Tabulated1D(grid, elastic_xs) + + # No energy distribution for elastic scattering + neutron = Product('neutron') + neutron.distribution.append(UncorrelatedAngleEnergy()) + rx.products.append(neutron) + + # ====================================================================== + # ANGLE DISTRIBUTION (FOR UNCORRELATED) + + if i_reaction < ace.nxs[5] + 1: + # Check if angular distribution data exist + loc = int(ace.xss[ace.jxs[8] + i_reaction]) + if loc <= 0: + # Angular distribution is either given as part of a product + # angle-energy distribution or is not given at all (in which + # case isotropic scattering is assumed) + angle_dist = None + else: + angle_dist = AngleDistribution.from_ace(ace, ace.jxs[9], loc) + + # Apply angular distribution to each uncorrelated angle-energy + # distribution + if angle_dist is not None: + for d in neutron.distribution: + d.angle = angle_dist + + # ====================================================================== + # PHOTON PRODUCTION + + rx.products += _get_photon_products(ace, mt) + + return rx diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index 663b6472c..63e573a01 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -1,9 +1,39 @@ from collections import Iterable +from difflib import get_close_matches from numbers import Real +from warnings import warn import numpy as np +import h5py import openmc.checkvalue as cv +from .ace import Table, get_table +from .container import Tabulated1D + + +_THERMAL_NAMES = {'al': 'c_Al27', 'al27': 'c_Al27', + 'be': 'c_Be', + 'bebeo': 'c_Be_in_BeO', 'be-o': 'c_Be_in_BeO', + 'benz': 'c_Benzine', + 'cah': 'c_Ca_in_CaH2', + 'dd2o': 'c_D_in_D2O', 'hwtr': 'c_D_in_D2O', + 'fe': 'c_Fe56', 'fe56': 'c_Fe56', + 'graph': 'c_Graphite', 'grph': 'c_Graphite', + 'hca': 'c_H_in_CaH2', + 'hch2': 'c_H_in_CH2', 'poly': 'c_H_in_CH2', + 'hh2o': 'c_H_in_H2O', 'lwtr': 'c_H_in_H2O', + 'hzrh': 'c_H_in_ZrH', 'h-zr': 'c_H_in_ZrH', + 'lch4': 'c_liquid_CH4', 'lmeth': 'c_liquid_CH4', + 'mg': 'c_Mg24', + 'obeo': 'c_O_in_BeO', 'o-be': 'c_O_in_BeO', + 'orthod': 'c_ortho_D', 'dortho': 'c_ortho_D', + 'orthoh': 'c_ortho_H', 'hortho': 'c_ortho_H', + 'ouo2': 'c_O_in_UO2', 'o2-u': 'c_O_in_UO2', + 'parad': 'c_para_D', 'dpara': 'c_para_D', + 'parah': 'c_para_H', 'hpara': 'c_para_H', + 'sch4': 'c_solid_CH4', 'smeth': 'c_solid_CH4', + 'uuo2': 'c_U_in_UO2', 'u-o2': 'c_U_in_UO2', + 'zrzrh': 'c_Zr_in_ZrH', 'zr-h': 'c_Zr_in_ZrH'} class CoherentElastic(object): @@ -87,6 +117,277 @@ class CoherentElastic(object): Coherent elastic scattering cross section """ - bragg_edges = dataset.value[0,:] - factors = dataset.value[1,:] + bragg_edges = dataset.value[0, :] + factors = dataset.value[1, :] return cls(bragg_edges, factors) + + +class ThermalScattering(object): + """A ThermalScattering object contains thermal scattering data as represented by + an S(alpha, beta) table. + + Parameters + ---------- + name : str + ZAID identifier of the table, e.g. lwtr.10t. + atomic_weight_ratio : float + Atomic mass ratio of the target nuclide. + temperature : float + Temperature of the target nuclide in eV. + + Attributes + ---------- + atomic_weight_ratio : float + Atomic mass ratio of the target nuclide. + elastic_xs : openmc.data.Tabulated1D or openmc.data.CoherentElastic + Elastic scattering cross section derived in the coherent or incoherent + approximation + inelastic_xs : openmc.data.Tabulated1D + Inelastic scattering cross section derived in the incoherent + approximation + name : str + Name of the table, e.g. lwtr.20t. + temperature : float + Temperature of the target nuclide in eV. + zaids : Iterable of int + ZAID identifiers that the thermal scattering data applies to + + """ + + def __init__(self, name, atomic_weight_ratio, temperature): + self.name = name + self.atomic_weight_ratio = atomic_weight_ratio + self.temperature = temperature + self.elastic_xs = None + self.elastic_mu_out = None + self.inelastic_xs = None + self.inelastic_e_out = None + self.inelastic_mu_out = None + self.secondary_mode = None + self.zaids = [] + + def __repr__(self): + if hasattr(self, 'name'): + return "".format(self.name) + else: + return "" + + def export_to_hdf5(self, path, mode='a'): + """Export table to an HDF5 file. + + Parameters + ---------- + path : str + Path to write HDF5 file to + mode : {'r', r+', 'w', 'x', 'a'} + Mode that is used to open the HDF5 file. This is the second argument + to the :class:`h5py.File` constructor. + + """ + + f = h5py.File(path, mode) + + # Write basic data + g = f.create_group(self.name) + g.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio + g.attrs['temperature'] = self.temperature + g.attrs['zaids'] = self.zaids + + # Write thermal elastic scattering + if self.elastic_xs is not None: + elastic_group = g.create_group('elastic') + self.elastic_xs.to_hdf5(elastic_group, 'xs') + if self.elastic_mu_out is not None: + elastic_group.create_dataset('mu_out', data=self.elastic_mu_out) + + # Write thermal inelastic scattering + if self.inelastic_xs is not None: + inelastic_group = g.create_group('inelastic') + self.inelastic_xs.to_hdf5(inelastic_group, 'xs') + inelastic_group.attrs['secondary_mode'] = np.string_(self.secondary_mode) + if self.secondary_mode in ('equal', 'skewed'): + inelastic_group.create_dataset('energy_out', data=self.inelastic_e_out) + inelastic_group.create_dataset('mu_out', data=self.inelastic_mu_out) + elif self.secondary_mode == 'continuous': + self.inelastic_dist.to_hdf5(inelastic_group) + + @classmethod + def from_hdf5(self, group): + """Generate thermal scattering data from HDF5 group + + Parameters + ---------- + group : h5py.Group + HDF5 group to read from + + Returns + ------- + openmc.data.ThermalScattering + Neutron thermal scattering data + + """ + name = group.name[1:] + atomic_weight_ratio = group.attrs['atomic_weight_ratio'] + temperature = group.attrs['temperature'] + table = ThermalScattering(name, atomic_weight_ratio, temperature) + table.zaids = group.attrs['zaids'] + + # Read thermal elastic scattering + if 'elastic' in group: + elastic_group = group['elastic'] + + # Cross section + elastic_xs_type = elastic_group['xs'].attrs['type'].decode() + if elastic_xs_type == 'tab1': + table.elastic_xs = Tabulated1D.from_hdf5(elastic_group['xs']) + elif elastic_xs_type == 'bragg': + table.elastic_xs = CoherentElastic.from_hdf5(elastic_group['xs']) + + # Angular distribution + if 'mu_out' in elastic_group: + table.elastic_mu_out = elastic_group['mu_out'].value + + # Read thermal inelastic scattering + if 'inelastic' in group: + inelastic_group = group['inelastic'] + table.secondary_mode = inelastic_group.attrs['secondary_mode'].decode() + table.inelastic_xs = Tabulated1D.from_hdf5(inelastic_group['xs']) + if table.secondary_mode in ('equal', 'skewed'): + table.inelastic_e_out = inelastic_group['energy_out'] + table.inelastic_mu_out = inelastic_group['mu_out'] + elif table.secondary_mode == 'continuous': + table.inelastic_dist = AngleEnergy.from_hdf5(inelastic_group) + + return table + + @classmethod + def from_ace(cls, ace_or_filename, name=None): + """Generate thermal scattering data from an ACE table + + Parameters + ---------- + ace : 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. + name : str + GND-conforming name of the material, e.g. c_H_in_H2O. If none is + passed, the appropriate name is guessed based on the name of the ACE + table. + + Returns + ------- + openmc.data.ThermalScattering + Thermal scattering data + + """ + if isinstance(ace_or_filename, Table): + ace = ace_or_filename + else: + ace = get_table(ace_or_filename) + + # Get new name that is GND-consistent + ace_name, xs = ace.name.split('.') + if name is None: + if ace_name.lower() in _THERMAL_NAMES: + name = _THERMAL_NAMES[ace_name.lower()] + '.' + xs + else: + # Make an educated guess?? This actually works well for JEFF-3.2 + # which stupidly uses names like lw00.32t, lw01.32t, etc. for + # different temperatures + matches = get_close_matches( + ace_name.lower(), _THERMAL_NAMES.keys(), cutoff=0.5) + if len(matches) > 0: + name = _THERMAL_NAMES[matches[0]] + '.' + xs + else: + # OK, we give up. Just use the ACE name. + name = 'c_' + ace.name + warn('Thermal scattering material "{}" is not recognized. ' + 'Assigning a name of {}.'.format(ace.name, name)) + + table = cls(name, ace.atomic_weight_ratio, ace.temperature) + + # Incoherent inelastic scattering cross section + idx = ace.jxs[1] + n_energy = int(ace.xss[idx]) + energy = ace.xss[idx+1 : idx+1+n_energy] + xs = ace.xss[idx+1+n_energy : idx+1+2*n_energy] + table.inelastic_xs = Tabulated1D(energy, xs) + + if ace.nxs[7] == 0: + table.secondary_mode = 'equal' + elif ace.nxs[7] == 1: + table.secondary_mode = 'skewed' + elif ace.nxs[7] == 2: + table.secondary_mode = 'continuous' + + n_energy_out = ace.nxs[4] + if table.secondary_mode in ('equal', 'skewed'): + n_mu = ace.nxs[3] + idx = ace.jxs[3] + table.inelastic_e_out = ace.xss[idx:idx+n_energy*n_energy_out*(n_mu+2):n_mu+2] + table.inelastic_e_out.shape = (n_energy, n_energy_out) + + table.inelastic_mu_out = ace.xss[idx:idx+n_energy*n_energy_out*(n_mu+2)] + table.inelastic_mu_out.shape = (n_energy, n_energy_out, n_mu+2) + table.inelastic_mu_out = table.inelastic_mu_out[:, :, 1:] + else: + n_mu = ace.nxs[3] - 1 + idx = ace.jxs[3] + locc = ace.xss[idx:idx + n_energy].astype(int) + n_energy_out = ace.xss[idx + n_energy:idx + 2*n_energy].astype(int) + energy_out = [] + mu_out = [] + for i in range(n_energy): + idx = locc[i] + + # Outgoing energy distribution for incoming energy i + e = ace.xss[idx + 1:idx + 1 + n_energy_out[i]*(n_mu + 3):n_mu + 3] + p = ace.xss[idx + 2:idx + 2 + n_energy_out[i]*(n_mu + 3):n_mu + 3] + c = ace.xss[idx + 3:idx + 3 + n_energy_out[i]*(n_mu + 3):n_mu + 3] + eout_i = Tabular(e, p, 'linear-linear', ignore_negative=True) + eout_i.c = c + + # Outgoing angle distribution for each (incoming, outgoing) energy pair + mu_i = [] + for j in range(n_energy_out[i]): + mu = ace.xss[idx + 4:idx + 4 + n_mu] + p_mu = 1./n_mu*np.ones(n_mu) + mu_ij = Discrete(mu, p_mu) + mu_ij.c = np.cumsum(p_mu) + mu_i.append(mu_ij) + idx += 3 + n_mu + + energy_out.append(eout_i) + mu_out.append(mu_i) + + # Create correlated angle-energy distribution + breakpoints = [n_energy] + interpolation = [2] + energy = inelastic_xs.x + table.inelastic_dist = CorrelatedAngleEnergy( + breakpoints, interpolation, energy, energy_out, mu_out) + + # Incoherent/coherent elastic scattering cross section + idx = ace.jxs[4] + if idx != 0: + n_energy = int(ace.xss[idx]) + energy = ace.xss[idx+1 : idx+1+n_energy] + P = ace.xss[idx+1+n_energy : idx+1+2*n_energy] + + if ace.nxs[5] == 4: + table.elastic_xs = CoherentElastic(energy, P) + else: + table.elastic_xs = Tabulated1D(energy, P) + + # Angular distribution + n_mu = ace.nxs[6] + if n_mu != -1: + idx = ace.jxs[6] + table.elastic_mu_out = ace.xss[idx:idx + n_energy*n_mu] + table.elastic_mu_out.shape = (n_energy, n_mu) + + # Get relevant ZAIDs + pairs = np.fromiter(map(lambda p: p[0], ace.pairs), int) + table.zaids = pairs[np.nonzero(pairs)] + + return table diff --git a/openmc/data/urr.py b/openmc/data/urr.py index 308a45917..b010a6c47 100644 --- a/openmc/data/urr.py +++ b/openmc/data/urr.py @@ -169,3 +169,42 @@ class ProbabilityTables(object): return cls(energy, table, interpolation, inelastic_flag, absorption_flag, multiply_smooth) + + @classmethod + def from_ace(cls, ace): + """Generate probability tables from an ACE table + + Parameters + ---------- + ace : openmc.data.ace.Table + ACE table to read from + + Returns + ------- + openmc.data.ProbabilityTables + Unresolved resonance region probability tables + + """ + # Check if URR probability tables are present + idx = ace.jxs[23] + if idx == 0: + return None + + N = int(ace.xss[idx]) # Number of incident energies + M = int(ace.xss[idx+1]) # Length of probability table + interpolation = int(ace.xss[idx+2]) + inelastic_flag = int(ace.xss[idx+3]) + absorption_flag = int(ace.xss[idx+4]) + multiply_smooth = (int(ace.xss[idx+5]) == 1) + idx += 6 + + # Get energies at which tables exist + energy = ace.xss[idx : idx+N] + idx += N + + # Get probability tables + table = ace.xss[idx : idx+N*6*M] + table.shape = (N, 6, M) + + return cls(energy, table, interpolation, inelastic_flag, + absorption_flag, multiply_smooth) diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index 9df9d150d..179b47330 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -625,7 +625,7 @@ class Library(object): in the report. Defaults to 'all'. nuclides : {'all', 'sum'} The nuclides of the cross-sections to include in the report. This - may be a list of nuclide name strings (e.g., ['U-235', 'U-238']). + may be a list of nuclide name strings (e.g., ['U235', 'U238']). The special string 'all' will report the cross sections for all nuclides in the spatial domain. The special string 'sum' will report the cross sections summed over all nuclides. Defaults to 'all'. @@ -758,7 +758,7 @@ class Library(object): xsdata_name : str Name to apply to the "xsdata" entry produced by this method nuclide : str - A nuclide name string (e.g., 'U-235'). Defaults to 'total' to + A nuclide name string (e.g., 'U235'). Defaults to 'total' to obtain a material-wise macroscopic cross section. xs_type: {'macro', 'micro'} Provide the macro or micro cross section in units of cm^-1 or diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index 2c0cb21c6..8e6d123f7 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -120,7 +120,7 @@ class MGXS(object): being tracked. This is unity if the by_nuclide attribute is False. nuclides : Iterable of str or 'sum' The optional user-specified nuclides for which to compute cross - sections (e.g., 'U-238', 'O-16'). If by_nuclide is True but nuclides + sections (e.g., 'U238', 'O16'). If by_nuclide is True but nuclides are not specified by the user, all nuclides in the spatial domain are included. This attribute is 'sum' if by_nuclide is false. sparse : bool @@ -494,7 +494,7 @@ class MGXS(object): ------- list of str A list of the string names for each nuclide in the spatial domain - (e.g., ['U-235', 'U-238', 'O-16']) + (e.g., ['U235', 'U238', 'O16']) Raises ------ @@ -522,7 +522,7 @@ class MGXS(object): Parameters ---------- nuclide : str - A nuclide name string (e.g., 'U-235') + A nuclide name string (e.g., 'U235') Returns ------- @@ -557,7 +557,7 @@ class MGXS(object): Parameters ---------- nuclides : Iterable of str or 'all' or 'sum' - A list of nuclide name strings (e.g., ['U-235', 'U-238']). The + A list of nuclide name strings (e.g., ['U235', 'U238']). The special string 'all' will return the atom densities for all nuclides in the spatial domain. The special string 'sum' will return the atom density summed across all nuclides in the spatial domain. Defaults @@ -716,7 +716,7 @@ class MGXS(object): subdomains : Iterable of Integral or 'all' Subdomain IDs of interest. Defaults to 'all'. nuclides : Iterable of str or 'all' or 'sum' - A list of nuclide name strings (e.g., ['U-235', 'U-238']). The + A list of nuclide name strings (e.g., ['U235', 'U238']). The special string 'all' will return the cross sections for all nuclides in the spatial domain. The special string 'sum' will return the cross section summed over all nuclides. Defaults to 'all'. @@ -954,7 +954,7 @@ class MGXS(object): ---------- nuclides : list of str A list of nuclide name strings - (e.g., ['U-235', 'U-238']; default is []) + (e.g., ['U235', 'U238']; default is []) groups : list of int A list of energy group indices starting at 1 for the high energies (e.g., [1, 2, 3]; default is []) @@ -1111,7 +1111,7 @@ class MGXS(object): Defaults to 'all'. nuclides : Iterable of str or 'all' or 'sum' The nuclides of the cross-sections to include in the report. This - may be a list of nuclide name strings (e.g., ['U-235', 'U-238']). + may be a list of nuclide name strings (e.g., ['U235', 'U238']). The special string 'all' will report the cross sections for all nuclides in the spatial domain. The special string 'sum' will report the cross sections summed over all nuclides. Defaults to 'all'. @@ -1215,7 +1215,7 @@ class MGXS(object): Defaults to 'all'. nuclides : Iterable of str or 'all' or 'sum' The nuclides of the cross-sections to include in the report. This - may be a list of nuclide name strings (e.g., ['U-235', 'U-238']). + may be a list of nuclide name strings (e.g., ['U235', 'U238']). The special string 'all' will report the cross sections for all nuclides in the spatial domain. The special string 'sum' will report the cross sections summed over all nuclides. Defaults to 'all'. @@ -1414,7 +1414,7 @@ class MGXS(object): Energy groups of interest. Defaults to 'all'. nuclides : Iterable of str or 'all' or 'sum' The nuclides of the cross-sections to include in the dataframe. This - may be a list of nuclide name strings (e.g., ['U-235', 'U-238']). + may be a list of nuclide name strings (e.g., ['U235', 'U238']). The special string 'all' will include the cross sections for all nuclides in the spatial domain. The special string 'sum' will include the cross sections summed over all nuclides. Defaults @@ -1603,7 +1603,7 @@ class MatrixMGXS(MGXS): being tracked. This is unity if the by_nuclide attribute is False. nuclides : Iterable of str or 'sum' The optional user-specified nuclides for which to compute cross - sections (e.g., 'U-238', 'O-16'). If by_nuclide is True but nuclides + sections (e.g., 'U238', 'O16'). If by_nuclide is True but nuclides are not specified by the user, all nuclides in the spatial domain are included. This attribute is 'sum' if by_nuclide is false. sparse : bool @@ -1652,7 +1652,7 @@ class MatrixMGXS(MGXS): subdomains : Iterable of Integral or 'all' Subdomain IDs of interest. Defaults to 'all'. nuclides : Iterable of str or 'all' or 'sum' - A list of nuclide name strings (e.g., ['U-235', 'U-238']). The + A list of nuclide name strings (e.g., ['U235', 'U238']). The special string 'all' will return the cross sections for all nuclides in the spatial domain. The special string 'sum' will return the cross section summed over all nuclides. Defaults to @@ -1790,7 +1790,7 @@ class MatrixMGXS(MGXS): ---------- nuclides : list of str A list of nuclide name strings - (e.g., ['U-235', 'U-238']; default is []) + (e.g., ['U235', 'U238']; default is []) in_groups : list of int A list of incoming energy group indices starting at 1 for the high energies (e.g., [1, 2, 3]; default is []) @@ -1840,7 +1840,7 @@ class MatrixMGXS(MGXS): Defaults to 'all'. nuclides : Iterable of str or 'all' or 'sum' The nuclides of the cross-sections to include in the report. This - may be a list of nuclide name strings (e.g., ['U-235', 'U-238']). + may be a list of nuclide name strings (e.g., ['U235', 'U238']). The special string 'all' will report the cross sections for all nuclides in the spatial domain. The special string 'sum' will report the cross sections summed over all nuclides. Defaults to @@ -2024,7 +2024,7 @@ class TotalXS(MGXS): being tracked. This is unity if the by_nuclide attribute is False. nuclides : Iterable of str or 'sum' The optional user-specified nuclides for which to compute cross - sections (e.g., 'U-238', 'O-16'). If by_nuclide is True but nuclides + sections (e.g., 'U238', 'O16'). If by_nuclide is True but nuclides are not specified by the user, all nuclides in the spatial domain are included. This attribute is 'sum' if by_nuclide is false. sparse : bool @@ -2142,7 +2142,7 @@ class TransportXS(MGXS): being tracked. This is unity if the by_nuclide attribute is False. nuclides : Iterable of str or 'sum' The optional user-specified nuclides for which to compute cross - sections (e.g., 'U-238', 'O-16'). If by_nuclide is True but nuclides + sections (e.g., 'U238', 'O16'). If by_nuclide is True but nuclides are not specified by the user, all nuclides in the spatial domain are included. This attribute is 'sum' if by_nuclide is false. sparse : bool @@ -2272,7 +2272,7 @@ class NuTransportXS(TransportXS): being tracked. This is unity if the by_nuclide attribute is False. nuclides : Iterable of str or 'sum' The optional user-specified nuclides for which to compute cross - sections (e.g., 'U-238', 'O-16'). If by_nuclide is True but nuclides + sections (e.g., 'U238', 'O16'). If by_nuclide is True but nuclides are not specified by the user, all nuclides in the spatial domain are included. This attribute is 'sum' if by_nuclide is false. sparse : bool @@ -2393,7 +2393,7 @@ class AbsorptionXS(MGXS): being tracked. This is unity if the by_nuclide attribute is False. nuclides : Iterable of str or 'sum' The optional user-specified nuclides for which to compute cross - sections (e.g., 'U-238', 'O-16'). If by_nuclide is True but nuclides + sections (e.g., 'U238', 'O16'). If by_nuclide is True but nuclides are not specified by the user, all nuclides in the spatial domain are included. This attribute is 'sum' if by_nuclide is false. sparse : bool @@ -2509,7 +2509,7 @@ class CaptureXS(MGXS): being tracked. This is unity if the by_nuclide attribute is False. nuclides : Iterable of str or 'sum' The optional user-specified nuclides for which to compute cross - sections (e.g., 'U-238', 'O-16'). If by_nuclide is True but nuclides + sections (e.g., 'U238', 'O16'). If by_nuclide is True but nuclides are not specified by the user, all nuclides in the spatial domain are included. This attribute is 'sum' if by_nuclide is false. sparse : bool @@ -2631,7 +2631,7 @@ class FissionXS(MGXS): being tracked. This is unity if the by_nuclide attribute is False. nuclides : Iterable of str or 'sum' The optional user-specified nuclides for which to compute cross - sections (e.g., 'U-238', 'O-16'). If by_nuclide is True but nuclides + sections (e.g., 'U238', 'O16'). If by_nuclide is True but nuclides are not specified by the user, all nuclides in the spatial domain are included. This attribute is 'sum' if by_nuclide is false. sparse : bool @@ -2742,7 +2742,7 @@ class NuFissionXS(MGXS): being tracked. This is unity if the by_nuclide attribute is False. nuclides : Iterable of str or 'sum' The optional user-specified nuclides for which to compute cross - sections (e.g., 'U-238', 'O-16'). If by_nuclide is True but nuclides + sections (e.g., 'U238', 'O16'). If by_nuclide is True but nuclides are not specified by the user, all nuclides in the spatial domain are included. This attribute is 'sum' if by_nuclide is false. sparse : bool @@ -2858,7 +2858,7 @@ class KappaFissionXS(MGXS): being tracked. This is unity if the by_nuclide attribute is False. nuclides : Iterable of str or 'sum' The optional user-specified nuclides for which to compute cross - sections (e.g., 'U-238', 'O-16'). If by_nuclide is True but nuclides + sections (e.g., 'U238', 'O16'). If by_nuclide is True but nuclides are not specified by the user, all nuclides in the spatial domain are included. This attribute is 'sum' if by_nuclide is false. sparse : bool @@ -2971,7 +2971,7 @@ class ScatterXS(MGXS): being tracked. This is unity if the by_nuclide attribute is False. nuclides : Iterable of str or 'sum' The optional user-specified nuclides for which to compute cross - sections (e.g., 'U-238', 'O-16'). If by_nuclide is True but nuclides + sections (e.g., 'U238', 'O16'). If by_nuclide is True but nuclides are not specified by the user, all nuclides in the spatial domain are included. This attribute is 'sum' if by_nuclide is false. sparse : bool @@ -3086,7 +3086,7 @@ class NuScatterXS(MGXS): being tracked. This is unity if the by_nuclide attribute is False. nuclides : Iterable of str or 'sum' The optional user-specified nuclides for which to compute cross - sections (e.g., 'U-238', 'O-16'). If by_nuclide is True but nuclides + sections (e.g., 'U238', 'O16'). If by_nuclide is True but nuclides are not specified by the user, all nuclides in the spatial domain are included. This attribute is 'sum' if by_nuclide is false. sparse : bool @@ -3220,7 +3220,7 @@ class ScatterMatrixXS(MatrixMGXS): being tracked. This is unity if the by_nuclide attribute is False. nuclides : Iterable of str or 'sum' The optional user-specified nuclides for which to compute cross - sections (e.g., 'U-238', 'O-16'). If by_nuclide is True but nuclides + sections (e.g., 'U238', 'O16'). If by_nuclide is True but nuclides are not specified by the user, all nuclides in the spatial domain are included. This attribute is 'sum' if by_nuclide is false. sparse : bool @@ -3383,7 +3383,7 @@ class ScatterMatrixXS(MatrixMGXS): ---------- nuclides : list of str A list of nuclide name strings - (e.g., ['U-235', 'U-238']; default is []) + (e.g., ['U235', 'U238']; default is []) in_groups : list of int A list of incoming energy group indices starting at 1 for the high energies (e.g., [1, 2, 3]; default is []) @@ -3465,7 +3465,7 @@ class ScatterMatrixXS(MatrixMGXS): subdomains : Iterable of Integral or 'all' Subdomain IDs of interest. Defaults to 'all'. nuclides : Iterable of str or 'all' or 'sum' - A list of nuclide name strings (e.g., ['U-235', 'U-238']). The + A list of nuclide name strings (e.g., ['U235', 'U238']). The special string 'all' will return the cross sections for all nuclides in the spatial domain. The special string 'sum' will return the cross section summed over all nuclides. Defaults to 'all'. @@ -3612,7 +3612,7 @@ class ScatterMatrixXS(MatrixMGXS): Energy groups of interest. Defaults to 'all'. nuclides : Iterable of str or 'all' or 'sum' The nuclides of the cross-sections to include in the dataframe. This - may be a list of nuclide name strings (e.g., ['U-235', 'U-238']). + may be a list of nuclide name strings (e.g., ['U235', 'U238']). The special string 'all' will include the cross sections for all nuclides in the spatial domain. The special string 'sum' will include the cross sections summed over all nuclides. Defaults @@ -3679,7 +3679,7 @@ class ScatterMatrixXS(MatrixMGXS): Defaults to 'all'. nuclides : Iterable of str or 'all' or 'sum' The nuclides of the cross-sections to include in the report. This - may be a list of nuclide name strings (e.g., ['U-235', 'U-238']). + may be a list of nuclide name strings (e.g., ['U235', 'U238']). The special string 'all' will report the cross sections for all nuclides in the spatial domain. The special string 'sum' will report the cross sections summed over all nuclides. Defaults to 'all'. @@ -3869,7 +3869,7 @@ class NuScatterMatrixXS(ScatterMatrixXS): being tracked. This is unity if the by_nuclide attribute is False. nuclides : Iterable of str or 'sum' The optional user-specified nuclides for which to compute cross - sections (e.g., 'U-238', 'O-16'). If by_nuclide is True but nuclides + sections (e.g., 'U238', 'O16'). If by_nuclide is True but nuclides are not specified by the user, all nuclides in the spatial domain are included. This attribute is 'sum' if by_nuclide is false. sparse : bool @@ -3991,7 +3991,7 @@ class MultiplicityMatrixXS(MatrixMGXS): being tracked. This is unity if the by_nuclide attribute is False. nuclides : Iterable of str or 'sum' The optional user-specified nuclides for which to compute cross - sections (e.g., 'U-238', 'O-16'). If by_nuclide is True but nuclides + sections (e.g., 'U238', 'O16'). If by_nuclide is True but nuclides are not specified by the user, all nuclides in the spatial domain are included. This attribute is 'sum' if by_nuclide is false. sparse : bool @@ -4138,7 +4138,7 @@ class NuFissionMatrixXS(MatrixMGXS): being tracked. This is unity if the by_nuclide attribute is False. nuclides : Iterable of str or 'sum' The optional user-specified nuclides for which to compute cross - sections (e.g., 'U-238', 'O-16'). If by_nuclide is True but nuclides + sections (e.g., 'U238', 'O16'). If by_nuclide is True but nuclides are not specified by the user, all nuclides in the spatial domain are included. This attribute is 'sum' if by_nuclide is false. sparse : bool @@ -4253,7 +4253,7 @@ class Chi(MGXS): being tracked. This is unity if the by_nuclide attribute is False. nuclides : Iterable of str or 'sum' The optional user-specified nuclides for which to compute cross - sections (e.g., 'U-238', 'O-16'). If by_nuclide is True but nuclides + sections (e.g., 'U238', 'O16'). If by_nuclide is True but nuclides are not specified by the user, all nuclides in the spatial domain are included. This attribute is 'sum' if by_nuclide is false. sparse : bool @@ -4331,7 +4331,7 @@ class Chi(MGXS): ---------- nuclides : list of str A list of nuclide name strings - (e.g., ['U-235', 'U-238']; default is []) + (e.g., ['U235', 'U238']; default is []) groups : list of Integral A list of energy group indices starting at 1 for the high energies (e.g., [1, 2, 3]; default is []) @@ -4442,7 +4442,7 @@ class Chi(MGXS): subdomains : Iterable of Integral or 'all' Subdomain IDs of interest. Defaults to 'all'. nuclides : Iterable of str or 'all' or 'sum' - A list of nuclide name strings (e.g., ['U-235', 'U-238']). The + A list of nuclide name strings (e.g., ['U235', 'U238']). The special string 'all' will return the cross sections for all nuclides in the spatial domain. The special string 'sum' will return the cross section summed over all nuclides. Defaults to 'all'. @@ -4575,7 +4575,7 @@ class Chi(MGXS): Energy groups of interest. Defaults to 'all'. nuclides : Iterable of str or 'all' or 'sum' The nuclides of the cross-sections to include in the dataframe. This - may be a list of nuclide name strings (e.g., ['U-235', 'U-238']). + may be a list of nuclide name strings (e.g., ['U235', 'U238']). The special string 'all' will include the cross sections for all nuclides in the spatial domain. The special string 'sum' will include the cross sections summed over all nuclides. Defaults to diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index 8a572c874..6f8188b10 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -100,7 +100,7 @@ class XSdata(object): alias : str Separate unique identifier for the xsdata object zaid : int - 1000*(atomic number) + mass number. As an example, the zaid of U-235 + 1000*(atomic number) + mass number. As an example, the zaid of U235 would be 92235. awr : float Atomic weight ratio of an isotope. That is, the ratio of the mass diff --git a/openmc/nuclide.py b/openmc/nuclide.py index bfeff9311..14609161f 100644 --- a/openmc/nuclide.py +++ b/openmc/nuclide.py @@ -13,18 +13,18 @@ class Nuclide(object): Parameters ---------- name : str - Name of the nuclide, e.g. U-235 + Name of the nuclide, e.g. U235 xs : str Cross section identifier, e.g. 71c Attributes ---------- name : str - Name of the nuclide, e.g. U-235 + Name of the nuclide, e.g. U235 xs : str Cross section identifier, e.g. 71c zaid : int - 1000*(atomic number) + mass number. As an example, the zaid of U-235 + 1000*(atomic number) + mass number. As an example, the zaid of U235 would be 92235. scattering : 'data' or 'iso-in-lab' or None The type of angular scattering distribution to use diff --git a/openmc/tallies.py b/openmc/tallies.py index fab88e2f7..79e7c56bc 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -520,7 +520,7 @@ class Tally(object): Nuclide to add to the tally. The nuclide should be a Nuclide object when a user is adding nuclides to a Tally for input file generation. The nuclide is a str when a Tally is created from a StatePoint file - (e.g., 'H-1', 'U-235') unless a Summary has been linked with the + (e.g., 'H1', 'U235') unless a Summary has been linked with the StatePoint. The nuclide may be a CrossNuclide or AggregateNuclide for derived tallies created by tally arithmetic. @@ -1166,7 +1166,7 @@ class Tally(object): Parameters ---------- nuclide : str - The name of the Nuclide (e.g., 'H-1', 'U-238') + The name of the Nuclide (e.g., 'H1', 'U238') Returns ------- @@ -1342,7 +1342,7 @@ class Tally(object): ---------- nuclides : list of str A list of nuclide name strings - (e.g., ['U-235', 'U-238']; default is []) + (e.g., ['U235', 'U238']; default is []) Returns ------- @@ -1435,7 +1435,7 @@ class Tally(object): the filter_types parameter. nuclides : list of str A list of nuclide name strings - (e.g., ['U-235', 'U-238']; default is []) + (e.g., ['U235', 'U238']; default is []) value : str A string for the type of value to return - 'mean' (default), 'std_dev', 'rel_err', 'sum', or 'sum_sq' are accepted @@ -2887,7 +2887,7 @@ class Tally(object): correspond to the filter_types parameter. nuclides : list of str A list of nuclide name strings - (e.g., ['U-235', 'U-238']; default is []) + (e.g., ['U235', 'U238']; default is []) Returns ------- @@ -3025,7 +3025,7 @@ class Tally(object): interest. nuclides : list of str A list of nuclide name strings to sum across - (e.g., ['U-235', 'U-238']; default is []) + (e.g., ['U235', 'U238']; default is []) remove_filter : bool If a filter is being summed over, this bool indicates whether to remove that filter in the returned tally. Default is False. @@ -3173,7 +3173,7 @@ class Tally(object): interest. nuclides : list of str A list of nuclide name strings to average across - (e.g., ['U-235', 'U-238']; default is []) + (e.g., ['U235', 'U238']; default is []) remove_filter : bool If a filter is being averaged over, this bool indicates whether to remove that filter in the returned tally. Default is False. diff --git a/scripts/openmc-ace-to-hdf5 b/scripts/openmc-ace-to-hdf5 new file mode 100755 index 000000000..743cf4a50 --- /dev/null +++ b/scripts/openmc-ace-to-hdf5 @@ -0,0 +1,139 @@ +#!/usr/bin/env python + +import argparse +import os +import xml.etree.ElementTree as ET +import warnings + +import openmc.data + +description = """ +This script can be used to create HDF5 nuclear data libraries used by +OpenMC. There are four different ways you can specify ACE libraries that are to +be converted: + +1. List each ACE library as a positional argument. This is very useful in + conjunction with the usual shell utilities (ls, find, etc.). +2. Use the --xml option to specify a pre-v0.9 cross_sections.xml file. +3. Use the --xsdir option to specify a MCNP xsdir file. +4. Use the --xsdata option to specify a Serpent xsdata file. + +The script does not use any extra information from cross_sections.xml/ xsdir/ +xsdata files to determine whether the nuclide is metastable. Instead, the +--metastable argument can be used to specify whether the ZAID naming convention +follows the NNDC data convention (1000*Z + A + 300 + 100*m), or the MCNP data +convention (essentially the same as NNDC, except that the first metastable state +of Am242 is 95242 and the ground state is 95642). + +""" + +class CustomFormatter(argparse.ArgumentDefaultsHelpFormatter, + argparse.RawDescriptionHelpFormatter): + pass + +parser = argparse.ArgumentParser( + description=description, + formatter_class=CustomFormatter +) +parser.add_argument('libraries', nargs='*', + help='ACE libraries to convert to HDF5') +parser.add_argument('-d', '--destination', default='.', + help='Directory to create new library in') +parser.add_argument('-m', '--metastable', choices=['mcnp', 'nndc'], default='nndc', + help='How to interpret ZAIDs for metastable nuclides') +parser.add_argument('--xml', help='Old-style cross_sections.xml that ' + 'lists ACE libraries') +parser.add_argument('--xsdir', help='MCNP xsdir file that lists ' + 'ACE libraries') +parser.add_argument('--xsdata', help='Serpent xsdata file that lists ' + 'ACE libraries') +args = parser.parse_args() + +if not os.path.isdir(args.destination): + os.mkdir(args.destination) + +# If the --xml argument was given, get the list of ACE libraries directory from +# elements within the specified cross_sections.xml file +ace_libraries = [] +if args.xml is not None: + tree = ET.parse(args.xml) + root = tree.getroot() + if root.find('directory') is not None: + directory = root.find('directory').text + else: + directory = os.path.dirname(args.xml) + + for ace_table in root.findall('ace_table'): + ace_libraries.append(os.path.join(directory, ace_table.attrib['path'])) + +elif args.xsdir is not None: + # Find 'directory' section + lines = open(args.xsdir, 'r').readlines() + for index, line in enumerate(lines): + if line.strip().lower() == 'directory': + break + else: + raise IOError("Could not find 'directory' section in MCNP xsdir file") + + # Create list of ACE libraries + for line in lines[index + 1:]: + words = line.split() + if len(words) < 3: + continue + + path = os.path.join(os.path.dirname(args.xsdir), words[2]) + if path not in ace_libraries: + ace_libraries.append(path) + +elif args.xsdata is not None: + with open(args.xsdata, 'r') as xsdata: + for line in xsdata: + words = line.split() + if len(words) >= 9: + path = os.path.join(os.path.dirname(args.xsdata, words[8])) + if path not in ace_libraries: + ace_libraries.append(path) + +else: + ace_libraries = args.libraries + +library = openmc.data.DataLibrary() + +for filename in ace_libraries: + # Check that ACE library exists + if not os.path.exists(filename): + warnings.warn("ACE library '{}' does not exist.".format(filename)) + continue + + lib = openmc.data.ace.Library(filename) + for table in lib.tables: + if table.name.endswith('c'): + # Continuous-energy neutron data + neutron = openmc.data.IncidentNeutron.from_ace( + table, args.metastable) + print(neutron.name) + + # Determine filename + outfile = os.path.join(args.destination, + neutron.name.replace('.', '_') + '.h5') + neutron.export_to_hdf5(outfile) + + # Register with library + library.register_file(outfile) + + elif table.name.endswith('t'): + # Thermal scattering data + thermal = openmc.data.ThermalScattering.from_ace(table) + print(thermal.name) + + # Determine filename + outfile = os.path.join(args.destination, + thermal.name.replace('.', '_') + '.h5') + thermal.export_to_hdf5(outfile) + + # Register with library + library.register_file(outfile, 'thermal') + +# Write cross_sections.xml +libpath = os.path.join(args.destination, 'cross_sections.xml') +library.export_to_xml(libpath) diff --git a/scripts/openmc-ascii-to-binary b/scripts/openmc-ascii-to-binary deleted file mode 100755 index e707e2a85..000000000 --- a/scripts/openmc-ascii-to-binary +++ /dev/null @@ -1,13 +0,0 @@ -#!/usr/bin/env python - -from openmc.ace import ascii_to_binary -import sys - - -if __name__ == '__main__': - # Check for proper number of arguments - if len(sys.argv) < 3: - sys.exit('Usage: {0} ascii_file binary_file'.format(sys.argv[0])) - - # Convert ASCII file - ascii_to_binary(sys.argv[1], sys.argv[2]) diff --git a/scripts/openmc-update-inputs b/scripts/openmc-update-inputs index 2a8097854..5b5bf0978 100755 --- a/scripts/openmc-update-inputs +++ b/scripts/openmc-update-inputs @@ -22,11 +22,14 @@ optional arguments: from __future__ import print_function import argparse +from difflib import get_close_matches from itertools import chain from random import randint from shutil import move import xml.etree.ElementTree as ET +import openmc.data +from openmc.data.thermal import _THERMAL_NAMES description = "Update OpenMC's input XML files to the latest format." epilog = """\ @@ -40,6 +43,11 @@ geometry.xml: Lattices containing 'outside' attributes/tags will be replaced with lattices containing 'outer' attributes, and the appropriate cells/universes will be added. Any 'surfaces' attributes/elements on a cell will be renamed 'region'. + +materials.xml: Nuclide names will be changed from ACE aliases (e.g., Am-242m) to + HDF5/GND names (e.g., Am242_m1). Thermal scattering table names will be + changed from ACE aliases (e.g., HH2O) to HDF5/GND names (e.g., c_H_in_H2O). + """ @@ -245,6 +253,62 @@ def update_geometry(geometry_root): return was_updated +def get_thermal_name(name): + """Get proper S(a,b) table name, e.g. 'HH2O' -> 'c_H_in_H2O'""" + if name.lower() in _THERMAL_NAMES: + return _THERMAL_NAMES[name.lower()] + else: + # Make an educated guess?? This actually works well for + # JEFF-3.2 which stupidly uses names like lw00.32t, + # lw01.32t, etc. for different temperatures + matches = get_close_matches( + name.lower(), _THERMAL_NAMES.keys(), cutoff=0.5) + if len(matches) > 0: + return _THERMAL_NAMES[matches[0]] + '.' + xs + else: + # OK, we give up. Just use the ACE name. + return 'c_' + name + return name + +def update_materials(root): + """Update the given XML materials tree. Return True if changes were made.""" + was_updated = False + + for material in root.findall('material'): + for nuclide in material.findall('nuclide'): + if 'name' in nuclide.attrib: + nucname = nuclide.attrib['name'] + nucname = nucname.replace('-', '') + nucname = nucname.replace('Nat', '0') + if nucname.endswith('m'): + nucname = nucname[:-1] + '_m1' + nuclide.set('name', nucname) + was_updated = True + + elif nuclide.find('name') is not None: + name_elem = nuclide.find('name') + nucname = name_elem.text + nucname = nucname.replace('-', '') + nucname = nucname.replace('Nat', '0') + if nucname.endswith('m'): + nucname = nucname[:-1] + '_m1' + name_elem.text = nucname + was_updated = True + + for sab in material.findall('sab'): + if 'name' in sab.attrib: + sabname = sab.attrib['name'] + sab.set('name', get_thermal_name(sabname)) + was_updated = True + + elif sab.find('name') is not None: + name_elem = sab.find('name') + sabname = name_elem.text + name_elem.text = get_thermal(sabname) + was_updated = True + + return was_updated + if __name__ == '__main__': args = parse_args() @@ -256,6 +320,8 @@ if __name__ == '__main__': if root.tag == 'geometry': was_updated = update_geometry(root) + elif root.tag == 'materials': + was_updated = update_materials(root) if was_updated: # Move the original geometry file to preserve it. diff --git a/scripts/openmc-xsdata-to-xml b/scripts/openmc-xsdata-to-xml deleted file mode 100755 index 708a196f9..000000000 --- a/scripts/openmc-xsdata-to-xml +++ /dev/null @@ -1,148 +0,0 @@ -#!/usr/bin/env python - -import os -import sys -from xml.dom.minidom import getDOMImplementation - -types = {1: "neutron", 2: "dosimetry", 3: "thermal"} - - -class Xsdata(object): - - def __init__(self, filename): - self._table_dict = {} - self.tables = [] - - for line in open(filename, 'r'): - words = line.split() - - # If this listing is just an alias listing, only assign the alias - # attribute - name = words[1] - alias = words[0] - table = self.find_table(name) - if table: - if name not in table.alias: - table.alias.append(alias) - continue - - table = XsdataTable() - table.name = name - table.type = types[int(words[2])] - table.zaid = int(words[3]) - table.metastable = int(words[4]) - table.awr = float(words[5]) - table.temperature = 8.6173423e-11 * float(words[6]) - table.binary = int(words[7]) - table.path = words[8] - - self.tables.append(table) - self._table_dict[name] = table - - # Check for common directory - self.directory = os.path.dirname(self.tables[0].path) - for table in self.tables: - if not table.path.startswith(self.directory): - self.directory = None - break - - def to_xml(self): - # Create XML document - impl = getDOMImplementation() - doc = impl.createDocument(None, "cross_sections", None) - - # Get root element - root = doc.documentElement - - # Add a directory node - if self.directory: - directoryNode = doc.createElement("directory") - text = doc.createTextNode(self.directory) - directoryNode.appendChild(text) - root.appendChild(directoryNode) - - for table in self.tables: - table.path = os.path.basename(table.path) - - # Add a node for each table - for table in self.tables: - node = table.to_xml_node(doc) - root.appendChild(node) - - return doc - - def find_table(self, name): - if name in self._table_dict: - return self._table_dict[name] - else: - return None - - -class XsdataTable(object): - - def __init__(self): - self.alias = [] - - def to_xml_node(self, doc): - node = doc.createElement("ace_table") - node.setAttribute("name", self.name) - for attribute in ["alias", "zaid", "type", "metastable", - "awr", "temperature", "binary", "path"]: - if hasattr(self, attribute): - # Join string for alias attribute - if attribute == "alias": - if not self.alias: - continue - string = " ".join(self.alias) - else: - string = "{0}".format(getattr(self, attribute)) - - # Skip metastable and binary if 0 - if attribute == "metastable" and self.metastable == 0: - continue - if attribute == "binary" and self.binary == 0: - continue - - # Create attribute node - # nodeAttr = doc.createElement(attribute) - # text = doc.createTextNode(string) - # nodeAttr.appendChild(text) - # node.appendChild(nodeAttr) - node.setAttribute(attribute, string) - return node - - -if __name__ == '__main__': - # Read command line arguments - if len(sys.argv) < 3: - sys.exit("Usage: convert_xsdata.py xsdataFile xmlFile") - xsdataFile = sys.argv[1] - xmlFile = sys.argv[2] - - # Read xsdata and create XML document object - xsdataObject = Xsdata(xsdataFile) - doc = xsdataObject.to_xml() - - # Reduce number of lines - lines = doc.toprettyxml(indent=' ') - lines = lines.replace('\n ', '') - lines = lines.replace('\n ', '') - lines = lines.replace('\n ', '') - lines = lines.replace('\n ', '') - lines = lines.replace('\n ', '') - lines = lines.replace('\n ', '') - lines = lines.replace('\n ', '') - lines = lines.replace('\n ', '') - lines = lines.replace('\n ', '') - lines = lines.replace('\n ', '') - lines = lines.replace('\n ', '') - lines = lines.replace('\n ', '') - lines = lines.replace('\n ', '') - lines = lines.replace('\n ', '') - lines = lines.replace('\n ', '') - lines = lines.replace('\n ', '') - - # Write document in pretty XML to specified file - f = open(xmlFile, 'w') - f.write(lines) - f.close() diff --git a/scripts/openmc-xsdir-to-xml b/scripts/openmc-xsdir-to-xml deleted file mode 100755 index 7e1606fc4..000000000 --- a/scripts/openmc-xsdir-to-xml +++ /dev/null @@ -1,288 +0,0 @@ -#!/usr/bin/env python - -import os -import sys -from xml.dom.minidom import getDOMImplementation - -elements = [None, "H", "He", "Li", "Be", "B", "C", "N", "O", "F", "Ne", "Na", - "Mg", "Al", "Si", "P", "S", "Cl", "Ar", "K", "Ca", "Sc", "Ti", "V", - "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga", "Ge", "As", "Se", - "Br", "Kr", "Rb", "Sr", "Y", "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", - "Pd", "Ag", "Cd", "In", "Sn", "Sb", "Te", "I", "Xe", "Cs", "Ba", - "La", "Ce", "Pr", "Nd", "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", - "Er", "Tm", "Yb", "Lu", "Hf", "Ta", "W", "Re", "Os", "Ir", "Pt", - "Au", "Hg", "Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac", - "Th", "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm", - "Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds", "Rg", - "Cn"] - - -class Xsdir(object): - - def __init__(self, filename): - self.f = open(filename, 'r') - self.filename = os.path.abspath(filename) - self.directory = os.path.dirname(filename) - self.awr = {} - self.tables = [] - - self.filetype = set() - self.recordlength = set() - self.entries = set() - - # Read first section (DATAPATH) - line = self.f.readline() - words = line.split() - if words: - if words[0].lower().startswith('datapath'): - if '=' in words[0]: - index = line.index('=') - self.datapath = line[index+1:].strip() - else: - if len(line.strip()) > 8: - self.datapath = line[8:].strip() - else: - self.f.seek(0) - - # Read second section - line = self.f.readline() - words = line.split() - assert len(words) == 3 - assert words[0].lower() == 'atomic' - assert words[1].lower() == 'weight' - assert words[2].lower() == 'ratios' - - while True: - line = self.f.readline() - words = line.split() - - # Check for end of second section - if len(words) % 2 != 0 or words[0] == 'directory': - break - - for zaid, awr in zip(words[::2], words[1::2]): - self.awr[zaid] = awr - - # Read third section - while words[0] != 'directory': - words = self.f.readline().split() - - while True: - words = self.f.readline().split() - if not words: - break - - # Handle continuation lines - while words[-1] == '+': - extraWords = self.f.readline().split() - words = words[:-1] + extraWords - assert len(words) >= 7 - - # Create XsdirTable object and add to line - table = XsdirTable(self.directory) - self.tables.append(table) - - # All tables have at least 7 attributes - table.name = words[0] - table.awr = float(words[1]) - table.filename = words[2] - table.access = words[3] - table.filetype = int(words[4]) - table.location = int(words[5]) - table.length = int(words[6]) - - self.filetype.add(table.filetype) - - if len(words) > 7: - table.recordlength = int(words[7]) - self.recordlength.add(table.recordlength) - if len(words) > 8: - table.entries = int(words[8]) - self.entries.add(table.entries) - if len(words) > 9: - table.temperature = float(words[9]) - if len(words) > 10: - table.ptable = (words[10] == 'ptable') - - if len(self.filetype) == 1: - if 1 in self.filetype: - self.filetype = 'ascii' - elif 2 in self.filetype: - self.filetype = 'binary' - else: - self.filetype = None - - if len(self.recordlength) == 1: - self.recordlength = list(self.recordlength)[0] - else: - self.recordlength = None - if len(self.entries) == 1: - self.entries = list(self.entries)[0] - else: - self.recordlength = None - - def to_xml(self): - # Create XML document - impl = getDOMImplementation() - doc = impl.createDocument(None, "cross_sections", None) - - # Get root element - root = doc.documentElement - - # Add a directory node - if self.directory: - directoryNode = doc.createElement("directory") - text = doc.createTextNode(self.directory) - directoryNode.appendChild(text) - root.appendChild(directoryNode) - - for table in self.tables: - table.path = os.path.basename(table.path) - - # Add filetype, record_length and entries nodes - if self.filetype: - node = doc.createElement("filetype") - text = doc.createTextNode(self.filetype) - node.appendChild(text) - root.appendChild(node) - if self.recordlength: - node = doc.createElement("record_length") - text = doc.createTextNode(str(self.recordlength)) - node.appendChild(text) - root.appendChild(node) - if self.entries: - node = doc.createElement("entries") - text = doc.createTextNode(str(self.entries)) - node.appendChild(text) - root.appendChild(node) - - # Add a node for each table - for table in self.tables: - if table.name[-1] in ['e', 'p', 'u', 'h', 'g', 'm', 'd']: - continue - node = table.to_xml_node(doc) - root.appendChild(node) - - return doc - - -class XsdirTable(object): - - def __init__(self, directory=None): - self.directory = None - self.name = None - self.awr = None - self.filename = None - self.access = None - self.filetype = None - self.location = None - self.length = None - self.recordlength = None - self.entries = None - self.temperature = None - self.ptable = False - - @property - def path(self): - if self.directory: - return os.path.join(self.directory, self.filename) - else: - return self.filename - - @path.setter - def path(self, value): - self.diretory = '' - self.filename = value - - @property - def metastable(self): - # Only valid for neutron cross-sections - if not self.name.endswith('c'): - return - - # Handle special case of Am-242 and Am-242m - if self.zaid == '95242': - return 1 - elif self.zaid == '95642': - return 0 - - # All other cases - A = int(self.zaid) % 1000 - if A > 300: - return 1 - else: - return 0 - - @property - def alias(self): - zaid = self.zaid - if zaid: - Z = int(zaid[:-3]) - A = zaid[-3:] - - if A == '000': - s = 'Nat' - elif zaid == '95242': - s = '242m' - elif zaid == '95642': - s = '242' - elif int(A) > 300: - s = str(int(A) - 400) + "m" - else: - s = str(int(A)) - - return "{0}-{1}.{2}".format(elements[Z], s, self.xs) - else: - return None - - @property - def zaid(self): - if self.name.endswith('c'): - return self.name[:self.name.find('.')] - else: - return 0 - - @property - def xs(self): - return self.name[self.name.find('.')+1:] - - def to_xml_node(self, doc): - node = doc.createElement("ace_table") - node.setAttribute("name", self.name) - for attribute in ["alias", "zaid", "type", "metastable", "awr", - "temperature", "path", "location"]: - if hasattr(self, attribute): - string = str(getattr(self, attribute)) - - # Skip metastable and binary if 0 - if attribute == "metastable" and self.metastable == 0: - continue - - # Skip any attribute that is none - if getattr(self, attribute) is None: - continue - - # Create attribute node - node.setAttribute(attribute, string) - - return node - - -if __name__ == '__main__': - # Read command line arguments - if len(sys.argv) < 3: - sys.exit("Usage: convert_xsdir.py xsdirFile xmlFile") - xsdirFile = sys.argv[1] - xmlFile = sys.argv[2] - - # Read xsdata and create XML document object - xsdirObject = Xsdir(xsdirFile) - doc = xsdirObject.to_xml() - - # Reduce number of lines - lines = doc.toprettyxml(indent=' ') - - # Write document in pretty XML to specified file - f = open(xmlFile, 'w') - f.write(lines) - f.close() diff --git a/src/reaction_header.F90 b/src/reaction_header.F90 index 8af75b3cf..74198dd78 100644 --- a/src/reaction_header.F90 +++ b/src/reaction_header.F90 @@ -38,7 +38,7 @@ contains integer(HSIZE_T) :: dims(1) call read_attribute(this % Q_value, group_id, 'Q_value') - call read_attribute(this % MT, group_id, 'MT') + call read_attribute(this % MT, group_id, 'mt') call read_attribute(this % threshold, group_id, 'threshold_idx') call read_attribute(cm, group_id, 'center_of_mass') this % scatter_in_cm = (cm == 1) diff --git a/src/relaxng/cross_sections.rnc b/src/relaxng/cross_sections.rnc index 75c0ea29c..7fbc610a2 100644 --- a/src/relaxng/cross_sections.rnc +++ b/src/relaxng/cross_sections.rnc @@ -1,25 +1,12 @@ element cross_sections { - element ace_table { - (element name { xsd:string { maxLength = "15" } } | - attribute name { xsd:string { maxLength = "15" } }) & - (element alias { xsd:string { maxLength = "15" } } | - attribute alias { xsd:string { maxLength = "15" } })? & - (element zaid { xsd:int } | attribute zaid { xsd:int }) & - (element metastable { xsd:int } | attribute metastable { xsd:int })? & - (element awr { xsd:double } | attribute awr { xsd:double }) & - (element temperature { xsd:double } | attribute temperature { xsd:double }) & - (element path { xsd:string { maxLength = "255" } } | - attribute path { xsd:string { maxLength = "255" } }) & - (element location { xsd:int } | attribute location { xsd:int })? & - (element filetype { ( "ascii" | "binary" ) } | - attribute filetype { ( "ascii" | "binary" ) })? + element library { + (element materials { xsd:string } | + attribute materials { xsd:string }) & + (element type { xsd:string } | + attribute type { xsd:string }) & + (element path { xsd:string } | + attribute path { xsd:string }) }* & - element directory { xsd:string { maxLength = "255" } }? & - - element filetype { ( "ascii" | "binary" ) } & - - element record_length { xsd:int }? & - - element entries { xsd:int }? + element directory { xsd:string { maxLength = "255" } }? } \ No newline at end of file diff --git a/src/relaxng/cross_sections.rng b/src/relaxng/cross_sections.rng index 5e531e013..435f7fa84 100644 --- a/src/relaxng/cross_sections.rng +++ b/src/relaxng/cross_sections.rng @@ -2,106 +2,32 @@ - + - - - 15 - + + - - - 15 - - - - - - - - 15 - - - - - 15 - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + - - + + - - + + - - 255 - + - - 255 - + - - - - - - - - - - - - - - - ascii - binary - - - - - ascii - binary - - - - @@ -112,21 +38,5 @@ - - - ascii - binary - - - - - - - - - - - -