diff --git a/openmc/data/__init__.py b/openmc/data/__init__.py index 6dd6a9218..aaf98c574 100644 --- a/openmc/data/__init__.py +++ b/openmc/data/__init__.py @@ -12,10 +12,10 @@ from .neutron import * from .photon import * from .decay import * from .reaction import * -from .ace import * +from . import ace from .angle_distribution import * from .function import * -from .endf import * +from . import endf from .energy_distribution import * from .product import * from .angle_energy import * diff --git a/openmc/data/ace.py b/openmc/data/ace.py index 385408bd4..cbbc4e06d 100644 --- a/openmc/data/ace.py +++ b/openmc/data/ace.py @@ -16,13 +16,84 @@ generates ACE-format cross sections. """ from os import SEEK_CUR +from pathlib import PurePath import struct import sys import numpy as np from openmc.mixin import EqualityMixin -from openmc.data.endf import ENDF_FLOAT_RE +import openmc.checkvalue as cv +from .data import ATOMIC_SYMBOL +from .endf import ENDF_FLOAT_RE + + +def get_metadata(zaid, metastable_scheme='nndc'): + """Return basic identifying data for a nuclide with a given ZAID. + + Parameters + ---------- + zaid : int + ZAID (1000*Z + A) obtained from a library + 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 + ------- + name : str + Name of the table + element : str + The atomic symbol of the isotope in the table; e.g., Zr. + Z : 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. + + """ + + cv.check_type('zaid', zaid, int) + cv.check_value('metastable_scheme', metastable_scheme, ['nndc', 'mcnp']) + + 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 + element = ATOMIC_SYMBOL[Z] + name = '{}{}'.format(element, mass_number) + if metastable > 0: + name += '_m{}'.format(metastable) + + return (name, element, Z, mass_number, metastable) + def ascii_to_binary(ascii_file, binary_file): """Convert an ACE file in ASCII format (type 1) to binary format (type 2). @@ -160,7 +231,7 @@ class Library(EqualityMixin): # Determine whether file is ASCII or binary try: - fh = open(filename, 'rb') + fh = open(str(filename), 'rb') # Grab 10 lines of the library sb = b''.join([fh.readline() for i in range(10)]) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 99847be44..3a76cfe2a 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -14,7 +14,7 @@ import numpy as np import h5py from . import HDF5_VERSION, HDF5_VERSION_MAJOR -from .ace import Library, Table, get_table +from .ace import Library, Table, get_table, get_metadata from .data import ATOMIC_SYMBOL, K_BOLTZMANN, EV_PER_MEV from .endf import Evaluation, SUM_RULES, get_head_record, get_tab1_record from .fission_energy import FissionEnergyRelease @@ -33,73 +33,6 @@ from openmc.mixin import EqualityMixin _RESONANCE_ENERGY_GRID = np.logspace(-3, 3, 61) -def _get_metadata(zaid, metastable_scheme='nndc'): - """Return basic identifying data for a nuclide with a given ZAID. - - Parameters - ---------- - zaid : int - ZAID (1000*Z + A) obtained from a library - 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 - ------- - name : str - Name of the table - element : str - The atomic symbol of the isotope in the table; e.g., Zr. - Z : 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. - - """ - - cv.check_type('zaid', zaid, int) - cv.check_value('metastable_scheme', metastable_scheme, ['nndc', 'mcnp']) - - 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 - element = ATOMIC_SYMBOL[Z] - name = '{}{}'.format(element, mass_number) - if metastable > 0: - name += '_m{}'.format(metastable) - - return (name, element, Z, mass_number, metastable) - - class IncidentNeutron(EqualityMixin): """Continuous-energy neutron interaction data. @@ -676,7 +609,7 @@ class IncidentNeutron(EqualityMixin): # If mass number hasn't been specified, make an educated guess zaid, xs = ace.name.split('.') name, element, Z, mass_number, metastable = \ - _get_metadata(int(zaid), metastable_scheme) + get_metadata(int(zaid), metastable_scheme) # Assign temperature to the running list kTs = [ace.temperature*EV_PER_MEV] diff --git a/openmc/data/photon.py b/openmc/data/photon.py index be04ec911..4fe529900 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -12,6 +12,7 @@ from scipy.interpolate import CubicSpline from openmc.mixin import EqualityMixin import openmc.checkvalue as cv from . import HDF5_VERSION +from .ace import Table, get_metadata, get_table from .data import ATOMIC_SYMBOL, EV_PER_MEV from .endf import Evaluation, get_head_record, get_tab1_record, get_list_record from .function import Tabulated1D @@ -96,6 +97,7 @@ _STOPPING_POWERS = {} # for each element are in a 2D array with shape (n, k) stored on the key 'Z'. _BREMSSTRAHLUNG = {} + class AtomicRelaxation(EqualityMixin): """Atomic relaxation data. @@ -272,7 +274,7 @@ class AtomicRelaxation(EqualityMixin): class IncidentPhoton(EqualityMixin): """Photon interaction data. - This class stores photo-atomic, photo-nuclear, atomic relaxation, + This class stores photo-atomic, photo-nuclear, atomic relaxation, Compton profile, stopping power, and bremsstrahlung data assembled from different sources. To create an instance, the factory method :meth:`IncidentPhoton.from_endf` can be used. To add atomic relaxation or @@ -369,6 +371,68 @@ class IncidentPhoton(EqualityMixin): AtomicRelaxation) self._atomic_relaxation = atomic_relaxation + @classmethod + def from_ace(cls, ace_or_filename): + """Generate incident photon data from an ACE table + + Parameters + ---------- + ace_or_filename : str or openmc.data.ace.Table + ACE table to read from. If given as a string, it is assumed to be + the filename for the ACE file. + + Returns + ------- + openmc.data.IncidentPhoton + Photon interaction data + + """ + # First obtain the data for the first provided ACE table/file + if isinstance(ace_or_filename, Table): + ace = ace_or_filename + else: + ace = get_table(ace_or_filename) + + # Get atomic number based on name of ACE table + zaid = ace.name.split('.')[0] + Z = get_metadata(int(zaid))[2] + + # Read each reaction + data = cls(Z) + for mt in (502, 504, 515, 522): + data.reactions[mt] = PhotonReaction.from_ace(ace, mt) + + # Compton profiles + n_shell = ace.nxs[5] + if n_shell != 0: + # Get number of electrons in each shell + idx = ace.jxs[6] + data.compton_profiles['num_electrons'] = ace.xss[idx : idx+n_shell] + + # Get binding energy for each shell + idx = ace.jxs[7] + data.compton_profiles['binding_energy'] = ace.xss[idx : idx+n_shell] + + # Create Compton profile for each electron shell + profiles = [] + for k in range(n_shell): + # Get number of momentum values and interpolation scheme + loca = int(ace.xss[ace.jxs[9] + k]) + jj = int(ace.xss[ace.jxs[10] + loca - 1]) + m = int(ace.xss[ace.jxs[10] + loca]) + + # Read momentum and PDF + idx = ace.jxs[10] + loca + 1 + pz = ace.xss[idx : idx+m] + pdf = ace.xss[idx+m : idx+2*m] + + # Create proflie function + J_k = Tabulated1D(pz, pdf, [m], [jj]) + profiles.append(J_k) + data.compton_profiles['J'] = profiles + + return data + @classmethod def from_endf(cls, photoatomic, relaxation=None): """Generate incident photon data from an ENDF evaluation @@ -548,15 +612,15 @@ class IncidentPhoton(EqualityMixin): if rx.scattering_factor is not None: rx.scattering_factor.to_hdf5(incoh_group, 'scattering_factor') - # Write pair production cross section + # Write electron-field pair production cross section if 515 in self: - pair_group = group.create_group('pair_production') + pair_group = group.create_group('pair_production_electron') pair_group.create_dataset('xs', data=self[515].xs(union_grid)) - # Write triplet production cross section + # Write nuclear-field pair production cross section if 517 in self: - triplet_group = group.create_group('triplet_production') - triplet_group.create_dataset('xs', data=self[517].xs(union_grid)) + pair_group = group.create_group('pair_production_nuclear') + pair_group.create_dataset('xs', data=self[517].xs(union_grid)) # Write photoelectric cross section photoelec_group = group.create_group('photoelectric') @@ -713,7 +777,89 @@ class PhotonReaction(EqualityMixin): self._xs = xs @classmethod - def from_endf(self, ev, mt): + def from_ace(cls, ace, mt): + """Generate photon reaction from an ACE table + + Parameters + ---------- + ace : openmc.data.ace.Table + ACE table to read from + mt : int + The MT value of the reaction to get data for + + Returns + ------- + openmc.data.PhotonReaction + Photon reaction data + + """ + # Create instance + rx = cls(mt) + + # Get energy grid (stored as logarithms) + n = ace.nxs[3] + idx = ace.jxs[1] + energy = np.exp(ace.xss[idx : idx+n])*EV_PER_MEV + + # Get index for appropriate reaction + if mt == 502: + # Coherent scattering + idx = ace.jxs[1] + 2*n + elif mt == 504: + # Incoherent scattering + idx = ace.jxs[1] + n + elif mt == 515: + # Pair production + idx = ace.jxs[1] + 4*n + elif mt == 522: + # Photoelectric + idx = ace.jxs[1] + 3*n + else: + raise ValueError('ACE photoatomic cross sections do not have ' + 'data for MT={}.'.format(mt)) + + # Store cross section + xs = np.exp(ace.xss[idx : idx+n]) + rx.xs = Tabulated1D(energy, xs, [n], [5]) + + # Get form factors for incoherent/coherent scattering + if mt == 502: + idx = ace.jxs[3] + if ace.nxs[6] > 0: + n = (ace.jxs[4] - ace.jxs[3]) // 2 + x = ace.xss[idx : idx+n] + idx += n + else: + x = np.array([ + 0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.08, 0.1, 0.12, + 0.15, 0.18, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, + 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, + 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, + 3.6, 3.8, 4.0, 4.2, 4.4, 4.6, 4.8, 5.0, 5.2, 5.4, 5.6, + 5.8, 6.0]) + n = x.size + ff = ace.xss[idx+n : idx+2*n] + rx.scattering_factor = Tabulated1D(x, ff) + + elif mt == 504: + idx = ace.jxs[2] + if ace.nxs[6] > 0: + n = (ace.jxs[3] - ace.jxs[2]) // 2 + x = ace.xss[idx : idx+n] + idx += n + else: + x = np.array([ + 0.0, 0.005, 0.01, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6, + 0.7, 0.8, 0.9, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 8.0 + ]) + n = x.size + ff = ace.xss[idx : idx+n] + rx.scattering_factor = Tabulated1D(x, ff) + + return rx + + @classmethod + def from_endf(cls, ev, mt): """Generate photon reaction from an ENDF evaluation Parameters @@ -729,7 +875,7 @@ class PhotonReaction(EqualityMixin): Photon reaction data """ - rx = PhotonReaction(mt) + rx = cls(mt) # Read photon cross section if (23, mt) in ev.section: diff --git a/scripts/openmc-update-inputs b/scripts/openmc-update-inputs index ef7cdf47b..0d25b5b97 100755 --- a/scripts/openmc-update-inputs +++ b/scripts/openmc-update-inputs @@ -248,8 +248,7 @@ def update_materials(root): # If a nuclide name is in the ZAID notation (e.g., a number), # convert it to the proper nuclide name. if nucname.strip().isnumeric(): - nucname = \ - openmc.data.neutron._get_metadata(int(nucname))[0] + nucname = openmc.data.ace.get_metadata(int(nucname))[0] nucname = nucname.replace('Nat', '0') if nucname.endswith('m'): nucname = nucname[:-1] + '_m1' diff --git a/src/photon_header.F90 b/src/photon_header.F90 index 27fe62f18..e50a871d2 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -179,14 +179,18 @@ contains call close_group(rgroup) ! Read pair production - rgroup = open_group(group_id, 'pair_production') - call read_dataset(this % pair_production_nuclear, rgroup, 'xs') + rgroup = open_group(group_id, 'pair_production_electron') + call read_dataset(this % pair_production_electron, rgroup, 'xs') call close_group(rgroup) ! Read pair production - rgroup = open_group(group_id, 'triplet_production') - call read_dataset(this % pair_production_electron, rgroup, 'xs') - call close_group(rgroup) + if (object_exists(group_id, 'pair_production_nuclear')) then + rgroup = open_group(group_id, 'pair_production_nuclear') + call read_dataset(this % pair_production_nuclear, rgroup, 'xs') + call close_group(rgroup) + else + this % pair_production_nuclear(:) = ZERO + end if ! Read photoelectric rgroup = open_group(group_id, 'photoelectric')