From 04308b7677aa3ec775d90dd54e38f18c32762c9f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 10 Jan 2020 12:59:30 -0600 Subject: [PATCH] Rewrite of openmc-ace-to-hdf5, moving some functionality to ace.py --- docs/source/pythonapi/data.rst | 3 + openmc/data/ace.py | 127 +++++++++++++- scripts/openmc-ace-to-hdf5 | 231 +++++++++----------------- tests/unit_tests/test_data_neutron.py | 13 ++ 4 files changed, 218 insertions(+), 156 deletions(-) diff --git a/docs/source/pythonapi/data.rst b/docs/source/pythonapi/data.rst index be29d7c0e..48797ef5d 100644 --- a/docs/source/pythonapi/data.rst +++ b/docs/source/pythonapi/data.rst @@ -148,6 +148,7 @@ Classes ace.Library ace.Table + ace.TableType Functions +++++++++ @@ -158,6 +159,8 @@ Functions :template: myfunction.rst ace.ascii_to_binary + ace.get_libraries_from_xsdir + ace.get_libraries_from_xsdata ENDF Format ----------- diff --git a/openmc/data/ace.py b/openmc/data/ace.py index ee194560b..8ebf2ddb7 100644 --- a/openmc/data/ace.py +++ b/openmc/data/ace.py @@ -15,7 +15,8 @@ generates ACE-format cross sections. """ -from pathlib import PurePath +import enum +from pathlib import PurePath, Path import struct import sys @@ -358,7 +359,7 @@ class Library(EqualityMixin): atomic_weight_ratio = float(words[0]) temperature = float(words[1]) commentlines = int(words[3]) - for i in range(commentlines): + for _ in range(commentlines): lines.pop(0) lines.append(ace_file.readline()) else: @@ -383,7 +384,7 @@ class Library(EqualityMixin): # verify that we are supposed to read this table in if (table_names is not None) and (name not in table_names): - for i in range(n_lines - 1): + for _ in range(n_lines - 1): ace_file.readline() lines = [ace_file.readline() for i in range(13)] continue @@ -423,6 +424,42 @@ class Library(EqualityMixin): lines = [ace_file.readline() for i in range(13)] +class TableType(enum.Enum): + """Type of ACE data table.""" + NEUTRON_CONTINUOUS = 'c' + NEUTRON_DISCRETE = 'd' + THERMAL_SCATTERING = 't' + DOSIMETRY = 'y' + PHOTOATOMIC = 'p' + PHOTONUCLEAR = 'u' + PROTON = 'h' + DEUTERON = 'o' + TRITON = 'r' + HELIUM3 = 's' + ALPHA = 'a' + + @classmethod + def from_suffix(cls, suffix): + """Determine ACE table type from a suffix. + + Parameters + ---------- + suffix : str + Single letter ACE table designator, e.g., 'c' + + Returns + ------- + TableType + ACE table type + + """ + for member in cls: + if suffix.endswith(member.value): + return member + raise ValueError("Suffix '{}' has no corresponding ACE table type." + .format(suffix)) + + class Table(EqualityMixin): """ACE cross section table @@ -445,6 +482,11 @@ class Table(EqualityMixin): xss : numpy.ndarray Raw data for the ACE table + Attributes + ---------- + data_type : TableType + Type of the ACE data + """ def __init__(self, name, atomic_weight_ratio, temperature, pairs, nxs, jxs, xss): @@ -456,5 +498,84 @@ class Table(EqualityMixin): self.jxs = jxs self.xss = xss + @property + def zaid(self): + return self.name.split('.')[0] + + @property + def data_type(self): + xs = self.name.split('.')[1] + return TableType.from_suffix(xs[-1]) + def __repr__(self): return "".format(self.name) + + +def get_libraries_from_xsdir(path): + """Determine paths to ACE files from an MCNP xsdir file. + + Parameters + ---------- + path : str or path-like + Path to xsdir file + + Returns + ------- + list + List of paths to ACE libraries + """ + xsdir = Path(path) + + # Find 'directory' section + with open(path, 'r') as fh: + lines = fh.readlines() + for index, line in enumerate(lines): + if line.strip().lower() == 'directory': + break + else: + raise OSError("Could not find 'directory' section in MCNP xsdir file") + + # Handle continuation lines indicated by '+' at end of line + lines = lines[index + 1:] + continue_lines = [i for i, line in enumerate(lines) + if line.strip().endswith('+')] + for i in reversed(continue_lines): + lines[i] += lines[i].strip()[:-1] + lines.pop(i + 1) + + # Create list of ACE libraries + libraries = [] + for line in lines: + words = line.split() + if len(words) < 3: + continue + + lib = (xsdir.parent / words[2]).resolve() + if lib not in libraries: + libraries.append(lib) + + return libraries + + +def get_libraries_from_xsdata(path): + """Determine paths to ACE files from a Serpent xsdata file. + + Parameters + ---------- + path : str or path-like + Path to xsdata file + + Returns + ------- + list + List of paths to ACE libraries + """ + xsdata = Path(path) + with open(xsdata, 'r') as xsdata: + libraries = [] + for line in xsdata: + words = line.split() + if len(words) >= 9: + lib = (xsdata.parent / words[8]).resolve() + if lib not in libraries: + libraries.append(lib) + return libraries diff --git a/scripts/openmc-ace-to-hdf5 b/scripts/openmc-ace-to-hdf5 index 9953efa8c..3aeb95198 100755 --- a/scripts/openmc-ace-to-hdf5 +++ b/scripts/openmc-ace-to-hdf5 @@ -1,49 +1,49 @@ #!/usr/bin/env python3 -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 +"""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. +2. Use the --xsdir option to specify a MCNP xsdir file. +3. 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). +The script does not use any extra information from 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). """ +import argparse +from functools import partial +import os +from pathlib import Path +import warnings + +import openmc.data +from openmc.data.ace import TableType + + class CustomFormatter(argparse.ArgumentDefaultsHelpFormatter, argparse.RawDescriptionHelpFormatter): pass + parser = argparse.ArgumentParser( - description=description, + description=__doc__, formatter_class=CustomFormatter ) parser.add_argument('libraries', nargs='*', help='ACE libraries to convert to HDF5') -parser.add_argument('-d', '--destination', default='.', +parser.add_argument('-d', '--destination', type=Path, default=Path.cwd(), 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 ' @@ -54,152 +54,77 @@ parser.add_argument('--libver', choices=['earliest', 'latest'], "performance") args = parser.parse_args() -if not os.path.isdir(args.destination): - os.mkdir(args.destination) +if not args.destination.is_dir(): + args.destination.mkdir(parents=True, exist_ok=True) -# 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'): - path = os.path.join(directory, ace_table.attrib['path']) - if path not in ace_libraries: - ace_libraries.append(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") - - # Handle continuation lines indicated by '+' at end of line - lines = lines[index + 1:] - continue_lines = [i for i, line in enumerate(lines) - if line.strip().endswith('+')] - for i in reversed(continue_lines): - lines[i] += lines[i].strip()[:-1] + lines.pop(i + 1) - - # Create list of ACE libraries - for line in lines: - 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) - +if args.xsdir is not None: + ace_libraries.extend(openmc.data.ace.get_libraries_from_xsdir(args.xsdir)) 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) - + ace_libraries.extend(openmc.data.ace.get_libraries_from_xsdata(args.xsdata)) else: - ace_libraries = args.libraries + ace_libraries = [Path(lib) for lib in args.libraries] -nuclides = {} +converted = {} library = openmc.data.DataLibrary() -for filename in ace_libraries: +for path in ace_libraries: # Check that ACE library exists - if not os.path.exists(filename): - warnings.warn("ACE library '{}' does not exist.".format(filename)) + if not os.path.exists(path): + warnings.warn("ACE library '{}' does not exist.".format(path)) continue - lib = openmc.data.ace.Library(filename) + lib = openmc.data.ace.Library(path) for table in lib.tables: - name, xs = table.name.split('.') - if xs.endswith('c'): - # Continuous-energy neutron data - if name not in nuclides: - try: - neutron = openmc.data.IncidentNeutron.from_ace( - table, args.metastable) - except Exception as e: - print('Failed to convert {}: {}'.format(table.name, e)) - continue - - print('Converting {} (ACE) to {} (HDF5)'.format(table.name, - neutron.name)) - - # Determine filename - outfile = os.path.join(args.destination, - neutron.name.replace('.', '_') + '.h5') - neutron.export_to_hdf5(outfile, 'w', libver=args.libver) - - # Register with library - library.register_file(outfile) - - # Add nuclide to list - nuclides[name] = outfile - else: - # Then we only need to append the data - try: - neutron = \ - openmc.data.IncidentNeutron.from_hdf5(nuclides[name]) - print('Converting {} (ACE) to {} (HDF5)' - .format(table.name, neutron.name)) - neutron.add_temperature_from_ace(table, args.metastable) - neutron.export_to_hdf5(nuclides[name] + '_1', 'w', - libver=args.libver) - os.rename(nuclides[name] + '_1', nuclides[name]) - except Exception as e: - print('Failed to convert {}: {}'.format(table.name, e)) - continue - - elif xs.endswith('t'): + # Check type of the ACE table and determine appropriate class / + # conversion function + if table.data_type == TableType.NEUTRON_CONTINUOUS: + name = table.zaid + cls = openmc.data.IncidentNeutron + converter = partial(cls.from_ace, metastable_scheme=args.metastable) + elif table.data_type == TableType.THERMAL_SCATTERING: # Adjust name to be the new thermal scattering name - name = openmc.data.get_thermal_name(name) - # Thermal scattering data - if name not in nuclides: - try: - thermal = openmc.data.ThermalScattering.from_ace(table) - except Exception as e: - print('Failed to convert {}: {}'.format(table.name, e)) - continue - print('Converting {} (ACE) to {} (HDF5)'.format(table.name, - thermal.name)) + name = openmc.data.get_thermal_name(table.zaid) + cls = openmc.data.ThermalScattering + converter = cls.from_ace + else: + print("Can't convert ACE table {}".format(table.name)) + continue - # Determine filename - outfile = os.path.join(args.destination, - thermal.name.replace('.', '_') + '.h5') - thermal.export_to_hdf5(outfile, 'w', libver=args.libver) + if name not in converted: + try: + data = converter(table) + except Exception as e: + print('Failed to convert {}: {}'.format(table.name, e)) + continue - # Register with library - library.register_file(outfile) + print('Converting {} (ACE) to {} (HDF5)'.format(table.name, data.name)) - # Add data to list - nuclides[name] = outfile + # Determine output filename + outfile = args.destination / (data.name.replace('.', '_') + '.h5') + data.export_to_hdf5(outfile, 'w', libver=args.libver) - else: - # Then we only need to append the data - try: - thermal = openmc.data.ThermalScattering.from_hdf5( - nuclides[name]) - print('Converting {} (ACE) to {} (HDF5)' - .format(table.name, thermal.name)) - thermal.add_temperature_from_ace(table) - thermal.export_to_hdf5(nuclides[name] + '_1', 'w', - libver=args.libver) - os.rename(nuclides[name] + '_1', nuclides[name]) - except Exception as e: - print('Failed to convert {}: {}'.format(table.name, e)) - continue + # Register with library + library.register_file(outfile) + + # Add nuclide to list + converted[name] = outfile + else: + # Read existing HDF5 file + data = cls.from_hdf5(converted[name]) + + # Add data for new temperature + try: + print('Converting {} (ACE) to {} (HDF5)' + .format(table.name, data.name)) + data.add_temperature_from_ace(table, args.metastable) + except Exception as e: + print('Failed to convert {}: {}'.format(table.name, e)) + continue + + # Re-export + data.export_to_hdf5(converted[name] + '_1', 'w', libver=args.libver) + os.rename(converted[name] + '_1', converted[name]) # Write cross_sections.xml -libpath = os.path.join(args.destination, 'cross_sections.xml') -library.export_to_xml(libpath) +library.export_to_xml(args.destination / 'cross_sections.xml') diff --git a/tests/unit_tests/test_data_neutron.py b/tests/unit_tests/test_data_neutron.py index ac6164e5b..48ea117d8 100644 --- a/tests/unit_tests/test_data_neutron.py +++ b/tests/unit_tests/test_data_neutron.py @@ -517,3 +517,16 @@ def test_ace_convert(run_in_tmpdir): assert tab_a.temperature == pytest.approx(tab_b.temperature) assert np.all(tab_a.nxs == tab_b.nxs) assert np.all(tab_a.jxs == tab_b.jxs) + assert tab_a.zaid == tab_b.zaid + assert tab_a.data_type == tab_b.data_type + + +def test_ace_table_types(): + TT = openmc.data.ace.TableType + assert TT.from_suffix('c') == TT.NEUTRON_CONTINUOUS + assert TT.from_suffix('nc') == TT.NEUTRON_CONTINUOUS + assert TT.from_suffix('80c') == TT.NEUTRON_CONTINUOUS + assert TT.from_suffix('t') == TT.THERMAL_SCATTERING + assert TT.from_suffix('20t') == TT.THERMAL_SCATTERING + with pytest.raises(ValueError): + TT.from_suffix('z')