diff --git a/docs/source/pythonapi/data.rst b/docs/source/pythonapi/data.rst index 7f75ddaaa..1c62db487 100644 --- a/docs/source/pythonapi/data.rst +++ b/docs/source/pythonapi/data.rst @@ -19,6 +19,9 @@ Core Classes openmc.data.CoherentElastic openmc.data.FissionEnergyRelease openmc.data.DataLibrary + openmc.data.IncidentPhoton + openmc.data.PhotonReaction + openmc.data.AtomicRelaxation openmc.data.Decay openmc.data.FissionProductYields openmc.data.WindowedMultipole diff --git a/openmc/data/__init__.py b/openmc/data/__init__.py index 7158e9fe3..6dd6a9218 100644 --- a/openmc/data/__init__.py +++ b/openmc/data/__init__.py @@ -9,6 +9,7 @@ WMP_VERSION = 'v0.2' from .data import * from .neutron import * +from .photon import * from .decay import * from .reaction import * from .ace import * diff --git a/openmc/data/compton_profiles.h5 b/openmc/data/compton_profiles.h5 new file mode 100644 index 000000000..e4fa2e56d Binary files /dev/null and b/openmc/data/compton_profiles.h5 differ diff --git a/openmc/data/endf.py b/openmc/data/endf.py index 0fa75ded8..37dd435e0 100644 --- a/openmc/data/endf.py +++ b/openmc/data/endf.py @@ -23,10 +23,30 @@ from .function import Tabulated1D, INTERPOLATION_SCHEME from openmc.stats.univariate import Uniform, Tabular, Legendre -LIBRARIES = {0: 'ENDF/B', 1: 'ENDF/A', 2: 'JEFF', 3: 'EFF', - 4: 'ENDF/B High Energy', 5: 'CENDL', 6: 'JENDL', - 31: 'INDL/V', 32: 'INDL/A', 33: 'FENDL', 34: 'IRDF', - 35: 'BROND', 36: 'INGDB-90', 37: 'FENDL/A', 41: 'BROND'} +_LIBRARY = {0: 'ENDF/B', 1: 'ENDF/A', 2: 'JEFF', 3: 'EFF', + 4: 'ENDF/B High Energy', 5: 'CENDL', 6: 'JENDL', + 31: 'INDL/V', 32: 'INDL/A', 33: 'FENDL', 34: 'IRDF', + 35: 'BROND', 36: 'INGDB-90', 37: 'FENDL/A', 41: 'BROND'} + +_SUBLIBRARY = { + 0: 'Photo-nuclear data', + 1: 'Photo-induced fission product yields', + 3: 'Photo-atomic data', + 4: 'Radioactive decay data', + 5: 'Spontaneous fission product yields', + 6: 'Atomic relaxation data', + 10: 'Incident-neutron data', + 11: 'Neutron-induced fission product yields', + 12: 'Thermal neutron scattering data', + 19: 'Neutron standards', + 113: 'Electro-atomic data', + 10010: 'Incident-proton data', + 10011: 'Proton-induced fission product yields', + 10020: 'Incident-deuteron data', + 10030: 'Incident-triton data', + 20030: 'Incident-helion (3He) data', + 20040: 'Incident-alpha data' +} SUM_RULES = {1: [2, 3], 3: [4, 5, 11, 16, 17, 22, 23, 24, 25, 27, 28, 29, 30, 32, 33, 34, 35, @@ -348,6 +368,14 @@ class Evaluation(object): self._read_header() + def __repr__(self): + if 'zsymam' in self.target: + name = self.target['zsymam'].replace(' ', '') + else: + name = 'Unknown' + return '<{} for {} {}>'.format(self.info['sublibrary'], name, + self.info['library']) + def _read_header(self): file_obj = io.StringIO(self.section[1, 451]) @@ -360,8 +388,7 @@ class Evaluation(object): self._LRP = items[2] self.target['fissionable'] = (items[3] == 1) try: - global LIBRARIES - library = LIBRARIES[items[4]] + library = _LIBRARY[items[4]] except KeyError: library = 'Unknown' self.info['modification'] = items[5] @@ -384,7 +411,7 @@ class Evaluation(object): self.projectile['mass'] = items[0] self.info['energy_max'] = items[1] library_release = items[2] - self.info['sublibrary'] = items[4] + self.info['sublibrary'] = _SUBLIBRARY[items[4]] library_version = items[5] self.info['library'] = (library, library_version, library_release) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 1d3a0176d..95f919369 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -104,20 +104,23 @@ def _get_metadata(zaid, metastable_scheme='nndc'): class IncidentNeutron(EqualityMixin): """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`. + This class stores data derived from an ENDF-6 format neutron interaction + sublibrary. Instances of this class are not normally instantiated by the + user but rather created using the factory methods + :meth:`IncidentNeutron.from_hdf5`, :meth:`IncidentNeutron.from_ace`, and + :math:`IncidentNeutron.from_endf`. Parameters ---------- name : str Name of the nuclide using the GND naming convention atomic_number : int - Number of protons in the nucleus + Number of protons in the target nucleus mass_number : int - Number of nucleons in the nucleus + Number of nucleons in the target nucleus metastable : int - Metastable state of the nucleus. A value of zero indicates ground state. + Metastable state of the target nucleus. A value of zero indicates ground + state. atomic_weight_ratio : float Atomic mass ratio of the target nuclide. kTs : Iterable of float @@ -127,22 +130,19 @@ class IncidentNeutron(EqualityMixin): Attributes ---------- atomic_number : int - Number of protons in the nucleus + Number of protons in the target nucleus atomic_symbol : str Atomic symbol of the nuclide, e.g., 'Zr' atomic_weight_ratio : float Atomic weight ratio of the target nuclide. - energy : dict of numpy.ndarray - The energy values (eV) at which reaction cross-sections are tabulated. - They keys of the dict are the temperature string ('294K') for each - set of energies fission_energy : None or openmc.data.FissionEnergyRelease The energy released by fission, tabulated by component (e.g. prompt neutrons or beta particles) and dependent on incident neutron energy mass_number : int - Number of nucleons in the nucleus + Number of nucleons in the target nucleus metastable : int - Metastable state of the nucleus. A value of zero indicates ground state. + Metastable state of the target nucleus. A value of zero indicates ground + state. name : str Name of the nuclide using the GND naming convention reactions : collections.OrderedDict diff --git a/openmc/data/photon.py b/openmc/data/photon.py new file mode 100644 index 000000000..e87a8792b --- /dev/null +++ b/openmc/data/photon.py @@ -0,0 +1,642 @@ +from collections import OrderedDict, Mapping, Callable +from copy import deepcopy +from io import StringIO +from numbers import Integral, Real +import os + +import h5py +import numpy as np +import pandas as pd + +from openmc.mixin import EqualityMixin +import openmc.checkvalue as cv +from . import HDF5_VERSION +from .data import ATOMIC_SYMBOL +from .endf import Evaluation, get_head_record, get_tab1_record, get_list_record +from .function import Tabulated1D + + +_SUBSHELLS = ['K', 'L1', 'L2', 'L3', 'M1', 'M2', 'M3', 'M4', 'M5', + 'N1', 'N2', 'N3', 'N4', 'N5', 'N6', 'N7', 'O1', 'O2', + 'O3', 'O4', 'O5', 'O6', 'O7', 'O8', 'O9', 'P1', 'P2', + 'P3', 'P4', 'P5', 'P6', 'P7', 'P8', 'P9', 'P10', 'P11', + 'Q1', 'Q2', 'Q3'] + +_REACTION_NAME = { + 501: 'Total photon interaction', + 502: 'Photon coherent scattering', + 504: 'Photon incoherent scattering', + 515: 'Pair production, electron field', + 516: 'Total pair production', + 517: 'Pair production, nuclear field', + 522: 'Photoelectric absorption', + 526: 'Electro-atomic scattering', + 527: 'Electro-atomic bremsstrahlung', + 528: 'Electro-atomic excitation', + 534: 'K (1s1/2) subshell photoelectric', + 535: 'L1 (2s1/2) subshell photoelectric', + 536: 'L2 (2p1/2) subshell photoelectric', + 537: 'L3 (2p3/2) subshell photoelectric', + 538: 'M1 (3s1/2) subshell photoelectric', + 539: 'M2 (3p1/2) subshell photoelectric', + 540: 'M3 (3p3/2) subshell photoelectric', + 541: 'M4 (3d3/2) subshell photoelectric', + 542: 'M5 (3d5/2) subshell photoelectric', + 543: 'N1 (4s1/2) subshell photoelectric', + 544: 'N2 (4p1/2) subshell photoelectric', + 545: 'N3 (4p3/2) subshell photoelectric', + 546: 'N4 (4d3/2) subshell photoelectric', + 547: 'N5 (4d5/2) subshell photoelectric', + 548: 'N6 (4f5/2) subshell photoelectric', + 549: 'N7 (4f7/2) subshell photoelectric', + 550: 'O1 (5s1/2) subshell photoelectric', + 551: 'O2 (5p1/2) subshell photoelectric', + 552: 'O3 (5p3/2) subshell photoelectric', + 553: 'O4 (5d3/2) subshell photoelectric', + 554: 'O5 (5d5/2) subshell photoelectric', + 555: 'O6 (5f5/2) subshell photoelectric', + 556: 'O7 (5f7/2) subshell photoelectric', + 557: 'O8 (5g7/2) subshell photoelectric', + 558: 'O9 (5g9/2) subshell photoelectric', + 559: 'P1 (6s1/2) subshell photoelectric', + 560: 'P2 (6p1/2) subshell photoelectric', + 561: 'P3 (6p3/2) subshell photoelectric', + 562: 'P4 (6d3/2) subshell photoelectric', + 563: 'P5 (6d5/2) subshell photoelectric', + 564: 'P6 (6f5/2) subshell photoelectric', + 565: 'P7 (6f7/2) subshell photoelectric', + 566: 'P8 (6g7/2) subshell photoelectric', + 567: 'P9 (6g9/2) subshell photoelectric', + 568: 'P10 (6h9/2) subshell photoelectric', + 569: 'P11 (6h11/2) subshell photoelectric', + 570: 'Q1 (7s1/2) subshell photoelectric', + 571: 'Q2 (7p1/2) subshell photoelectric', + 572: 'Q3 (7p3/2) subshell photoelectric' +} + +# Compton profiles are read from a pre-generated HDF5 file when they are first +# needed. The dictionary stores an array of electron momentum values (at which +# the profiles are tabulated) with the key 'pz' and the profile for each element +# is a 2D array with shape (n_shells, n_momentum_values) stored on the key Z +_COMPTON_PROFILES = {} + + +class AtomicRelaxation(EqualityMixin): + """Atomic relaxation data. + + This class stores the binding energy, number of electrons, and electron + transitions possible from ioniziation for each subshell with an atom. All of + the data originates from an ENDF-6 atomic relaxation sub-library + (NSUB=6). Instances of this class are not normally instantiated directly but + rather created using the factory method :math:`AtomicRelaxation.from_endf`. + + Parameters + ---------- + binding_energy : dict + Dictionary indicating the binding energy in eV (values) for given + subshells (keys). The subshells should be given as strings, e.g., 'K', + 'L1', 'L2', etc. + num_electrons : dict + Dictionary indicating the number of electrons in a subshell when neutral + (values) for given subshells (keys). The subshells should be given as + strings, e.g., 'K', 'L1', 'L2', etc. + transitions : pandas.DataFrame + Dictionary indicating allowed transitions and their probabilities + (values) for given subshells (keys). The subshells should be given as + strings, e.g., 'K', 'L1', 'L2', etc. The transitions are represented as + a DataFrame with columns indicating the secondary and tertiary subshell, + the energy of the transition in eV, and the fractional probability of + the transition. + + Attributes + ---------- + binding_energy : dict + Dictionary indicating the binding energy in eV (values) for given + subshells (keys). The subshells should be given as strings, e.g., 'K', + 'L1', 'L2', etc. + num_electrons : dict + Dictionary indicating the number of electrons in a subshell when neutral + (values) for given subshells (keys). The subshells should be given as + strings, e.g., 'K', 'L1', 'L2', etc. + transitions : pandas.DataFrame + Dictionary indicating allowed transitions and their probabilities + (values) for given subshells (keys). The subshells should be given as + strings, e.g., 'K', 'L1', 'L2', etc. The transitions are represented as + a DataFrame with columns indicating the secondary and tertiary subshell, + the energy of the transition in eV, and the fractional probability of + the transition. + + See Also + -------- + IncidentPhoton + + """ + def __init__(self, binding_energy, num_electrons, transitions): + self.binding_energy = binding_energy + self.num_electrons = num_electrons + self.transitions = transitions + self.compton_profile = OrderedDict() + + @property + def binding_energy(self): + return self._binding_energy + + @property + def num_electrons(self): + return self._num_electrons + + @property + def subshells(self): + return list(sorted(self.binding_energy.keys())) + + @property + def transitions(self): + return self._transitions + + @binding_energy.setter + def binding_energy(self, binding_energy): + cv.check_type('binding energies', binding_energy, Mapping) + for subshell, energy in binding_energy.items(): + cv.check_value('subshell', subshell, _SUBSHELLS) + cv.check_type('binding energy', energy, Real) + cv.check_greater_than('binding energy', energy, 0.0, True) + self._binding_energy = binding_energy + + @num_electrons.setter + def num_electrons(self, num_electrons): + cv.check_type('number of electrons', num_electrons, Mapping) + for subshell, num in num_electrons.items(): + cv.check_value('subshell', subshell, _SUBSHELLS) + cv.check_type('number of electrons', num, Real) + cv.check_greater_than('number of electrons', num, 0.0, True) + self._num_electrons = num_electrons + + @transitions.setter + def transitions(self, transitions): + cv.check_type('transitions', transitions, Mapping) + for subshell, df in transitions.items(): + cv.check_value('subshell', subshell, _SUBSHELLS) + cv.check_type('transitions', df, pd.DataFrame) + self._transitions = transitions + + @classmethod + def from_endf(cls, ev_or_filename): + """Generate atomic relaxation data from an ENDF evaluation + + Parameters + ---------- + ev_or_filename : str or openmc.data.endf.Evaluation + ENDF atomic relaxation evaluation to read from. If given as a + string, it is assumed to be the filename for the ENDF file. + + Returns + ------- + openmc.data.AtomicRelaxation + Atomic relaxation data + + """ + if isinstance(ev_or_filename, Evaluation): + ev = ev_or_filename + else: + ev = Evaluation(ev_or_filename) + + # Atomic relaxation data is always MF=28, MT=533 + if (28, 533) not in ev.section: + raise IOError('{} does not appear to be an atomic relaxation ' + 'sublibrary.'.format(ev)) + + # Determine number of subshells + file_obj = StringIO(ev.section[28, 533]) + params = get_head_record(file_obj) + n_subshells = params[4] + + # Helper function to map designator to subshell string or None + def subshell(i): + if i == 0: + return None + else: + return _SUBSHELLS[i - 1] + + # Create data dictionaries + binding_energy = {} + num_electrons = {} + transitions = {} + columns = ['secondary', 'tertiary', 'energy (eV)', 'probability'] + + # Read data for each subshell + for i in range(n_subshells): + params, list_items = get_list_record(file_obj) + subi = subshell(int(params[0])) + n_transitions = int(params[5]) + binding_energy[subi] = list_items[0] + num_electrons[subi] = list_items[1] + + if n_transitions > 0: + # Read transition data + records = [] + for j in range(n_transitions): + subj = subshell(int(list_items[6*(j+1)])) + subk = subshell(int(list_items[6*(j+1) + 1])) + etr = list_items[6*(j+1) + 2] + ftr = list_items[6*(j+1) + 3] + records.append((subj, subk, etr, ftr)) + + # Create dataframe for transitions + transitions[subi] = pd.DataFrame.from_records( + records, columns=columns) + + # Return instance of class + data = cls(binding_energy, num_electrons, transitions) + + return data + + def to_hdf5(self, group): + raise NotImplementedError + + +class IncidentPhoton(EqualityMixin): + """Photon interaction data. + + This class stores photo-atomic, photo-nuclear, atomic relaxation, and + Compton profile data assembled from different sources. To create an + instance, the factory method :meth:`IncidentPhoton.from_endf` can be + used. To add atomic relaxation or Compton profile data, set the + :attr:`IncidentPhoton.atomic_relaxation` and + :attr:`IncidentPhoton.compton_profiles` attributes directly. + + Parameters + ---------- + atomic_number : int + Number of protons in the target nucleus + + Attributes + ---------- + atomic_number : int + Number of protons in the target nucleus + atomic_relaxation : openmc.data.AtomicRelaxation or None + Atomic relaxation data + compton_profiles : dict + Dictionary of Compton profile data with keys 'num_electrons' (number of + electrons in each subshell), 'binding_energy' (ionization potential of + each subshell), and 'J' (Hartree-Fock Compton profile as a function of + the projection of the electron momentum on the scattering vector, + :math:`p_z` for each subshell). Note that subshell occupancies may not + match the atomic relaxation data. + reactions : collections.OrderedDict + Contains the cross sections for each photon reaction. The keys are MT + values and the values are instances of :class:`PhotonReaction`. + summed_reactions : collections.OrderedDict + Contains summed cross sections. The keys are MT values and the values + are instances of :class:`PhotonReaction`. + + """ + + def __init__(self, atomic_number): + self.atomic_number = atomic_number + self._atomic_relaxation = None + self.reactions = OrderedDict() + self.summed_reactions = OrderedDict() + self.compton_profiles = {} + + def __contains__(self, mt): + return mt in self.reactions or mt in self.summed_reactions + + def __getitem__(self, mt): + if mt in self.reactions: + return self.reactions[mt] + elif mt in self.summed_reactions: + return self.summed_reactions[mt] + else: + raise KeyError('No reaction with MT={}.'.format(mt)) + + def __repr__(self): + return "".format(self.name) + + def __iter__(self): + return iter(self.reactions.values()) + + @property + def atomic_number(self): + return self._atomic_number + + @property + def atomic_relaxation(self): + return self._atomic_relaxation + + @property + def name(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, True) + self._atomic_number = atomic_number + + @atomic_relaxation.setter + def atomic_relaxation(self, atomic_relaxation): + cv.check_type('atomic relaxation data', atomic_relaxation, + AtomicRelaxation) + self._atomic_relaxation = atomic_relaxation + + @classmethod + def from_endf(cls, photoatomic, relaxation=None): + """Generate incident photon data from an ENDF evaluation + + Parameters + ---------- + photoatomic : str or openmc.data.endf.Evaluation + ENDF photoatomic data evaluation to read from. If given as a string, + it is assumed to be the filename for the ENDF file. + relaxation : str or openmc.data.endf.Evaluation, optional + ENDF atomic relaxation data evaluation to read from. If given as a + string, it is assumed to be the filename for the ENDF file. + + Returns + ------- + openmc.data.IncidentPhoton + Photon interaction data + + """ + if isinstance(photoatomic, Evaluation): + ev = photoatomic + else: + ev = Evaluation(photoatomic) + + Z = ev.target['atomic_number'] + data = cls(Z) + + # Read each reaction + for mf, mt, nc, mod in ev.reaction_list: + if mf == 23: + data.reactions[mt] = PhotonReaction.from_endf(ev, mt) + + # Add atomic relaxation data if it hasn't been added already + if relaxation is not None: + data.atomic_relaxation = AtomicRelaxation.from_endf(relaxation) + + # If Compton profile data hasn't been loaded, do so + if not _COMPTON_PROFILES: + filename = os.path.join(os.path.dirname(__file__), 'compton_profiles.h5') + with h5py.File(filename, 'r') as f: + _COMPTON_PROFILES['pz'] = f['pz'].value + for i in range(1, 101): + group = f['{:03}'.format(i)] + num_electrons = group['num_electrons'].value + binding_energy = group['binding_energy'].value + J = group['J'].value + _COMPTON_PROFILES[i] = {'num_electrons': num_electrons, + 'binding_energy': binding_energy, + 'J': J} + + # Add Compton profile data + pz = _COMPTON_PROFILES['pz'] + profile = _COMPTON_PROFILES[Z] + data.compton_profiles['num_electrons'] = profile['num_electrons'] + data.compton_profiles['binding_energy'] = profile['binding_energy'] + data.compton_profiles['J'] = [Tabulated1D(pz, J_k) for J_k in profile['J']] + + return data + + def export_to_hdf5(self, path, mode='a'): + """Export incident photon data 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. + + """ + # Open file and write version + f = h5py.File(path, mode, libver='latest') + if 'version' not in f.attrs: + f.attrs['version'] = np.array(HDF5_VERSION) + + group = f.create_group(self.name) + group.attrs['Z'] = Z = self.atomic_number + + # Determine union energy grid + union_grid = np.array([]) + for rx in self: + union_grid = np.union1d(union_grid, rx.xs.x) + group.create_dataset('energy', data=union_grid) + + # Write coherent scattering cross section + rx = self.reactions[502] + coh_group = group.create_group('coherent') + coh_group.create_dataset('xs', data=rx.xs(union_grid)) + if rx.scattering_factor is not None: + # Create integrated form factor + ff = deepcopy(rx.scattering_factor) + ff.x *= ff.x + ff.y *= ff.y/Z**2 + int_ff = Tabulated1D(ff.x, ff.integral()) + int_ff.to_hdf5(coh_group, 'integrated_scattering_factor') + if rx.anomalous_real is not None: + rx.anomalous_real.to_hdf5(coh_group, 'anomalous_real') + if rx.anomalous_imag is not None: + rx.anomalous_imag.to_hdf5(coh_group, 'anomalous_imag') + + # Write incoherent scattering cross section + rx = self[504] + incoh_group = group.create_group('incoherent') + incoh_group.create_dataset('xs', data=rx.xs(union_grid)) + if rx.scattering_factor is not None: + rx.scattering_factor.to_hdf5(incoh_group, 'scattering_factor') + + # Write pair production cross section + if 515 in self: + pair_group = group.create_group('pair_production') + pair_group.create_dataset('xs', data=self[515].xs(union_grid)) + + # Write triplet 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)) + + # Write photoelectric cross section + photoelec_group = group.create_group('photoelectric') + photoelec_group.create_dataset('xs', data=self[522].xs(union_grid)) + + # Write photoionization cross sections + shell_group = group.create_group('subshells') + designators = [] + for mt, rx in self.reactions.items(): + if mt >= 534 and mt <= 572: + # Get name of subshell + shell = _SUBSHELLS[mt - 534] + designators.append(shell) + sub_group = shell_group.create_group(shell) + + if self.atomic_relaxation is not None: + relax = self.atomic_relaxation + # Write subshell binding energy and number of electrons + sub_group.attrs['binding_energy'] = relax.binding_energy[shell] + sub_group.attrs['num_electrons'] = relax.num_electrons[shell] + + # Write transition data with replacements + if shell in relax.transitions: + shell_values = _SUBSHELLS.copy() + shell_values.insert(0, None) + df = relax.transitions[shell].replace( + shell_values, range(len(shell_values))) + sub_group.create_dataset('transitions', data=df.as_matrix()) + + # Determine threshold + threshold = rx.xs.x[0] + idx = np.searchsorted(union_grid, threshold, side='right') - 1 + + # Interpolate cross section onto union grid and write + photoionization = rx.xs(union_grid[idx:]) + sub_group.create_dataset('xs', data=photoionization) + assert len(union_grid) == len(photoionization) + idx + sub_group['xs'].attrs['threshold_idx'] = idx + + shell_group.attrs['designators'] = np.array(designators, dtype='S') + + # Write Compton profiles + if self.compton_profiles: + compton_group = group.create_group('compton_profiles') + + profile = self.compton_profiles + compton_group.create_dataset('num_electrons', + data=profile['num_electrons']) + compton_group.create_dataset('binding_energy', + data=profile['binding_energy']) + + # Get electron momentum values + compton_group.create_dataset('pz', data=profile['J'][0].x) + + # Create/write 2D array of profiles + J = np.array([Jk.y for Jk in profile['J']]) + compton_group.create_dataset('J', data=J) + + +class PhotonReaction(EqualityMixin): + """Photon-induced reaction + + Parameters + ---------- + mt : int + The ENDF MT number for this reaction. + + Attributes + ---------- + anomalous_real : openmc.data.Tabulated1D + Real part of the anomalous scattering factor + anomlaous_imag : openmc.data.Tabulated1D + Imaginary part of the anomalous scatttering factor + mt : int + The ENDF MT number for this reaction. + scattering_factor : openmc.data.Tabulated1D + Coherent or incoherent form factor. + xs : Callable + Cross section as a function of incident photon energy + + """ + + def __init__(self, mt): + self.mt = mt + self._xs = None + self._scattering_factor = None + self._anomalous_real = None + self._anomalous_imag = None + + def __repr__(self): + if self.mt in _REACTION_NAME: + return "".format( + self.mt, _REACTION_NAME[self.mt]) + else: + return "".format(self.mt) + + @property + def anomalous_real(self): + return self._anomalous_real + + @property + def anomalous_imag(self): + return self._anomalous_imag + + @property + def scattering_factor(self): + return self._scattering_factor + + @property + def xs(self): + return self._xs + + @anomalous_real.setter + def anomalous_real(self, anomalous_real): + cv.check_type('real part of anomalous scattering factor', + anomalous_real, Callable) + self._anomalous_real = anomalous_real + + @anomalous_imag.setter + def anomalous_imag(self, anomalous_imag): + cv.check_type('imaginary part of anomalous scattering factor', + anomalous_imag, Callable) + self._anomalous_imag = anomalous_imag + + @scattering_factor.setter + def scattering_factor(self, scattering_factor): + cv.check_type('scattering factor', scattering_factor, Callable) + self._scattering_factor = scattering_factor + + @xs.setter + def xs(self, xs): + cv.check_type('reaction cross section', xs, Callable) + self._xs = xs + + @classmethod + def from_endf(self, ev, mt): + """Generate photon reaction from an ENDF evaluation + + Parameters + ---------- + ev : openmc.data.endf.Evaluation + ENDF photo-atomic interaction data evaluation + mt : int + The MT value of the reaction to get data for + + Returns + ------- + openmc.data.PhotonReaction + Photon reaction data + + """ + rx = PhotonReaction(mt) + + # Read photon cross section + if (23, mt) in ev.section: + file_obj = StringIO(ev.section[23, mt]) + get_head_record(file_obj) + params, rx.xs = get_tab1_record(file_obj) + + # Set subshell binding energy and/or fluorescence yield + if mt >= 534 and mt <= 599: + rx.subshell_binding_energy = params[0] + if mt >= 534 and mt <= 572: + rx.fluorescence_yield = params[1] + + # Read form factors / scattering functions + if (27, mt) in ev.section: + file_obj = StringIO(ev.section[27, mt]) + get_head_record(file_obj) + params, rx.scattering_factor = get_tab1_record(file_obj) + + # Check for anomalous scattering factor + if mt == 502: + if (27, 506) in ev.section: + file_obj = StringIO(ev.section[27, 506]) + get_head_record(file_obj) + params, rx.anomalous_real = get_tab1_record(file_obj) + + if (27, 505) in ev.section: + file_obj = StringIO(ev.section[27, 505]) + get_head_record(file_obj) + params, rx.anomalous_imag = get_tab1_record(file_obj) + + return rx diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index ac3a049ab..e3e1094e9 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -1095,7 +1095,7 @@ class Reaction(EqualityMixin): ev : openmc.data.endf.Evaluation ENDF evaluation mt : int - The MT value of the reaction to get angular distributions for + The MT value of the reaction to get data for Returns ------- diff --git a/scripts/openmc-make-compton b/scripts/openmc-make-compton new file mode 100755 index 000000000..1f920f366 --- /dev/null +++ b/scripts/openmc-make-compton @@ -0,0 +1,104 @@ +#!/usr/bin/env python + +from __future__ import print_function, division +import os +import sys +import tarfile + +from six.moves import input +from six.moves.urllib.request import urlopen +import numpy as np +import h5py + + +base_url = 'http://geant4.cern.ch/support/source/' +filename = 'G4EMLOW.6.48.tar.gz' +block_size = 16384 + +# ============================================================================== +# DOWNLOAD FILES FROM GEANT4 SITE + +# Establish connection to URL +req = urlopen(base_url + filename) + +# Get file size from header +if sys.version_info[0] < 3: + file_size = int(req.info().getheaders('Content-Length')[0]) +else: + file_size = req.length +downloaded = 0 + +# Check if file already downloaded +download = True +if os.path.exists(filename): + if os.path.getsize(filename) == file_size: + print('Already downloaded ' + filename) + download = False + else: + overwrite = input('Overwrite {}? ([y]/n) '.format(filename)) + if overwrite.lower().startswith('n'): + download = False + +if download: + # Copy file to disk + print('Downloading {}... '.format(filename), end='') + with open(filename, 'wb') as fh: + while True: + chunk = req.read(block_size) + if not chunk: break + fh.write(chunk) + downloaded += len(chunk) + status = '{0:10} [{1:3.2f}%]'.format( + downloaded, downloaded * 100. / file_size) + print(status + chr(8)*len(status), end='') + print('') + +# ============================================================================== +# EXTRACT FILES FROM TGZ + +if not os.path.isdir('G4EMLOW6.48'): + with tarfile.open(filename, 'r') as tgz: + print('Extracting {0}...'.format(filename)) + tgz.extractall() + +# ============================================================================== +# GENERATE COMPTON PROFILE HDF5 FILE + +print('Generating compton_profiles.h5...') + +shell_file = os.path.join('G4EMLOW6.48', 'doppler', 'shell-doppler.dat') + +with open(shell_file, 'r') as shell: + with h5py.File('compton_profiles.h5', 'w') as f: + # Read/write electron momentum values + pz = np.loadtxt(os.path.join('G4EMLOW6.48', 'doppler', 'p-biggs.dat')) + f.create_dataset('pz', data=pz) + + for Z in range(1, 101): + # Create group for this element + group = f.create_group('{:03}'.format(Z)) + + # Read data into one long array + path = os.path.join('G4EMLOW6.48', 'doppler', 'profile-{}.dat'.format(Z)) + J = np.fromstring(open(path, 'r').read(), sep=' ') + + # Determine number of electron shells and reshape + n_shells = J.size // 31 + J.shape = (n_shells, 31) + + # Write Compton profile for this Z + group.create_dataset('J', data=J) + + # Determine binding energies and number of electrons for each shell + num_electrons = [] + binding_energy = [] + while True: + words = shell.readline().split() + if words[0] == '-1': + break + num_electrons.append(float(words[0])) + binding_energy.append(float(words[1])) + + # Write binding energies and number of electrons + group.create_dataset('num_electrons', data=num_electrons) + group.create_dataset('binding_energy', data=binding_energy) diff --git a/setup.py b/setup.py index 67bc73403..d0c6bd5b2 100755 --- a/setup.py +++ b/setup.py @@ -68,6 +68,14 @@ if have_setuptools: 'validate': ['lxml'] }, + # Data files + 'package_data': { + 'openmc.data': [ + 'mass.mas12', + 'fission_Q_data_endfb71.h5', + 'compton_profiles.h5' + ] + }, }) # If Cython is present, add resonance reconstruction capability