From a38f53d26f5585cd2d312567a3f69baf369dbe9a Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 29 Jul 2016 10:50:08 -0500 Subject: [PATCH 01/33] Add Madland fission-Q support to openmc.data --- openmc/data/__init__.py | 1 + openmc/data/endf_utils.py | 43 +++++ openmc/data/fission_energy.py | 285 ++++++++++++++++++++++++++++++++++ openmc/data/neutron.py | 22 +++ 4 files changed, 351 insertions(+) create mode 100644 openmc/data/endf_utils.py create mode 100644 openmc/data/fission_energy.py diff --git a/openmc/data/__init__.py b/openmc/data/__init__.py index ae8ea8191..e60878d60 100644 --- a/openmc/data/__init__.py +++ b/openmc/data/__init__.py @@ -14,3 +14,4 @@ from .nbody import * from .thermal import * from .urr import * from .library import * +from .fission_energy import * diff --git a/openmc/data/endf_utils.py b/openmc/data/endf_utils.py new file mode 100644 index 000000000..6db3c611c --- /dev/null +++ b/openmc/data/endf_utils.py @@ -0,0 +1,43 @@ +"""This module contains a few utility functions for reading ENDF_ data. It is by +no means enough to read an entire ENDF file. For a more complete ENDF reader, +see Pyne_. + +.. _ENDF: http://www.nndc.bnl.gov/endf +.. _Pyne: http://www.pyne.io + +""" + +import re + +def read_float(float_string): + """Parse ENDF 6E11.0 formatted string into a float.""" + assert len(float_string) == 11 + pattern = '([\s\\-]\d+\\.\d+)([\\+\\-]\d+)' + mantissa, exponent = re.match(pattern, float_string).groups() + return float(mantissa + 'e' + exponent) + + +def read_CONT_line(line): + """Parse 80-column line from ENDF CONT record into floats and ints.""" + return (read_float(line[0:11]), read_float(line[11:22]), int(line[22:33]), + int(line[33:44]), int(line[44:55]), int(line[55:66]), + int(line[66:70]), int(line[70:72]), int(line[72:75]), + int(line[75:80])) + +def identify_nuclide(fname): + """Read the header of an ENDF file and extract identifying information.""" + with open(fname, 'r') as fh: + # Skip the tape id (TPID). + line = fh.readline() + + # Read the first HEAD and CONT info. + line = fh.readline() + ZA, AW, LRP, LFI, NLIB, NMOD, MAT, MF, MT, NS = read_CONT_line(line) + line = fh.readline() + ELIS, STA, LIS, LISO, junk, NFOR, MAT, MF, MT, NS = read_CONT_line(line) + + # Return dictionary of the most important identifying information. + return {'Z': int(ZA) // 1000, + 'A': int(ZA) % 1000, + 'LIS': LIS, + 'LISO': LISO} diff --git a/openmc/data/fission_energy.py b/openmc/data/fission_energy.py new file mode 100644 index 000000000..5716a3e85 --- /dev/null +++ b/openmc/data/fission_energy.py @@ -0,0 +1,285 @@ +from collections import Callable +import sys +#from warnings import warn + +import numpy as np +from numpy.polynomial.polynomial import Polynomial + +from .function import Tabulated1D, Sum +from .endf_utils import read_float, read_CONT_line, identify_nuclide +import openmc.checkvalue as cv + +if sys.version_info[0] >= 3: + basestring = str + + +class FissionEnergyRelease(object): + def __init__(self): + self._fragments = None + self._prompt_neutrons = None + self._delayed_neutrons = None + self._prompt_photons = None + self._delayed_photons = None + self._betas = None + self._neutrinos = None + self._form = None + + @property + def fragments(self): + return self._fragments + + @property + def prompt_neutrons(self): + return self._prompt_neutrons + + @property + def delayed_neutrons(self): + return self._delayed_neutrons + + @property + def prompt_photons(self): + return self._prompt_photons + + @property + def delayed_photons(self): + return self._delayed_photons + + @property + def betas(self): + return self._betas + + @property + def neutrinos(self): + return self._neutrinos + + @property + def recoverable(self): + return Sum([self.fragments, self.prompt_neutrons, self.delayed_neutrons, + self.prompt_photons, self.delayed_photons, self.betas]) + + @property + def total(self): + return Sum([self.fragments, self.prompt_neutrons, self.delayed_neutrons, + self.prompt_photons, self.delayed_photons, self.betas, + self.neutrinos]) + + @property + def form(self): + return self._form + + @fragments.setter + def fragments(self, energy_release): + cv.check_type('fragments', energy_release, Callable) + self._fragments = energy_release + + @prompt_neutrons.setter + def prompt_neutrons(self, energy_release): + cv.check_type('prompt_neutrons', energy_release, Callable) + self._prompt_neutrons = energy_release + + @delayed_neutrons.setter + def delayed_neutrons(self, energy_release): + cv.check_type('delayed_neutrons', energy_release, Callable) + self._delayed_neutrons = energy_release + + @prompt_photons.setter + def prompt_photons(self, energy_release): + cv.check_type('prompt_photons', energy_release, Callable) + self._prompt_photons = energy_release + + @delayed_photons.setter + def delayed_photons(self, energy_release): + cv.check_type('delayed_photons', energy_release, Callable) + self._delayed_photons = energy_release + + @betas.setter + def betas(self, energy_release): + cv.check_type('betas', energy_release, Callable) + self._betas = energy_release + + @neutrinos.setter + def neutrinos(self, energy_release): + cv.check_type('neutrinos', energy_release, Callable) + self._neutrinos = energy_release + + @form.setter + def form(self, form): + cv.check_value('format', form, ('Madland', 'Sher-Beck')) + self._form = form + + @classmethod + def from_endf(cls, filename, incident_neutron): + """Generate fission energy release data from an ENDF file. + + Parameters + ---------- + filename : str + Name of the ENDF file containing fission energy release data + + incident_neutron : openmc.data.IncidentNeutron + Corresponding incident neutron dataset + + Returns + ------- + openmc.data.FissionEnergyRelease + Fission energy release data + + """ + + # Check to make sure this ENDF file matches the expected isomer. + ident = identify_nuclide(filename) + if ident['Z'] != incident_neutron.atomic_number: + pass + if ident['A'] != incident_neutron.mass_number: + pass + if ident['LISO'] != incident_neutron.metastable: + pass + + # Extract the MF=1, MT=458 section. + lines = [] + with open(filename, 'r') as fh: + line = fh.readline() + while line != '': + if line[70:75] == ' 1458': + lines.append(line) + line = fh.readline() + + # Read the number of coefficients in this LIST record. + NPL = read_CONT_line(lines[1])[4] + + # Parse the ENDF LIST into an array. + data = [] + for i in range(NPL): + row, column = divmod(i, 6) + data.append(read_float(lines[2 + row][11*column:11*(column+1)])) + + # Declare the coefficient names and the order they are given in. The + # LIST contains a value followed immediately by an uncertainty for each + # of these components, times the polynomial order + 1. If we only find + # one value for each of these components, then we need to use the + # Sher-Beck formula for energy dependence. Otherwise, it is a + # polynomial. + labels = ('EFR', 'ENP', 'END', 'EGP', 'EGD', 'EB', 'ENU', 'ER', 'ET') + + # Associate each set of values and uncertainties with its label. + value = dict() + uncertainty = dict() + for i in range(len(labels)): + value[labels[i]] = data[2*i::18] + uncertainty[labels[i]] = data[2*i + 1::18] + + # In ENDF/B-7.1, data for 2nd-order coefficients were mistakenly not + # converted from MeV to eV. Check for this error and fix it if present. + n_coeffs = len(value['EFR']) + if n_coeffs == 3: # Only check 2nd-order data. + # Check each energy component for the error. If a 1 MeV neutron + # causes a change of more than 100 MeV, we know something is wrong. + error_present = False + for coeffs in value.values(): + second_order = coeffs[2] + if abs(second_order) * 1e12 > 1e8: + error_present = True + break + + # If we found the error, reduce all 2nd-order coeffs by 10**6. + if error_present: + for coeffs in value.values(): coeffs[2] *= 1e-6 + for coeffs in uncertainty.values(): coeffs[2] *= 1e-6 + + # Perform the sanity check again... just in case. + for coeffs in value.values(): + second_order = coeffs[2] + if abs(second_order) * 1e12 > 1e8: + raise ValueError("Encountered a ludicrously large second-" + "order polynomial coefficient.") + + # Convert eV to MeV. + for coeffs in value.values(): + for i in range(len(coeffs)): + coeffs[i] *= 10**(-6 + 6*i) + for coeffs in uncertainty.values(): + for i in range(len(coeffs)): + coeffs[i] *= 10**(-6 + 6*i) + + out = cls() + if n_coeffs > 1: + out.form = 'Madland' + out.fragments = Polynomial(value['EFR']) + out.prompt_neutrons = Polynomial(value['ENP']) + out.delayed_neutrons = Polynomial(value['END']) + out.prompt_photons = Polynomial(value['EGP']) + out.delayed_photons = Polynomial(value['EGD']) + out.betas = Polynomial(value['EB']) + out.neutrinos = Polynomial(value['ENU']) + else: + out.form = 'Sher-Beck' + raise NotImplemented + + return out + + @classmethod + def from_hdf5(cls, group): + """Generate fission energy release data from an HDF5 group. + + Parameters + ---------- + group : h5py.Group + HDF5 group to read from + + Returns + ------- + openmc.data.FissionEnergyRelease + Fission energy release data + + """ + + obj = cls() + if group.attrs['format'] == 'Madland': + obj.fragments = Polynomial(group['fragments'].value) + obj.prompt_neutrons = Polynomial(group['prompt_neutrons'].value) + obj.delayed_neutrons = Polynomial(group['delayed_neutrons'].value) + obj.prompt_photons = Polynomial(group['prompt_photons'].value) + obj.delayed_photons = Polynomial(group['delayed_photons'].value) + obj.betas = Polynomial(group['betas'].value) + obj.neutrinos = Polynomial(group['neutrinos'].value) + elif group.attrs['format'] == 'Sher-Beck': + raise NotImplemented + else: + raise ValueError('Unrecognized energy release format') + + return obj + + def to_hdf5(self, group): + """Write energy release data to an HDF5 group + + Parameters + ---------- + group : h5py.Group + HDF5 group to write to + + """ + + if self.form == 'Madland': + group.attrs['format'] = np.string_('Madland') + group.create_dataset('fragments', data=self.fragments.coef) + group.create_dataset('prompt_neutrons', + data=self.prompt_neutrons.coef) + group.create_dataset('delayed_neutrons', + data=self.delayed_neutrons.coef) + group.create_dataset('prompt_photons', + data=self.prompt_photons.coef) + group.create_dataset('delayed_photons', + data=self.delayed_photons.coef) + group.create_dataset('betas', data=self.betas.coef) + group.create_dataset('neutrinos', data=self.neutrinos.coef) + elif self.form == 'Sher-Beck': + group.attrs['format'] = np.string_('Sher-Beck') + self.fragments.to_hdf5(group, 'fragments') + self.prompt_neutrons.to_hdf5(group, 'prompt_neutrons') + self.delayed_neutrons.to_hdf5(group, 'delayed_neutrons') + self.prompt_photons.to_hdf5(group, 'prompt_photons') + self.delayed_photons.to_hdf5(group, 'delayed_photons') + self.betas.to_hdf5(group, 'betas') + self.neutrinos.to_hdf5(group, 'neutrinos') + else: + raise ValueError('Unrecognized energy release format') diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 63b2bab01..84d9bcb2c 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -9,6 +9,7 @@ import h5py from .data import ATOMIC_SYMBOL, SUM_RULES from .ace import Table, get_table +from .fission_energy import FissionEnergyRelease from .function import Tabulated1D, Sum from .product import Product from .reaction import Reaction, _get_photon_products @@ -81,6 +82,7 @@ class IncidentNeutron(object): self.temperature = temperature self._energy = None + self._fission_energy = None self.reactions = OrderedDict() self.summed_reactions = OrderedDict() self.urr = None @@ -126,6 +128,10 @@ class IncidentNeutron(object): def energy(self): return self._energy + @property + def fission_energy(self): + return self._fission_energy + @property def temperature(self): return self._temperature @@ -186,6 +192,12 @@ class IncidentNeutron(object): cv.check_type('energy grid', energy, Iterable, Real) self._energy = energy + @fission_energy.setter + def fission_energy(self, fission_energy): + cv.check_type('fission energy release', fission_energy, + FissionEnergyRelease) + self._fission_energy = fission_energy + @reactions.setter def reactions(self, reactions): cv.check_type('reactions', reactions, Mapping) @@ -276,6 +288,11 @@ class IncidentNeutron(object): urr_group = g.create_group('urr') self.urr.to_hdf5(urr_group) + # Write fission energy release data + if self.fission_energy is not None: + fer_group = g.create_group('fission_energy_release') + self.fission_energy.to_hdf5(fer_group) + f.close() @classmethod @@ -331,6 +348,11 @@ class IncidentNeutron(object): urr_group = group['urr'] data.urr = ProbabilityTables.from_hdf5(urr_group) + # Read fission energy release data + if 'fission_energy_release' in group: + fer_group = group['fission_energy_release'] + data.fission_energy = FissionEnergyRelease.from_hdf5(fer_group) + return data @classmethod From 9acac0725a26efc14906d692bc99f0dacabeccea Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 29 Jul 2016 13:54:30 -0500 Subject: [PATCH 02/33] Add Madland fission-Q support to F90 --- openmc/data/fission_energy.py | 27 ++++++ src/constants.F90 | 6 +- src/input_xml.F90 | 4 + src/nuclide_header.F90 | 34 ++++++- src/output.F90 | 2 + src/tally.F90 | 162 +++++++++++++++++++++++++++++----- 6 files changed, 209 insertions(+), 26 deletions(-) diff --git a/openmc/data/fission_energy.py b/openmc/data/fission_energy.py index 5716a3e85..7415a9b51 100644 --- a/openmc/data/fission_energy.py +++ b/openmc/data/fission_energy.py @@ -62,6 +62,19 @@ class FissionEnergyRelease(object): return Sum([self.fragments, self.prompt_neutrons, self.delayed_neutrons, self.prompt_photons, self.delayed_photons, self.betas, self.neutrinos]) + + @property + def prompt_q(self): + return Sum([self.fragments, self.prompt_neutrons, self.prompt_photons, + lambda E: -E]) + + @property + def recoverable_q(self): + return Sum([self.recoverable, lambda E: -E]) + + @property + def total_q(self): + return Sum([self.total, lambda E: -E]) @property def form(self): @@ -272,6 +285,20 @@ class FissionEnergyRelease(object): data=self.delayed_photons.coef) group.create_dataset('betas', data=self.betas.coef) group.create_dataset('neutrinos', data=self.neutrinos.coef) + + q_prompt = (self.fragments + self.prompt_neutrons + + self.prompt_photons + Polynomial((-1.0, 0.0))) + group.create_dataset('q_prompt', data=q_prompt.coef) + q_recoverable = (self.fragments + self.prompt_neutrons + + self.delayed_neutrons + self.prompt_photons + + self.delayed_photons + self.betas + + Polynomial((-1.0, 0.0))) + group.create_dataset('q_recoverable', data=q_recoverable.coef) + q_total = (self.fragments + self.prompt_neutrons + + self.delayed_neutrons + self.prompt_photons + + self.delayed_photons + self.betas + self.neutrinos + + Polynomial((-1.0, 0.0))) + group.create_dataset('q_total', data=q_total.coef) elif self.form == 'Sher-Beck': group.attrs['format'] = np.string_('Sher-Beck') self.fragments.to_hdf5(group, 'fragments') diff --git a/src/constants.F90 b/src/constants.F90 index d2bb6e1cf..5447b5cc2 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -289,7 +289,7 @@ module constants EVENT_ABSORB = 2 ! Tally score type - integer, parameter :: N_SCORE_TYPES = 20 + integer, parameter :: N_SCORE_TYPES = 22 integer, parameter :: & SCORE_FLUX = -1, & ! flux SCORE_TOTAL = -2, & ! total reaction rate @@ -310,7 +310,9 @@ module constants SCORE_NU_SCATTER_YN = -17, & ! angular flux-weighted nu-scattering moment (0:N) SCORE_EVENTS = -18, & ! number of events SCORE_DELAYED_NU_FISSION = -19, & ! delayed neutron production rate - SCORE_INVERSE_VELOCITY = -20 ! flux-weighted inverse velocity + SCORE_INVERSE_VELOCITY = -20, & ! flux-weighted inverse velocity + SCORE_FISS_Q_PROMPT = -21, & ! prompt fission Q-value + SCORE_FISS_Q_RECOV = -22 ! recoverable fission Q-value ! Maximum scattering order supported integer, parameter :: MAX_ANG_ORDER = 10 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 3c8291e9c..f661927d7 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3629,6 +3629,10 @@ contains t % score_bins(j) = SCORE_KAPPA_FISSION case ('inverse-velocity') t % score_bins(j) = SCORE_INVERSE_VELOCITY + case ('fission-q-prompt') + t % score_bins(j) = SCORE_FISS_Q_PROMPT + case ('fission-q-recoverable') + t % score_bins(j) = SCORE_FISS_Q_RECOV case ('current') t % score_bins(j) = SCORE_CURRENT t % type = TALLY_SURFACE_CURRENT diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 8ff83482e..31e069b68 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -88,6 +88,10 @@ module nuclide_header type(DictIntInt) :: reaction_index ! map MT values to index in reactions ! array; used at tally-time + ! Fission energy release + class(Function1D), allocatable :: fission_q_prompt ! prompt neutrons, gammas + class(Function1D), allocatable :: fission_q_recov ! neutrons, gammas, betas + contains procedure :: clear => nuclide_clear procedure :: print => nuclide_print @@ -192,6 +196,8 @@ module nuclide_header integer(HID_T) :: rxs_group integer(HID_T) :: rx_group integer(HID_T) :: total_nu + integer(HID_T) :: fer_group ! fission_energy_release group + integer(HID_T) :: fer_dset integer(SIZE_T) :: name_len, name_file_len integer(HSIZE_T) :: j integer(HSIZE_T) :: dims(1) @@ -251,8 +257,8 @@ module nuclide_header call this % urr_data % from_hdf5(urr_group) ! if the inelastic competition flag indicates that the inelastic cross - ! section should be determined from a normal reaction cross section, we need - ! to get the index of the reaction + ! section should be determined from a normal reaction cross section, we + ! need to get the index of the reaction if (this % urr_data % inelastic_flag > 0) then do i = 1, size(this % reactions) if (this % reactions(i) % MT == this % urr_data % inelastic_flag) then @@ -296,6 +302,30 @@ module nuclide_header call close_group(nu_group) end if + ! Read fission energy release data if present + call h5ltpath_valid_f(group_id, 'fission_energy_release', .true., exists, & + hdf5_err) + if (exists) then + fer_group = open_group(group_id, 'fission_energy_release') + call read_attribute(temp, fer_group, 'format') + if (temp == 'Madland') then + ! The data uses the Madland format, i.e. polynomials + + ! Read the prompt Q-value + allocate(Polynomial :: this % fission_q_prompt) + fer_dset = open_dataset(fer_group, 'q_prompt') + call this % fission_q_prompt % from_hdf5(fer_dset) + call close_dataset(fer_dset) + + ! Read the recoverable energy Q-value + allocate(Polynomial :: this % fission_q_recov) + fer_dset = open_dataset(fer_group, 'q_recoverable') + call this % fission_q_recov % from_hdf5(fer_dset) + call close_dataset(fer_dset) + end if + call close_group(fer_group) + end if + ! Create derived cross section data call this % create_derived() diff --git a/src/output.F90 b/src/output.F90 index 09df23f0d..35524a9a3 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -791,6 +791,8 @@ contains score_names(abs(SCORE_NU_SCATTER_YN)) = "Scattering Prod. Rate Moment" score_names(abs(SCORE_DELAYED_NU_FISSION)) = "Delayed-Nu-Fission Rate" score_names(abs(SCORE_INVERSE_VELOCITY)) = "Flux-Weighted Inverse Velocity" + score_names(abs(SCORE_FISS_Q_PROMPT)) = "Prompt fission power" + score_names(abs(SCORE_FISS_Q_RECOV)) = "Recoverable fission power" ! Create filename for tally output filename = trim(path_output) // "tallies.out" diff --git a/src/tally.F90 b/src/tally.F90 index b40a39ce0..1271da7da 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -625,14 +625,14 @@ contains if (survival_biasing) then ! No fission events occur if survival biasing is on -- need to ! calculate fraction of absorptions that would have resulted in - ! fission scale by kappa-fission - associate (nuc => nuclides(p%event_nuclide)) - if (micro_xs(p%event_nuclide)%absorption > ZERO .and. & - nuc%fissionable) then - score = p%absorb_wgt * & - nuc%reactions(nuc%index_fission(1))%Q_value * & - micro_xs(p%event_nuclide)%fission / & - micro_xs(p%event_nuclide)%absorption + ! fission scaled by kappa-fission + associate (nuc => nuclides(p % event_nuclide)) + if (micro_xs(p % event_nuclide) % absorption > ZERO .and. & + nuc % fissionable) then + score = p % absorb_wgt * & + nuc % reactions(nuc % index_fission(1)) % Q_value * & + micro_xs(p % event_nuclide) % fission / & + micro_xs(p % event_nuclide) % absorption end if end associate else @@ -641,12 +641,12 @@ contains ! All fission events will contribute, so again we can use ! particle's weight entering the collision as the estimate for ! the fission energy production rate - associate (nuc => nuclides(p%event_nuclide)) - if (nuc%fissionable) then - score = p%last_wgt * & - nuc%reactions(nuc%index_fission(1))%Q_value * & - micro_xs(p%event_nuclide)%fission / & - micro_xs(p%event_nuclide)%absorption + associate (nuc => nuclides(p % event_nuclide)) + if (nuc % fissionable) then + score = p % last_wgt * & + nuc % reactions(nuc % index_fission(1)) % Q_value * & + micro_xs(p % event_nuclide) % fission / & + micro_xs(p % event_nuclide) % absorption end if end associate end if @@ -654,22 +654,23 @@ contains else if (i_nuclide > 0) then associate (nuc => nuclides(i_nuclide)) - if (nuc%fissionable) then - score = nuc%reactions(nuc%index_fission(1))%Q_value * & - micro_xs(i_nuclide)%fission * atom_density * flux + if (nuc % fissionable) then + score = nuc % reactions(nuc % index_fission(1)) % Q_value * & + micro_xs(i_nuclide) % fission * atom_density * flux end if end associate else - do l = 1, materials(p%material)%n_nuclides + do l = 1, materials(p % material) % n_nuclides ! Determine atom density and index of nuclide - atom_density_ = materials(p%material)%atom_density(l) - i_nuc = materials(p%material)%nuclide(l) + atom_density_ = materials(p % material) % atom_density(l) + i_nuc = materials(p % material) % nuclide(l) ! If nuclide is fissionable, accumulate kappa fission associate(nuc => nuclides(i_nuc)) if (nuc % fissionable) then - score = score + nuc%reactions(nuc%index_fission(1))%Q_value * & - micro_xs(i_nuc)%fission * atom_density_ * flux + score = score + & + nuc % reactions(nuc % index_fission(1)) % Q_value * & + micro_xs(i_nuc) % fission * atom_density_ * flux end if end associate end do @@ -694,6 +695,123 @@ contains end if end if + case (SCORE_FISS_Q_PROMPT) + if (t % estimator == ESTIMATOR_ANALOG) then + if (survival_biasing) then + ! No fission events occur if survival biasing is on -- need to + ! calculate fraction of absorptions that would have resulted in + ! fission scaled by Q-value + associate (nuc => nuclides(p % event_nuclide)) + if (micro_xs(p % event_nuclide) % absorption > ZERO .and. & + allocated(nuc % fission_q_prompt)) then + score = p % absorb_wgt & + * nuc % fission_q_prompt % evaluate(p % last_E) & + * micro_xs(p % event_nuclide) % fission & + / micro_xs(p % event_nuclide) % absorption + end if + end associate + else + ! Skip any non-absorption events + if (p % event == EVENT_SCATTER) cycle SCORE_LOOP + ! All fission events will contribute, so again we can use + ! particle's weight entering the collision as the estimate for + ! the fission energy production rate + associate (nuc => nuclides(p % event_nuclide)) + if (allocated(nuc % fission_q_prompt)) then + score = p % last_wgt & + * nuc % fission_q_prompt % evaluate(p % last_E) & + * micro_xs(p % event_nuclide) % fission & + / micro_xs(p % event_nuclide) % absorption + end if + end associate + end if + + else + if (t % estimator == ESTIMATOR_COLLISION) then + E = p % last_E + else + E = p % E + end if + + if (i_nuclide > 0) then + if (allocated(nuclides(i_nuclide) % fission_q_prompt)) then + score = micro_xs(i_nuclide) % fission * atom_density * flux & + * nuclides(i_nuclide) % fission_q_prompt % evaluate(E) + else + score = ZERO + end if + else + score = ZERO + do l = 1, materials(p % material) % n_nuclides + atom_density_ = materials(p % material) % atom_density(l) + i_nuc = materials(p % material) % nuclide(l) + if (allocated(nuclides(i_nuc) % fission_q_prompt)) then + score = score + micro_xs(i_nuc) % fission * atom_density_ & + * flux & + * nuclides(i_nuc) % fission_q_prompt % evaluate(E) + end if + end do + end if + end if + + case (SCORE_FISS_Q_RECOV) + if (t % estimator == ESTIMATOR_ANALOG) then + if (survival_biasing) then + ! No fission events occur if survival biasing is on -- need to + ! calculate fraction of absorptions that would have resulted in + ! fission scaled by Q-value + associate (nuc => nuclides(p % event_nuclide)) + if (micro_xs(p % event_nuclide) % absorption > ZERO .and. & + allocated(nuc % fission_q_recov)) then + score = p % absorb_wgt & + * nuc % fission_q_recov % evaluate(p % last_E) & + * micro_xs(p % event_nuclide) % fission & + / micro_xs(p % event_nuclide) % absorption + end if + end associate + else + ! Skip any non-absorption events + if (p % event == EVENT_SCATTER) cycle SCORE_LOOP + ! All fission events will contribute, so again we can use + ! particle's weight entering the collision as the estimate for + ! the fission energy production rate + associate (nuc => nuclides(p % event_nuclide)) + if (allocated(nuc % fission_q_recov)) then + score = p % last_wgt & + * nuc % fission_q_recov % evaluate(p % last_E) & + * micro_xs(p % event_nuclide) % fission & + / micro_xs(p % event_nuclide) % absorption + end if + end associate + end if + + else + if (t % estimator == ESTIMATOR_COLLISION) then + E = p % last_E + else + E = p % E + end if + + if (i_nuclide > 0) then + if (allocated(nuclides(i_nuclide) % fission_q_recov)) then + score = micro_xs(i_nuclide) % fission * atom_density * flux & + * nuclides(i_nuclide) % fission_q_recov % evaluate(E) + else + score = ZERO + end if + else + score = ZERO + do l = 1, materials(p % material) % n_nuclides + atom_density_ = materials(p % material) % atom_density(l) + i_nuc = materials(p % material) % nuclide(l) + if (allocated(nuclides(i_nuc) % fission_q_recov)) then + score = score + micro_xs(i_nuc) % fission * atom_density_ & + * flux * nuclides(i_nuc) % fission_q_recov % evaluate(E) + end if + end do + end if + end if + case default if (t % estimator == ESTIMATOR_ANALOG) then ! Any other score is assumed to be a MT number. Thus, we just need From b64dcd24c83d05e6041f01f84e65b0d0355d97fa Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 29 Jul 2016 16:09:06 -0500 Subject: [PATCH 03/33] Add Sher-Beck fission-Q support --- openmc/data/fission_energy.py | 96 ++++++++++++++++++++++++----------- src/nuclide_header.F90 | 17 +++++++ src/tally.F90 | 4 +- 3 files changed, 84 insertions(+), 33 deletions(-) diff --git a/openmc/data/fission_energy.py b/openmc/data/fission_energy.py index 7415a9b51..b7641a237 100644 --- a/openmc/data/fission_energy.py +++ b/openmc/data/fission_energy.py @@ -1,4 +1,5 @@ from collections import Callable +from copy import deepcopy import sys #from warnings import warn @@ -64,16 +65,16 @@ class FissionEnergyRelease(object): self.neutrinos]) @property - def prompt_q(self): + def q_prompt(self): return Sum([self.fragments, self.prompt_neutrons, self.prompt_photons, lambda E: -E]) @property - def recoverable_q(self): + def q_recoverable(self): return Sum([self.recoverable, lambda E: -E]) @property - def total_q(self): + def q_total(self): return Sum([self.total, lambda E: -E]) @property @@ -226,7 +227,37 @@ class FissionEnergyRelease(object): out.neutrinos = Polynomial(value['ENU']) else: out.form = 'Sher-Beck' - raise NotImplemented + + # EFR and ENP are energy independent. Polynomial is used because it + # has a __call__ attribute that handles Iterable inputs. The + # energy-dependence of END is unspecified in ENDF-102 so assume it + # is independent. + out.fragments = Polynomial((value['EFR'][0])) + out.prompt_photons = Polynomial((value['EGP'][0])) + out.delayed_neutrons = Polynomial((value['END'][0])) + + # EDP, EB, and ENU are linear. + out.delayed_photons = Polynomial((value['EGD'][0], -0.075)) + out.betas = Polynomial((value['EB'][0], -0.075)) + out.neutrinos = Polynomial((value['ENU'][0], -0.105)) + + # Prompt neutrons require nu-data. It is not clear from ENDF-102 + # whether prompt or total nu values should be used, but the delayed + # neutron fraction is so small that the difference is negligible. + nu_prompt = [p for p in incident_neutron[18].products + if p.particle == 'neutron' + and p.emission_mode == 'prompt'] + if len(nu_prompt) == 0: + raise ValueError('Nu data is needed to compute fission energy ' + 'release with the Sher-Beck format.') + if len(nu_prompt) > 1: + raise ValueError('Ambiguous prompt nu value.') + if not isinstance(nu_prompt[0].yield_, Tabulated1D): + raise TypeError('Sher-Beck fission energy release currently ' + 'only supports Tabulated1D nu data.') + ENP = deepcopy(nu_prompt[0].yield_) + ENP.y = value['ENP'] + 1.307 * ENP.x - 8.07 * (ENP.y - ENP.y[0]) + out.prompt_neutrons = ENP return out @@ -247,16 +278,19 @@ class FissionEnergyRelease(object): """ obj = cls() + + obj.fragments = Polynomial(group['fragments'].value) + obj.delayed_neutrons = Polynomial(group['delayed_neutrons'].value) + obj.prompt_photons = Polynomial(group['prompt_photons'].value) + obj.delayed_photons = Polynomial(group['delayed_photons'].value) + obj.betas = Polynomial(group['betas'].value) + obj.neutrinos = Polynomial(group['neutrinos'].value) + if group.attrs['format'] == 'Madland': - obj.fragments = Polynomial(group['fragments'].value) obj.prompt_neutrons = Polynomial(group['prompt_neutrons'].value) - obj.delayed_neutrons = Polynomial(group['delayed_neutrons'].value) - obj.prompt_photons = Polynomial(group['prompt_photons'].value) - obj.delayed_photons = Polynomial(group['delayed_photons'].value) - obj.betas = Polynomial(group['betas'].value) - obj.neutrinos = Polynomial(group['neutrinos'].value) elif group.attrs['format'] == 'Sher-Beck': - raise NotImplemented + obj.prompt_neutrons = Tabulated1D.from_hdf5( + group['prompt_neutrons']) else: raise ValueError('Unrecognized energy release format') @@ -272,19 +306,20 @@ class FissionEnergyRelease(object): """ + group.create_dataset('fragments', data=self.fragments.coef) + group.create_dataset('delayed_neutrons', + data=self.delayed_neutrons.coef) + group.create_dataset('prompt_photons', + data=self.prompt_photons.coef) + group.create_dataset('delayed_photons', + data=self.delayed_photons.coef) + group.create_dataset('betas', data=self.betas.coef) + group.create_dataset('neutrinos', data=self.neutrinos.coef) + if self.form == 'Madland': group.attrs['format'] = np.string_('Madland') - group.create_dataset('fragments', data=self.fragments.coef) group.create_dataset('prompt_neutrons', data=self.prompt_neutrons.coef) - group.create_dataset('delayed_neutrons', - data=self.delayed_neutrons.coef) - group.create_dataset('prompt_photons', - data=self.prompt_photons.coef) - group.create_dataset('delayed_photons', - data=self.delayed_photons.coef) - group.create_dataset('betas', data=self.betas.coef) - group.create_dataset('neutrinos', data=self.neutrinos.coef) q_prompt = (self.fragments + self.prompt_neutrons + self.prompt_photons + Polynomial((-1.0, 0.0))) @@ -294,19 +329,18 @@ class FissionEnergyRelease(object): self.delayed_photons + self.betas + Polynomial((-1.0, 0.0))) group.create_dataset('q_recoverable', data=q_recoverable.coef) - q_total = (self.fragments + self.prompt_neutrons + - self.delayed_neutrons + self.prompt_photons + - self.delayed_photons + self.betas + self.neutrinos + - Polynomial((-1.0, 0.0))) - group.create_dataset('q_total', data=q_total.coef) elif self.form == 'Sher-Beck': group.attrs['format'] = np.string_('Sher-Beck') - self.fragments.to_hdf5(group, 'fragments') self.prompt_neutrons.to_hdf5(group, 'prompt_neutrons') - self.delayed_neutrons.to_hdf5(group, 'delayed_neutrons') - self.prompt_photons.to_hdf5(group, 'prompt_photons') - self.delayed_photons.to_hdf5(group, 'delayed_photons') - self.betas.to_hdf5(group, 'betas') - self.neutrinos.to_hdf5(group, 'neutrinos') + + q_prompt = deepcopy(self.prompt_neutrons) + q_prompt.y += self.fragments(q_prompt.x) + q_prompt.y += self.prompt_photons(q_prompt.x) + q_prompt.to_hdf5(group, 'q_prompt') + q_recoverable = q_prompt + q_recoverable.y += self.delayed_neutrons(q_recoverable.x) + q_recoverable.y += self.delayed_photons(q_recoverable.x) + q_recoverable.y += self.betas(q_recoverable.x) + q_recoverable.to_hdf5(group, 'q_recoverable') else: raise ValueError('Unrecognized energy release format') diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 31e069b68..2a2a376bd 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -322,6 +322,23 @@ module nuclide_header fer_dset = open_dataset(fer_group, 'q_recoverable') call this % fission_q_recov % from_hdf5(fer_dset) call close_dataset(fer_dset) + else if (temp == 'Sher-Beck') then + ! The data uses the Sher-Beck format. Python has handily converted this + ! format to Tabulated1Ds. + + ! Read the prompt Q-value + allocate(Tabulated1D :: this % fission_q_prompt) + fer_dset = open_dataset(fer_group, 'q_prompt') + call this % fission_q_prompt % from_hdf5(fer_dset) + call close_dataset(fer_dset) + + ! Read the recoverable energy Q-value + allocate(Tabulated1D :: this % fission_q_recov) + fer_dset = open_dataset(fer_group, 'q_recoverable') + call this % fission_q_recov % from_hdf5(fer_dset) + call close_dataset(fer_dset) + else + call fatal_error('Unrecognized fission energy release format.') end if call close_group(fer_group) end if diff --git a/src/tally.F90 b/src/tally.F90 index 1271da7da..8d14b1c2a 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -703,7 +703,7 @@ contains ! fission scaled by Q-value associate (nuc => nuclides(p % event_nuclide)) if (micro_xs(p % event_nuclide) % absorption > ZERO .and. & - allocated(nuc % fission_q_prompt)) then + allocated(nuc % fission_q_prompt)) then score = p % absorb_wgt & * nuc % fission_q_prompt % evaluate(p % last_E) & * micro_xs(p % event_nuclide) % fission & @@ -762,7 +762,7 @@ contains ! fission scaled by Q-value associate (nuc => nuclides(p % event_nuclide)) if (micro_xs(p % event_nuclide) % absorption > ZERO .and. & - allocated(nuc % fission_q_recov)) then + allocated(nuc % fission_q_recov)) then score = p % absorb_wgt & * nuc % fission_q_recov % evaluate(p % last_E) & * micro_xs(p % event_nuclide) % fission & From 725ef2d9238293c0863297ba9e987ed666f932c6 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 1 Aug 2016 22:34:04 -0500 Subject: [PATCH 04/33] Add function to create compact MT=458 data library --- openmc/data/endf_utils.py | 2 + openmc/data/fission_energy.py | 255 ++++++++++++++++++++++++++-------- 2 files changed, 196 insertions(+), 61 deletions(-) diff --git a/openmc/data/endf_utils.py b/openmc/data/endf_utils.py index 6db3c611c..bfd9ae5c8 100644 --- a/openmc/data/endf_utils.py +++ b/openmc/data/endf_utils.py @@ -24,6 +24,7 @@ def read_CONT_line(line): int(line[66:70]), int(line[70:72]), int(line[72:75]), int(line[75:80])) + def identify_nuclide(fname): """Read the header of an ENDF file and extract identifying information.""" with open(fname, 'r') as fh: @@ -39,5 +40,6 @@ def identify_nuclide(fname): # Return dictionary of the most important identifying information. return {'Z': int(ZA) // 1000, 'A': int(ZA) % 1000, + 'LFI': bool(LFI), 'LIS': LIS, 'LISO': LISO} diff --git a/openmc/data/fission_energy.py b/openmc/data/fission_energy.py index b7641a237..ccea3d96a 100644 --- a/openmc/data/fission_energy.py +++ b/openmc/data/fission_energy.py @@ -3,6 +3,7 @@ from copy import deepcopy import sys #from warnings import warn +import h5py import numpy as np from numpy.polynomial.polynomial import Polynomial @@ -14,6 +15,189 @@ if sys.version_info[0] >= 3: basestring = str +def _extract_458_data(filename): + """Read an ENDF file and extract the MF=1, MT=458 values. + + Parameters + ---------- + filename : str + Path to and ENDF file + + Returns + ------- + value : dict of str to list of float + Dictionary that gives lists of coefficients for each energy component. + The keys are the 2-3 letter strings used in ENDF-102, e.g. 'EFR' and + 'ET'. The list will have a length of 1 for Sher-Beck data, more for + polynomial data. + uncertainty : dict of str to list of float + A dictionary with the same format as above. This is probably a + one-standard deviation value, but that is not specified explicitly in + ENDF-102. Also, some evaluations will give zero uncertainty. Use with + caution. + + """ + ident = identify_nuclide(filename) + + if not ident['LFI']: + # This nuclide isn't fissionable. + return None + + # Extract the MF=1, MT=458 section. + lines = [] + with open(filename, 'r') as fh: + line = fh.readline() + while line != '': + if line[70:75] == ' 1458': + lines.append(line) + line = fh.readline() + + if len(lines) == 0: + # No 458 data here. + return None + + # Read the number of coefficients in this LIST record. + NPL = read_CONT_line(lines[1])[4] + + # Parse the ENDF LIST into an array. + data = [] + for i in range(NPL): + row, column = divmod(i, 6) + data.append(read_float(lines[2 + row][11*column:11*(column+1)])) + + # Declare the coefficient names and the order they are given in. The LIST + # contains a value followed immediately by an uncertainty for each of these + # components, times the polynomial order + 1. + labels = ('EFR', 'ENP', 'END', 'EGP', 'EGD', 'EB', 'ENU', 'ER', 'ET') + + # Associate each set of values and uncertainties with its label. + value = dict() + uncertainty = dict() + for i in range(len(labels)): + value[labels[i]] = data[2*i::18] + uncertainty[labels[i]] = data[2*i + 1::18] + + # In ENDF/B-7.1, data for 2nd-order coefficients were mistakenly not + # converted from MeV to eV. Check for this error and fix it if present. + n_coeffs = len(value['EFR']) + if n_coeffs == 3: # Only check 2nd-order data. + # Check each energy component for the error. If a 1 MeV neutron + # causes a change of more than 100 MeV, we know something is wrong. + error_present = False + for coeffs in value.values(): + second_order = coeffs[2] + if abs(second_order) * 1e12 > 1e8: + error_present = True + break + + # If we found the error, reduce all 2nd-order coeffs by 10**6. + if error_present: + for coeffs in value.values(): coeffs[2] *= 1e-6 + for coeffs in uncertainty.values(): coeffs[2] *= 1e-6 + + # Perform the sanity check again... just in case. + for coeffs in value.values(): + second_order = coeffs[2] + if abs(second_order) * 1e12 > 1e8: + raise ValueError("Encountered a ludicrously large second-" + "order polynomial coefficient.") + + # Convert eV to MeV. + for coeffs in value.values(): + for i in range(len(coeffs)): + coeffs[i] *= 10**(-6 + 6*i) + for coeffs in uncertainty.values(): + for i in range(len(coeffs)): + coeffs[i] *= 10**(-6 + 6*i) + + return value, uncertainty + + +def write_compact_458_library(endf_files, output_name=None, comment=None, + verbose=False): + """Read ENDF files, strip the MF=1 MT=458 data and write to small HDF5. + + Parameters + ---------- + endf_files : Collection of str + Strings giving the paths to the ENDF files that will be parsed for data. + output_name : str + Name of the output HDF5 file. Default is 'fission_Q_data.h5'. + comment : str + Comment to write in the output HDF5 file. Defaults to no comment. + verbose : bool + If True, print the name of each isomer as it is read. Defaults to + False. + + """ + # Open the output file. + if output_name is None: output_name = 'fission_Q_data.h5' + out = h5py.File(output_name, 'w', libver='latest') + + # Write comments, if given. This commented out comment is the one used for + # the library distributed with OpenMC. + #comment = ('This data is extracted from ENDF/B-VII.1 library. Thanks ' + # 'evaluators, for all your hard work :) Citation: ' + # 'M. B. Chadwick, M. Herman, P. Oblozinsky, ' + # 'M. E. Dunn, Y. Danon, A. C. Kahler, D. L. Smith, ' + # 'B. Pritychenko, G. Arbanas, R. Arcilla, R. Brewer, ' + # 'D. A. Brown, R. Capote, A. D. Carlson, Y. S. Cho, H. Derrien, ' + # 'K. Guber, G. M. Hale, S. Hoblit, S. Holloway, T. D. Johnson, ' + # 'T. Kawano, B. C. Kiedrowski, H. Kim, S. Kunieda, ' + # 'N. M. Larson, L. Leal, J. P. Lestone, R. C. Little, ' + # 'E. A. McCutchan, R. E. MacFarlane, M. MacInnes, ' + # 'C. M. Mattoon, R. D. McKnight, S. F. Mughabghab, ' + # 'G. P. A. Nobre, G. Palmiotti, A. Palumbo, M. T. Pigni, ' + # 'V. G. Pronyaev, R. O. Sayer, A. A. Sonzogni, N. C. Summers, ' + # 'P. Talou, I. J. Thompson, A. Trkov, R. L. Vogt, ' + # 'S. C. van der Marck, A. Wallner, M. C. White, D. Wiarda, ' + # 'and P. G. Young. ENDF/B-VII.1 nuclear data for science and ' + # 'technology: Cross sections, covariances, fission product ' + # 'yields and decay data", Nuclear Data Sheets, ' + # '112(12):2887-2996 (2011).') + if comment is not None: + out.attrs['comment'] = np.string_(comment) + + # Declare the order of the components. Use fixed-length numpy strings + # because they work well with h5py. + labels = np.array(('EFR', 'ENP', 'END', 'EGP', 'EGD', 'EB', 'ENU', 'ER', + 'ET'), dtype='S3') + out.attrs['component order'] = labels + + # Iterate over the given files. + if verbose: print('Reading ENDF files:') + for fname in endf_files: + if verbose: print(fname) + + ident = identify_nuclide(fname) + + # Skip non-fissionable nuclides. + if not ident['LFI']: continue + + # Get the important bits. + data = _extract_458_data(fname) + if data is None: continue + value, uncertainty = data + + # Make a group for this isomer. + name = str(ident['Z']) + str(ident['A']) + if ident['LISO'] != 0: + name += '_m' + str(ident['LISO']) + nuclide_group = out.create_group(name) + + # Write all the coefficients into one array. The first dimension gives + # the component (e.g. fragments or prompt neutrons); the second switches + # between value and uncertainty; the third gives the polynomial order. + n_coeffs = len(value['EFR']) + data_out = np.zeros((len(labels), 2, n_coeffs)) + for i, label in enumerate(labels): + data_out[i, 0, :] = value[label.decode()] + data_out[i, 1, :] = uncertainty[label.decode()] + nuclide_group.create_dataset('data', data=data_out) + + out.close() + + class FissionEnergyRelease(object): def __init__(self): self._fragments = None @@ -148,72 +332,21 @@ class FissionEnergyRelease(object): pass if ident['LISO'] != incident_neutron.metastable: pass + if not ident['LIF']: + pass - # Extract the MF=1, MT=458 section. - lines = [] - with open(filename, 'r') as fh: - line = fh.readline() - while line != '': - if line[70:75] == ' 1458': - lines.append(line) - line = fh.readline() + # Read the 458 data from the ENDF file. + value, uncertainty = _extract_458_data(filename) - # Read the number of coefficients in this LIST record. - NPL = read_CONT_line(lines[1])[4] - - # Parse the ENDF LIST into an array. - data = [] - for i in range(NPL): - row, column = divmod(i, 6) - data.append(read_float(lines[2 + row][11*column:11*(column+1)])) - - # Declare the coefficient names and the order they are given in. The - # LIST contains a value followed immediately by an uncertainty for each - # of these components, times the polynomial order + 1. If we only find - # one value for each of these components, then we need to use the - # Sher-Beck formula for energy dependence. Otherwise, it is a - # polynomial. + # Declare the coefficient names. If we only find one value for each of + # these components, then we need to use the Sher-Beck formula for energy + # dependence. Otherwise, it is a polynomial. labels = ('EFR', 'ENP', 'END', 'EGP', 'EGD', 'EB', 'ENU', 'ER', 'ET') - # Associate each set of values and uncertainties with its label. - value = dict() - uncertainty = dict() - for i in range(len(labels)): - value[labels[i]] = data[2*i::18] - uncertainty[labels[i]] = data[2*i + 1::18] - - # In ENDF/B-7.1, data for 2nd-order coefficients were mistakenly not - # converted from MeV to eV. Check for this error and fix it if present. + # How many coefficients are given for each coefficient? If we only find + # one value for each, then we need to use the Sher-Beck formula for + # energy dependence. Otherwise, it is a polynomial. n_coeffs = len(value['EFR']) - if n_coeffs == 3: # Only check 2nd-order data. - # Check each energy component for the error. If a 1 MeV neutron - # causes a change of more than 100 MeV, we know something is wrong. - error_present = False - for coeffs in value.values(): - second_order = coeffs[2] - if abs(second_order) * 1e12 > 1e8: - error_present = True - break - - # If we found the error, reduce all 2nd-order coeffs by 10**6. - if error_present: - for coeffs in value.values(): coeffs[2] *= 1e-6 - for coeffs in uncertainty.values(): coeffs[2] *= 1e-6 - - # Perform the sanity check again... just in case. - for coeffs in value.values(): - second_order = coeffs[2] - if abs(second_order) * 1e12 > 1e8: - raise ValueError("Encountered a ludicrously large second-" - "order polynomial coefficient.") - - # Convert eV to MeV. - for coeffs in value.values(): - for i in range(len(coeffs)): - coeffs[i] *= 10**(-6 + 6*i) - for coeffs in uncertainty.values(): - for i in range(len(coeffs)): - coeffs[i] *= 10**(-6 + 6*i) out = cls() if n_coeffs > 1: From 84cfd2d6a940920473283de15690160d6de59b58 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Tue, 2 Aug 2016 14:10:32 -0500 Subject: [PATCH 05/33] Add MT=458 data library and supporting methods --- data/fission_Q_data_endfb71.h5 | Bin 0 -> 67543 bytes data/get_nndc_data.py | 3 +- openmc/data/fission_energy.py | 195 ++++++++++++++++++++++----------- openmc/data/neutron.py | 3 + scripts/openmc-ace-to-hdf5 | 17 +++ 5 files changed, 155 insertions(+), 63 deletions(-) create mode 100644 data/fission_Q_data_endfb71.h5 diff --git a/data/fission_Q_data_endfb71.h5 b/data/fission_Q_data_endfb71.h5 new file mode 100644 index 0000000000000000000000000000000000000000..89e51b70783034b40ce4e207fdb3e34b5d5d33cc GIT binary patch literal 67543 zcmeHw2|QI@_y3VZloV1Vm&_z%(Qx-BWT?oPN}1;&V`)Uuq?FRAQX$QyXym3rrJ_Wl z3=NtmGmZZ<+;h&u`Memf&O(pgizCC` zJtMunSU%yQ!7Q_xrsm^J#yVMAYU!|o{5`_m!(+5qEC)aLkU*3;+C3=BJu);rLX+ha z8qRVL3Sz~CMuoHd+{3+Ci$cQ#SrgTSg{S&Qx<~qlhD;Qmp~W)MVomjP_gdud8K}uZ z%NE|@!R{fNEL$y>jYm-E6915hz!>}+GcA^BR0#T&3wq)n5{jM~qc>}@tlj;Byu&qF zrdq7&S}glu|42VgjIV9De`JiOpLa-Ls3vQg7Rxx?!#%_ut(zTo?&%*CkH`-0;t#~;yG8Ah!$_%aY z3{UG2e_ubMJ?3aR%Gb}`1N(BW&G+g8f4yBhhA{*U?f`um{>m zlrP$3TYujW^b04n=g>lUXh@8^cQpPW8?=wzW3bLd??V5#4-HuoiW8#^N9$!Dg^pWn zn4!;ea1RQN(qvhp-H8_bLW39LeS+TQ5FQvx@IyPuDb$y!G*y z7#`{sZKAHB;e@fe zMn>aVs=C@bI%-;jY!Ou8#ZZBlzy*FGIy6EeS)t)*58#52%01GLkB@I4_w7L*&CKo0 zX4;}6Z#K==Y?`T=31T+W(aa7pb3l+GTw40@p0le3e@C&f_tRhhE$%lABO_lx;pxnt1Oc$Yz|lPyH#maiX|)052BgHi%aMYkaSWHn+DoLyWkX@}WcMI|c}!2!ga? zoyk{>Amg!tm9@fVI`;TK0dM_K?8ykym*JJ(X`|OWc%4ZJ;1#@cJ2Z%X>%reY0q69EtncJ9t|}c?>Vi z_Kf{c!%IqOZ=^R?il`+ewA4vYYyY34XFnhJ;P1GLHyM}uocDz9z51VOY&pjv=^Z6b z<>a}@!~SIzYn8*o0eL;>&G#dXWg>8(%)a?*cW;B6-^0E)Y<1*FkE>A}5SIzQtq8ah zxpghj`KYb@Gv5?w?=znIyZ;WbcspxM`ZZ%te#&mmjGQ;@4fQKEs?HdL)fa2-E)Q1( zR`skvXWi#)k*Tj*T1Sn7=|3%oj&#~WJ$KT}`GZA&odG@*)jwU=R;I37WKBK%re-kI z!#`)K`u{$9cxCOTt>D7wgb%YWtf1%}hgAu>8Hos zV-?`w)nQAn@3RA5_gt-S@ArY%5)V5UH#Ji9JBL}e=KI9 z4?g+&t~vRwi9PYe)ATlL83tZSc|$*1FJ|Jk;`_vfB|^;px0^b7Q`=D{UgwQJosWLO zq(8qv*4n~|iI+1xu#bI>i5JPgJ6^#~%TMMtdBc8tE%zK<+Q_*lUU1^*jBqfj6xclyAz@h|mFC7UVw$oIv*F(*}F)`OUht_SNV zdc5Bbd&19-4k@Qz(s*SAD%S;B%foRhUskob(|FClWI3pxh6u!l&o1@b035rw?Qxv~ zr>dmg{$=JdP=?g@8@N6h6eQMPNZMls8dj`12=`}$C*?5?Kk8<2TwV)B#wmYdpa1;& zi0DNZaDD;XMQ7d=u%j_EOPtfuWBE#xcs->glG*?Eg0sh+YGmT|^U4v6gjquv@}=aG zeP-JZCSF#m%jI7Ul4Qsi$-g^Z88e}{`zLQ0b^q0#;E=DJs!>VHLk7!p77rNwQ+$~U z95O@y&=T_+Fx=1RYO}N@Az!IJP}FyYfuK+mMeo)R?I9yoD6~*+{<7$LiXM|@i2+`a zeBa@C9Fw5RKhU?|3=cO--jlDjs?JN6p}rqA$7@nuDNT>Q3`g{L5~seev)Q;KBt@Dk z|0-YnPWg@xRbCDZ=3mq7 zs!x9>(Xlyu2f{dt;T+ zSI-VGaf{p0UAd;5jh}MU+l#)kMIZk*y0>E**wrMTEuF6p5-+*jU%l{kkH?Ru@MYB{ zUxzZr1HXe@!`T%(81LKnuU~{N4Ijp!@4v7-akqUkgWf205aM`m1Y`V={JW2zmlCpf z8Zvz#@gDoixhHIPX}|qKkmz`2OH7Lj{JQOiMZ>d7&=l(uHsrYv*ogIGkS}~O=ac(5 zi6)AElu7FDDGyYk-XpJk+AQMveR+02W=h$?-oaZY&Fy7g!OI8R=F-6Zv$H<|Z@qlqhj-Y^Z!x`t zSDUnx@Cx1y>V?0|?JR5&OoscRiv!e$w9en;T1rUED%X6CAOpCh+~Z8@&j0xYU2G(6 zFWA~7`eC8Rg95}kAF?hfU%|w|!W0p}=tR)(yi^Uh6Fm74@+pt4E+(7Xqf|&W@fBr# zCoKWo%jn{!EwUJW2DZ9@eJr?dD^0wFAmnk2R)m3x!*p{*AU#8^KXMbL7DGCx=L;bp zP?lnd)0%iu#1XBM7&2KST@=a0PjBm34o50bR%p-N)|NuN(dLOEelL|-NG{IdhQI`x)wrz@-9g>L?C~sYduKeIv zT@&N%5hawbI1+bwkOLBhUp>RV>5Qb}99$khc1Fzb+p24J&OsdU20avBHW#VDug1+< zITuO5X`5`N=OL~5)z0;SZioTStZr0|JCcbLRBbkpLr?ybpWX@QbTw|>%y?ACky6Tc zn5-IGtXvXDd%@(mM2j63^!#$Z`uNvT=j36vWU7H)M=Uj8!KbsLL*PwZZOjL4BGK3L z%whmG?X}r>9*8X&e~j(C2keY1?r0zB3^dDkzm`3m4Y>LUXWYC&uciTWK3elx z?)^2xjSPA__aOP(@=}cXCCR`0{PMwBm)46(K5$K_d(GnWO`Ma8qg$7S#{j+c4cnc> zRpAc5>LJ&HpMsJDi83BG0fc^H<_lQ^zef+7)k4whx7aK^^FkG>EDsE;d{#@*+nn0w zy>_Y>9Jt@ESS*OfYePtWmSBZETz295(Eur0A3(L4lhfSAU{t@|LNfabfFEo7{ZT1P zKtGK&!>dZp15}$uXNx!uE^O3CI04MuVTL?iIV(4OImkyrX5rNe;;(}V9GM4{?FGwd1bZg00v&|qQQCD z0?hvR+1t20><*JYF{0~)SNGR!2U;PA_^t?q1A=Z8yoY~20y~_L z1#$D!u8ROyH?PrU)Y`Og%6et)_l*E63Mc!*`m~{$`Gt=Frp}&h_k*Yr47G$k1mT@S z-ksO;f;F<9hHvN5z}#m_;o)TiQ2fJPB`aqrY?q;MJ7 z@@W~Sse(aUkHpWi+za*&+LYq_W;S3MH+tT4y$(jszSTH2GmcZcC-qR?bpv2C>4$;K zp%Cy%V@R3(D_btT5Xc=ddHVOq(hU2Nm14+&;IYi{b2n~uz*8~icsM>Y_w1zCNsRIE zu>bnmH%>F@olJ)KsF*Ux&l(AS#oe*Y@k8?OK7LkRMSEzx51bHOrJ;f%3!_o%|zNeg6y zoGK{#T*>cuKWliwguFECFh?3M@m&4$CRgR*v1?B^#`Dv7t;Ec)BSfGq{zASFPTexK zIvW!NFg5XK`2}!AV3-6bN&~lJe#@y$pADQoj=QWf65SKKG0_C$WV=+UpZ6fZPTw72m#kk zPWab8SA|2)PClpq>@k?HekeZUmnS$l=HibvB|h+-=2@Ye+nXr*x}XTpCGl!-h{4s% z)+4GZ`VQ1ew(a1Qx9KwSH=mwrSLmDz|MX_WZeU`{-IW&Neo*d7;{w*& z7K+~UVPv$fpDHZcBtN&&?+Hc!HbA3o#|bak?^xAExhxv5(FV#l=IJWHgI7=42#u%l z!u9{>!=S>Zs4$tqw?U5Z9zLzx(Ewj}63Ye)1WQKFs@?$PSqn@A)=UHQWYmi~_iP4A z(>>-i*c)>uXR_|z$ZTX=oiNkE)^pymE5E=?yAycvx!uowpTW{qOPb)(-Lnxcx=Y z54QRpRcW?qq39n4SUmG|QH9CNtDBz6S5fq`hQ6l$nO^X9hOtcXZW^x*ofi+NX(&Mb z^)*VOnlxU7o<9^86iO-E3)}_APdr%=Z4&}4Y=4Z1s0`=&*Z6@U`1-60Al7LY0Y8`a zi(Oss!W|@jl{)l3xzI`m^Hv9GqX8*s9 zHT!lWh$&waUfc=s@nXtX#Tm7yrzbJxEBxD%0-uXayh#4t@v_wVE&o&37v6k%aE_@% z6Nl^{YlqOmL5`}h>_Xi1Wh?8b{=xk+1qnF%JcMfEv*x6zbDK%Ho)b&f82l4SLBr`dP2p1(NFq26UJ+ z*4rhIXy=JJFy(8lqMVzZJX5|NMu@nsb7SIVe(!F>xNA)LsxKEzSyaoUA9#644VcZ8 zuRRM--J2`G#EazL9WUG*q~r_Hx9B4=$0m+S`YV5jNh<*6jq%z-1q$UMV@7joz%1Fp zLT#=q2>EL8g+h_a!ACzgQ}ok2GFQ|}sKUW{3QJ24*HiSky-Uy=62DvZfRG!1csP_nkHwq7Ix5bwj$kg&tm1KfTlEWx$Mo#fi*=rUoalYs=?dG7tb@O^*m7w|q| zZ{zquf^RXve^*{o`akr+KjIMcLK~g%6cGIr2?#rYEmtl_bL;;V5OgEB<+t4LM^f$U zYX55__Lch6@o{J)$0cS;v)l z_p7sa)4Dn%Zg2BSY1)3tG~!Rd1+j%e^wXE&%9;~a+B>*5O(I&!voNvFfsY8^dHMgh z(xnZFhV*5)PAGO(_YSTpQ+N#5fys4G|Acha8~1&qIJ1>=-gI3N)j#<9(pO(7++KXeaYr*nkGqW? zSAvIYy}Zpv->2w>rtERtKg1gjHyN=`diWy_*>9T*EFYz>Re;}O@-E-JyP4{5Yo6~7r=yOW{M~31-zEPCMD`_>lpc|(ER|g6vvTo9NdEg4G1T_8j>BCi{yP&fDIZM8H z5x{ZuI{*F=a#o;%vW^yjxL)So%&N!PHD;6L-FN{ zfTD+$;={`-&~>{v`|-7M@F8G=wsy>c#-Yxj#tvfz=55be((}{NcTO&A=^dRU$29C$FmFd2Zd7g z$LUS|A?lPMlABdX(G&XmRgeNkEdD$tiK55(eWmqnLhl_0iMWh1uD-pRt8cpz`njGg zBNmKLh9r7fUxfpzINtmmC@Z0`Orl+6Um_X``iK8yZD&t+k$LS1gO`uX6?v*IAJ0*K>qF&n;z0 z#B*FB{DTVp;Azw=XVjE(fR7ig04|>NcqIG1`>}V^y6YBPe!S9g8)IaNbRRuXn#j3M z_honT`afBAk5wTWLU|%kC*3>j&y`*G=l60J9&W{BxJ*y)I?>ZvI0;K{;QBhg@57Zo zv!875q)Q&pEP}R_SJDO3YCru6xMC;teYl)J!r|V*^5_x?oO2@~#$bvD?jip8|A1F`$VoHjHPkgt`Re*ZJLIdLGYx zamweb5z}(8#l)rK+Y+Yd@yK=O15*z&J&$eVo!`5#hB7oEl7Dwi2=^3U;tNyrE)@BO zH*v`RN#2P18USr*S$ywu`nh_le{jzbem^+Q*3OM@Pzyyrb@QrWJ_ohoZdP)_#vLyx zdR)UZ@`gmjaXDw^1LZX~6T<+hiMTu1sidz<2g_qK*ef?r0+`Wp zAvH^aAFVl)no)OpQd8d6Ln=+~XBsHz#pQ2`Q z>@cc7(YWv814;pK)G{pR!!?SYsNW;tw$nwGDrHw{;c3(JG_bISJu7XJJ4mA|3=W!nZ3240tk%?|uLFoj`h zh;Bz}cGwSM&~GgD`aOB@S_b_?*&8Ji!r~0PNdDdN!aaGyeWB7^*O>!_8aZVD936P@ zGe1;<8o8=JORm;Z{iB41M2Yx8PD|UHg9n-@`cAiUvzW7ovz6g}?XHr5+r z@ojT;XuJjtJ>oBQNC7t5C{)MKr9L9HJVpstB~*wSKpHRHFdF0sUrbKB8abzlqL17;ZroKZ zWq53zjdgxUJ=LGM+}V?!AqXTqC#1NRnp(;egQ%AS*Ykpm#Q4g!%#fZJZ@Hc!WVxPe zNs^~5L%4-%gCF-ft%N7Wet+)w?*JaTv5%?`GGU#uTZ1hf$$^G}pti(i)O#|!$a=!M zLRm<@PiNp3CpvQoPM$_*jvz8gC{XnHcy08f%DAOTI)Ex8fAi_x?7C(Smp!j^*EPPP zL;l2c?VCN(kiI-!dp#~itaorF&Eqj#ccROBY--@e?E0u<--m1Lq^VB5gKN@!9>ewQ zh47c2N*6D5*Vt)9Lnw1Ab;506A2<|VJiWMgaJ^r^W4N-j`swr(E?!929A}~-efnjX z7{pC`)7|gtp9W;?a_2Ez$cZ}}{sdgP9(^CKU%qSRU^(5cAMu7~PI_;wFn+?}c0Uq^ zc_WC8fzC1g+z-!ArKU!BlB-k0A-)(K`T%%uT(Z`DQ7mO@AXMNaFk`xuN9!~Ta6d}g z|Hzf)pcc`qE3_HJiND}p936R${nhW~+T0EaFtl>%q(`Zp?0}gUWur=S*_vvL_P_KT z%~&LV_pD){d<1imyyE`Ny0-rO4Ao{>Sm4D6Z<(r1Bo@`ZiAk^PIkVYdPb@ro!1|Iq1i?K70)GGII((fN`n{2f&pHq7Om7&hS%rc&$?~?5`3@~ z-#7hi1lW0b`{p4|ez4B2L~7}&R*HVtSL>fwb}2))dHMd?CJhulo=oz9H|+Rq7O>fn z#_OcYC9NMqW1!EVrOhjBX}s{EFH*N;1v61 zd}x{}_xfruLimZO`A2b1^s2!vt@7pEAe2#Ho6x?bNna2!Wr^R7u!Jk@M?>B#Bq%UV zLshCuFVpy#rXf+UDY@T-`!n!zc5-{y0hmb&M|^5DnV80;KczN%-1n7n47^SZEi}Av zPK1FM$-g^Z_@eP+Us&N$}4N7BpUj6M+Pr#5ZI0V4G49Gi9fQ-2PZb(2!n0ffYM?8n|%u<7I`8WH5W9 z((#RFrgFA^Rm&Vy*v95hZTzCAH5t4#JzD#~iUkA?s(9N?eaarXVMO$ZmrN%VNx|DI z&a#;Dh3Ax+#FQ8DZPT6!p(e%drk;EA^Y-QyG3m`^iZo6aGv&+7ENx~^DpS5l{@wAy zQYkFHg$6a$9s0GztcA0%?EUVQM`MAXvwVY!q$*4idOOl-dIh*T=|;bS552(ZefN)+ z>i9uHtRw4MDSFr8;(L-es6ymKXXB3TPbhl)eVFMDSATRbpEHHVi(g4#_|-vUU~p5> zkUR$(FDwVn_al-}5s&681q~Mqh6oM}s; zJADT<#<_A(Ps;g+#WzA!yg02T`_F1yw$#SMD7 zN|`>$Ae7XVV^7UX9t|zGwhvz_@fwU;ocnCou92{$W&55f!skKx_KCMY%#x0m7xEacjQn{?e*&)E{|n%Ht2B0H@1(1LIFI2v zJuZ{T&C>%1N?uroZu$}p>C3-e+tzGQ?;TucBX|r~@R)}=J(VtAz*QRA_u;y@r!^j5 zwEvBUyPRQ{2RYPGmZ6QiO!MpKj)BHLO1b`}yQmo^wh2k<)R8b0UyZ*H5<8pveLS=Z z_~H2qPlAtl8pJK&X0pzrvGO*c0m`Hex7?QVu*T@%X2_D?_;1wVva!RgmpIZr|-lzC>Z0+{Q7^gj!t04MRH0YD6&n=`pHcMqTK8xl_g9#=Kfc;UF@j)3(AqG@^uAAn=213BfTOS$db2NZs`Mb)RpfuIS66K3q60n>N6YzR^1q});R&abWIZdV`$-|f=l`!Q? zwd`kFJi6UTcc-D~GkfZ0)iLEOWKH7DY-=W7B>(Ps;b~Gu{0aFoY38VgT~iw|aV~dT z9Syi_ec3kqu*cwwy7=y2C9{a|q)%M=dfGzKp9{M&XW3CrC|PpnwZPdbiXL~)YV?L4 zSoa&!co9xpRX`?gu`9>zF!esP#q&XqfLIoq(61HLb<)XetKEZ7N*fPQ*SG9>K494j zKdQ|AxkJeGKV+Fa;!YWD(l8j`7Al{cPG~L%4~vng1%yRv4cAFcnHz4y#EY=ton*)t z5st&e>tE3mVLwc~$n+iE@xr&$c;*N3w9VDmI7B+6xt#Av#`v)CMsQ9}P2sxwCaV6A zFfWWke@JYbQcBSWu}=<`eAxt~c7M<>v)e_{-N&CVdzR8!!_iPyZ@hnYuw_#57)@rPA@{bl+c7ffBpAQ zm&vcj@))k%eflQ}xc*ul?EF8QkEkxD@53eVI`wt$;0lfBFE=D~%Hq+dH_Ht>H0Tnx$fMdkPmX40$`(5)J9gGs7EZP5W^Bo0t!bwp*}0I)GMO zg{>#22L+6Q>kE%cmA;@24cP89?eel^xX zk}44X_|aaUdyk^WJ+g$*eT*F_C%hp z55V`-k~O{q;{g!`lLKC@!hB#-fu?(K(%!Fj;Ar{BI48H$KzQtjUuVREIZ?^Cd2VjBb^( z6SNs`7Z6|R0I*$bMaLF{w{0m!ezFOi&`0~x9lDLcnuIxzuDpx{PL1ZTw}9DTQ~k_I z3*00?PhIjZO&jlgpPy;yNjs&|QkB3oysF)Q9X2A3Nq;?UzRqM9rr~wb$}7$%2FWrE zuO$EOcwyO>54HKh&?}R4o%XbF#(XlDYvzjsBeRoM?sr#%uP>ZzF19TL@?OI4N5#7X zJesV|9|}2D-2I~8M$zN@Be|$UEWLDFMg>KW2g7dlfxi8cFD{FE&f)$XA;-fHG@Ys$ z18pB%{$PEH&KpFy^vJ+T3$aYL}lA#bcUp%LU0XTweG2`OTOz>Kr346u$p$ zEK^>F*Iqf4RK=99>`jhVr4mfMyx%Rq=8BG-?k;?9Eq9eaaAL|A$-g^Zs8pPPc*!3c zTgC}!XSQ-q4B1@%W=K51f-HSA)Zn9w)1L0Ft^l3Wrmd}&_aWp9HH=liklufX-A>VG zo!weED@7eH&^uKX%dVp6af9hV+HgkX?HSF@`}=}xI9mt^!dMfgJ{JmFp&;UbnAz|H&f5#aj-U8m|p44+sjNRA`YwTmp{`#(95 zhCBqpwjCCmJ%pmib{F4sV=z_5b|GK0o!++?H~fdthDAJjQJgj`lF#{9$KFlXG`4X0 z^2(66-P-8RpO~&$tS1`Mm#1sLC;iCm9bDBZJcet0^l|tn;Of}W_u=A8OLy!YT*|MqGY&} zI8l|FtM%o-n-%|mx>+%?S#{F>Y=6WMwM&R0w<{|`k=f{}II?(7T_|!4znZ^%$wDNm zKSofMyd?~IieL3BE)GKsP}{B;axVPnD#QRMDA*FU8o7)Uh#D5HLmKd_+s7AdKrHa9 zfpMj&h$K3?#E^OMGj<~l`03~Y8+Ic~IBoE8k3EPtzIF4N8**ug5>6my#M*~Q3ShTI zzE<6jB;W)aL(0;T1iTIpIN}G93jB0|g~>6b6(=y(bj?A$@vGn0_~#)?_|??YyYi4! z{ArG6ob!kzPOwz;;6)@5uVjmi;bkNWfAE!I^KT%jIKhsz!W)Pqep~T^l)H#GPH^7Y z@*X096Ij1jDMfDLS4X6eeuy;SSA}g?lp~V(ZH+EVDv%BMBWAVjuS5j!Q^9v-&k#ra zwq1=gUm&@7C6m}&-XW6ssbWEUDFXMOEdJ@935mS^ zPx6e(EhBhP_Iv7t%e)U9)}Jrh*SpvCn?2lSnDMAN7h+FJ^%O32z~Gym|JQS zAA^cpQ|oGqJ%9tI5wT8#x_S{bfNl%m^be`IZ$AU9(#JOEN&#*WqWcEUu;aWokrQ#8 z(;pb`KI1oWk|k)gN_q2*tqCS?JpM9u(HnM8eVXoGI1=~z>~Myu(Z{Q~;rivh3{}}R z^YVwF*-{Mp-JhoKiZ9#Bp#M>u1-`q>F-##O|L&@_-obMdI`aG>wlR3Tbt~tE=|a6} zqdbB4kMBt$)5gL$_7)XO=L+y~fslgzdOLt8GRX;m{u3vS9M#Z9(G$_ST5ysFJHYka z6N>&njh9KW&C%r(H-S!zFIO{#(3gC8{9?zM9qg*tMJ=n9*ubLd)~OrWEO32W%D9}9 zPkW44qp8@N=jWJs{k|sSJ2;<-m+g^Z-k~l`yuS615t;vpi5I?gyeAVcl7Dx+PKTpj zl_*|$V~c|}Tefm8-_u{eKz|tj=|yCEb+jKD&Zxx*FHH|3}Rvrs<00-yJV>U#hJ~Hv7Xplf4aGOIs*Y zV6>#&3OfOPNNlO_l4>u$Wk7fU)PMY|q|L95qCfoBLhJWMefV@%+M!J)4HW%<@ghdh z{}#OPt>g~+Kyuup+m}g?bCsgXglW%>!k5@WoHp&@8R^k2_{jTl3*tw*jBkm?CZkIc z1pl3jrA>RdY1582?cqOzEp&a*v{erj=F-pWULQ1nI_hpo`#2KO|Ke1A@cO-?%P2S0 z>5(*rA_$(nVG#X%yd7`5mC45aUyi+-u5HicvgFltO*Try^iNFJ&h8@`(wC=eKk_?| z^$sqx!#svd-Q9XuPvIh!p6;`2^7=oSUC&M?8bawn)CsqNec(`5WS7-DxLUG#4A-$c zadUbK7pe4ghl{-aPvG)QCmPb1r3;RKZ_zur4jtt&TpPZoXZ#7cP=gTu>C15OU9KcN zzF6JeU;n8`uHA7S!=IXu*w;6lOh2T5s>fb3xOd2Q zrs^817W75J=_|sh;j6M2-KOy(RLjw@Tm$pF zsRSFE?%j27U(Vg_E)&FJd30{A1^j1CwtRkL3ec=pVxxnDz=FZY9zFAx<|IG;6klC- ziwaNd2z1_9-1Lq8bDRAGm4|uk^E>p-KNzwYUKi(2l*M1U&h)yNskk%1=LB;%H1kDY zGy`t#Vwif;+l$1I$xKyOZTZ|}b1kOo`ptBoqN^S=oRs9>9WOkLc3c2F?{-A1-n)%6 z%~n=r-nZr8!#9D;?jq`N7F66QbGsT`5s)~6#b|(!}s_5 zst!9>4B_u*QcKa}yK8NtRS!G`t_7_M5jm9Ua5X&oSAG}}$py}}(y?aS7 z5d^_gvhSCt_~W_BZRzLZs$M`K^*x^0i{BT1pC~Cg+?XEEcn)%?2s2A3OsFrT$1~xL z$D+bFbL-6(81Ihbdas=3?h3dSV413=TW!I6r&7x=GW&ssqK$?8N*j)8M!nEu_ix=o)Wv7 znY+tt)VTUiEn(8TN>tyB{=(c{elM00ZQ4kNsgva2eLTCgUaVfP5dha-%m}zVu#F=) z%gQSJa2&CB!9sOt-tJMmi&G6o+)8WR?ioRZZ?*Fk;`DGcP+r62c5UjHY$>)ZbYriRTk(?xp+*Mf^YhAVvM z9L+xg*Y!h0Lnwm~b;506A2?hpd+XdgxI{1W7%n&0u|&jV4_=0(6!k{B>@)g4Tpwc! z%z6h`{uLg>HR!k;`%g$$$?3ii*Zzt8dcA`ysDQ_C`RW(H`4e!>KhyW&+OYH|v6;fZ zZOD_k%44|ZSb=jrg^M(YcYhg@*Z;{><1-hozC1I0wrN{K@1(1^kjHS9Ur{;TQ@D8X z+qLg(-9>dk@AEez=xJW7LjdZ2nAR5w_rR(6sbyE~~{NB8x*#~BQ zTYYex9&M(=_d?fIgcbAmZ_6m6%~XU(z$kcz^Ksu~%_ktPBx6R!lz46s)j5DAMawVT z0h*#l{(3an0et#e>9ujmX)yiz*P%K}9-QpCr}YL3ih)s!tOP&&o(IfrVx@wrAlR6H z|KYBHAMA#Q8&791W}2xK>}A^mW0__u4XsVZGg_Et)att@JvN0ggA=PxXXU6SF+GNo zaubBkSu=wZk4)DQG;U{l43Yf1o2hV*Ao&33cltLG%0JR=(h-aiWqOwfTK2?ultbujG}kN_GbHl z25UP%Twf#*6Uikb$Fd7o;li%r|vc*|((t&z5C$4)|8Q z3D9WU0VIDO*rihG0OrZ47j^Er0B-NFEqpEz#>pI=@+B$tBYQd?yyyjv-JCuxbGZc& zYe-J{Xf4p=VA7y>4;^e!GgkLtZA1MPrdR9sGtA_lPmOLrF<>Nve(C*R@f9GIL4SDO zb62f07Q?Ja^6!qWCt!IJSZ2b!$JqWKs8nc7Eh&zzM(nep&o(m)pcowF3(8lVy!S>qQecxDwh$TtUIx#AM;Yy;D_5xL z+mx97pIUNf)zd}kjQ;;x7`S8nGba70g7c>oTbX$M)HYds&Rw2?7sDJePNYL@WkdVmUB9{HA8Mts!{$R@B2+csUp&^J`sfr#k2S zvKK;uMUT15P%2=>2Bbqvn9EQ_3HhKqT+lz`-Rp7k9(Begi4_YuPp&Iam z-Ix}=)@KwwZk9&3cH~YQ5W-AOFWYoCOXI=2FKM$h?vcBUHcJ!t(`J9dOik0{%VfoY z)bnr;z6U{68K-xkm#Oe9hKjVAy}L)=Ni2nbrY<~`29*eJg(jIx( zqC-gp{d~OscEJ=s;~!Cm%!-JH^yTSVjp)76y@RW=gvW5n zE%3hgCuY~JxB5O@?T=!%^$xDY2Rw!g#u(536L4kT>HBbvzh{27cW@b$@ffbt`_Fy; z6L6)L_kFmo^Cw;I9bBo;cnsIbadt!|gn!$%c&3KORa%+Y(x1!nUmRQ#4Ikf(@P&&m z?G)itqUXD7W2a`@9a4movb5hI^Xd66;Wi)#qiwLw?wO7Kt0Dtq}+FG{$-l4SUdI z{Ukq2P$2UHQ=eR9H(tt2`~RY!TyY4x7s%94GE$RHoJnQsC#rezSxUp1`U%Ou`+OJQ z-?2UbhN|1S2R(1&1idmIOF!-$MdO9fJJS?lndbRfU!7>Y z@a)=FEST+Ftb5~3Eq9;4rChh7bZ`yNrW6mH9<`cQCIEouuD{gevI+c9=5Jdw%9sOk zHI+U(H?oN&QidRF()j&FzZ=+Nou7$_s#LP|s(=lD8#6i4D}jXH%OaWjiKg;+xysi} z`-|SBo9pyWr7`qL>tpg^$q7t1NEMs#_lxA1ZjGmn1WQ+&k7nqTB>(PsVcVNNLiLl< zf+p4!0W`v{WcsY54okuQ$(rT6qt#*iGF8p#6Ym1fgz{wDqaH*AWlW@HJrQko|+mjB8zNI67#>(-U@cGOq`kBXSD0*Cv zkMf17H_u+$piJZCZ=m@3Ook#-n8>+cTUa$2w9v4v^P< z2Svn>k73H!tz)bkdli}T)pX`Ew0zEV15Ohd>$A+9Xtyj{Pc{g{dRUx^m&h^1_rOwgyp~Ti#PA8rRe{~&)R=4|HB^m z6Z1>6$3*4&va;hnj|jan=Fe&5`W`b8{({@W8>w>s|d2GZu-3a~PEeKX)8E+aoDSG3nzxzuCYCx}93XeXxJ)!8|B%$ev zpZUV*9ZT-DNzr(X4n_mDek;PqV%czhTNbsRjQere$wNmR=qErqeckEj?-ID-L}$Ut zCks$h&UTQpIA)^rHajqD=DA~=n^<(1p6KtaO^uDu*!NO z^2)WK$NJ>pTIU9nD5k!VB~x)&u|LyZvk6PNq|CI}ygZP5{plR0y@s_+^5 z?~WJWjMal}76ihrH&VjeJ3Ba6HRFw3J|_T?)5$Yk)imHw<@u++eX0N(kBKPme&j{y zC&58*?0n5ZTeg0q==CqHda}zz6S_6K8QW%5QS|&6U~hb($&G@W;bU0dTpBrM_AP-ajtV zyM92IQ-sS?9T%^0;9h<0hdDaK0I|DYNssYT=(J)fJ!ba*fY}pR6F)HV>d5BIci+O) z|6QY6)H^mY_5aiTmYt7z%CxVM{JZ1zEEiQMBLd;6IlI1EPig1qYb`3A(*XxXB>zTt5PMc-MiXWdU<6Uv}@l4NaaDSCWY z(@J0HRy($OtT>I=eCykSQ@$(0#u2ScRhww}!lIG}7|X*_JWF#u$ds(J?Fd@Veeck_ zWR7<^H2UDfDM0u6j&CvcVc@W8)5?Q`WvKViUZdEUY}L$R+G{4hTmCYzj!9pdHo<(?QD&Imh7c252vhoTXgh%10eec&!BWV2%he5=TRBmLD8pvd?IvXk^%f$ zpf07(c|+0v?GN043I4|(=;^eM7p|;Y=ZGp%!7|jz|5(-kr(ZfQv0L4{X`Sh3ZiRWJ znU9;hK z_un3LF~zXd7Y?(~%J37W=ey7L=04SGQG{(D)R5negDAJjTv`@_3&9KII6T5miYh(7*pbZ^Hr5IlA1 zF|RyzV3;Zqwzv5VTT#dv9{9yv?~;n`e;--GT<@}ZS^Fb#By;|Hpm9vUUj@whZh%C} zz4tXt`eo9_k#4^G3UD^|L*hMZEXV<>`V`Y@~^J^PDb^^taW3X2R+yZ ztc=lwGCmscnZwrYu`?fol1(l<&z^K4^#6N7@bYq}U20>0Q1tIpa|VT_X~GXnX3x*o zt)}RG&!Qg6w@`m_G4G3cd^BE*Fh|o*itw89#8-D-(0I+k@}k^lL2S#mmzmE2Yw9gg z^{`+nUW3>LdMc)Aj{BCR8Ed8b&@#5u-Cluo}u}*l}W$q>e_^eh0L@% zK{@7AOvflP){{y8-SNWL|C0mZ^;JHV`HpRz{u9GhgUW*emZfOIIt|#7b4W$z=3^lJ zwk&H+iz~7IUls)C$s%FSU)m{pEFbUSCz?>`N21n@@zoUl-wNP=3I5vy|LuYQ_P`(N zf!@u+Rey4;%_|-EP~ar7KQXN{Ya<%cmv!756Mqfw9bA>acnp_o#niz)g^L$DuH#*z zA(Y{WI^j034;(BP8EEwmu0({}SGubbyqHLg@dXFwo(A76EDeJg#4O|4hf$k8l1%+GI zD`@IMub=vZE*t3@=$mq19(>3E9c^u01B2$@1i}A5f4kz}4?*OP>q}t$qo=Dk-Ybai z541c%Jfut6@LFh>vQbzVUDikMXGL@=Pbgj1rEIujMVGSC*x0|8^|hbHbt&s?Sk&kL7yK(9D?p^5e zz5dcmUCIXPSGttPSLJspPgu)capUqcoKQ%7kN-8du>6bsqn9eN9mCMiL4NL$k2wD3 ze(cW=z1(%spv1w06Ma$khDN%&2IbuF3vPLQ?f{%0`eMhwpRi4cDjR9*>Q0Cd?ou|i z8`Pz2q$B#*vYxi~h%RNFA<|vSx}Rmbl=U{qbt&unsnca0U0ow7&A*oQw2Q`eDeFXy z>r&R8Zq%i$_kAK=*45S1K03KeS;vLlrL3!M(WR_cXW6B!zsb5w* 1: + out.form = 'Madland' + out.fragments = Polynomial(energy_release['EFR']) + out.prompt_neutrons = Polynomial(energy_release['ENP']) + out.delayed_neutrons = Polynomial(energy_release['END']) + out.prompt_photons = Polynomial(energy_release['EGP']) + out.delayed_photons = Polynomial(energy_release['EGD']) + out.betas = Polynomial(energy_release['EB']) + out.neutrinos = Polynomial(energy_release['ENU']) + else: + out.form = 'Sher-Beck' + + # EFR and ENP are energy independent. Polynomial is used because it + # has a __call__ attribute that handles Iterable inputs. The + # energy-dependence of END is unspecified in ENDF-102 so assume it + # is independent. + out.fragments = Polynomial((energy_release['EFR'][0])) + out.prompt_photons = Polynomial((energy_release['EGP'][0])) + out.delayed_neutrons = Polynomial((energy_release['END'][0])) + + # EDP, EB, and ENU are linear. + out.delayed_photons = Polynomial((energy_release['EGD'][0], -0.075)) + out.betas = Polynomial((energy_release['EB'][0], -0.075)) + out.neutrinos = Polynomial((energy_release['ENU'][0], -0.105)) + + # Prompt neutrons require nu-data. It is not clear from ENDF-102 + # whether prompt or total nu value should be used, but the delayed + # neutron fraction is so small that the difference is negligible. + # MT=18 (n, fission) might not be available so try MT=19 (n, f) as + # well. + if 18 in incident_neutron.reactions: + nu_prompt = [p for p in incident_neutron[18].products + if p.particle == 'neutron' + and p.emission_mode == 'prompt'] + elif 19 in incident_neutron.reactions: + nu_prompt = [p for p in incident_neutron[19].products + if p.particle == 'neutron' + and p.emission_mode == 'prompt'] + else: + raise ValueError('IncidentNeutron data has no fission ' + 'reaction.') + if len(nu_prompt) == 0: + raise ValueError('Nu data is needed to compute fission energy ' + 'release with the Sher-Beck format.') + if len(nu_prompt) > 1: + raise ValueError('Ambiguous prompt value.') + if not isinstance(nu_prompt[0].yield_, Tabulated1D): + raise TypeError('Sher-Beck fission energy release currently ' + 'only supports Tabulated1D nu data.') + ENP = deepcopy(nu_prompt[0].yield_) + ENP.y = (energy_release['ENP'] + 1.307 * ENP.x + - 8.07 * (ENP.y - ENP.y[0])) + out.prompt_neutrons = ENP + + return out + @classmethod def from_endf(cls, filename, incident_neutron): """Generate fission energy release data from an ENDF file. @@ -327,72 +410,22 @@ class FissionEnergyRelease(object): # Check to make sure this ENDF file matches the expected isomer. ident = identify_nuclide(filename) if ident['Z'] != incident_neutron.atomic_number: - pass + raise ValueError('The atomic number of the ENDF evaluation does ' + 'not match the given IncidentNeutron.') if ident['A'] != incident_neutron.mass_number: - pass + raise ValueError('The atomic mass of the ENDF evaluation does ' + 'not match the given IncidentNeutron.') if ident['LISO'] != incident_neutron.metastable: - pass - if not ident['LIF']: - pass + raise ValueError('The metastable state of the ENDF evaluation does ' + 'not match the given IncidentNeutron.') + if not ident['LFI']: + raise ValueError('The ENDF evaluation is not fissionable.') # Read the 458 data from the ENDF file. value, uncertainty = _extract_458_data(filename) - # Declare the coefficient names. If we only find one value for each of - # these components, then we need to use the Sher-Beck formula for energy - # dependence. Otherwise, it is a polynomial. - labels = ('EFR', 'ENP', 'END', 'EGP', 'EGD', 'EB', 'ENU', 'ER', 'ET') - - # How many coefficients are given for each coefficient? If we only find - # one value for each, then we need to use the Sher-Beck formula for - # energy dependence. Otherwise, it is a polynomial. - n_coeffs = len(value['EFR']) - - out = cls() - if n_coeffs > 1: - out.form = 'Madland' - out.fragments = Polynomial(value['EFR']) - out.prompt_neutrons = Polynomial(value['ENP']) - out.delayed_neutrons = Polynomial(value['END']) - out.prompt_photons = Polynomial(value['EGP']) - out.delayed_photons = Polynomial(value['EGD']) - out.betas = Polynomial(value['EB']) - out.neutrinos = Polynomial(value['ENU']) - else: - out.form = 'Sher-Beck' - - # EFR and ENP are energy independent. Polynomial is used because it - # has a __call__ attribute that handles Iterable inputs. The - # energy-dependence of END is unspecified in ENDF-102 so assume it - # is independent. - out.fragments = Polynomial((value['EFR'][0])) - out.prompt_photons = Polynomial((value['EGP'][0])) - out.delayed_neutrons = Polynomial((value['END'][0])) - - # EDP, EB, and ENU are linear. - out.delayed_photons = Polynomial((value['EGD'][0], -0.075)) - out.betas = Polynomial((value['EB'][0], -0.075)) - out.neutrinos = Polynomial((value['ENU'][0], -0.105)) - - # Prompt neutrons require nu-data. It is not clear from ENDF-102 - # whether prompt or total nu values should be used, but the delayed - # neutron fraction is so small that the difference is negligible. - nu_prompt = [p for p in incident_neutron[18].products - if p.particle == 'neutron' - and p.emission_mode == 'prompt'] - if len(nu_prompt) == 0: - raise ValueError('Nu data is needed to compute fission energy ' - 'release with the Sher-Beck format.') - if len(nu_prompt) > 1: - raise ValueError('Ambiguous prompt nu value.') - if not isinstance(nu_prompt[0].yield_, Tabulated1D): - raise TypeError('Sher-Beck fission energy release currently ' - 'only supports Tabulated1D nu data.') - ENP = deepcopy(nu_prompt[0].yield_) - ENP.y = value['ENP'] + 1.307 * ENP.x - 8.07 * (ENP.y - ENP.y[0]) - out.prompt_neutrons = ENP - - return out + # Build the object. + return cls._from_dictionary(value, incident_neutron) @classmethod def from_hdf5(cls, group): @@ -419,9 +452,9 @@ class FissionEnergyRelease(object): obj.betas = Polynomial(group['betas'].value) obj.neutrinos = Polynomial(group['neutrinos'].value) - if group.attrs['format'] == 'Madland': + if group.attrs['format'].decode() == 'Madland': obj.prompt_neutrons = Polynomial(group['prompt_neutrons'].value) - elif group.attrs['format'] == 'Sher-Beck': + elif group.attrs['format'].decode() == 'Sher-Beck': obj.prompt_neutrons = Tabulated1D.from_hdf5( group['prompt_neutrons']) else: @@ -429,6 +462,44 @@ class FissionEnergyRelease(object): return obj + @classmethod + def from_compact_hdf5(cls, fname, incident_neutron): + """Generate fission energy release data from a small HDF5 library. + + Parameters + ---------- + fname : str + Path to an HDF5 file containing fission energy release data. This + file should have been generated form the + openmc.data.write_compact_458_library function. + + incident_neutron : openmc.data.IncidentNeutron + Corresponding incident neutron dataset + + Returns + ------- + openmc.data.FissionEnergyRelease or None + Fission energy release data for the given nuclide if it is present + in the data file + + """ + + fin = h5py.File(fname, 'r') + + components = [s.decode() for s in fin.attrs['component order']] + + nuclide_name = str(incident_neutron.atomic_number) + nuclide_name += str(incident_neutron.mass_number) + if incident_neutron.metastable != 0: + nuclide_name += '_m' + str(incident_neutron.metastable) + + if nuclide_name not in fin: return None + + data = {c : fin[nuclide_name + '/data'][i, 0, :] + for i, c in enumerate(components)} + + return cls._from_dictionary(data, incident_neutron) + def to_hdf5(self, group): """Write energy release data to an HDF5 group diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 84d9bcb2c..d2c9290b6 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -52,6 +52,9 @@ class IncidentNeutron(object): Atomic weight ratio of the target nuclide. energy : numpy.ndarray The energy values (MeV) at which reaction cross-sections are tabulated. + 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 metastable : int diff --git a/scripts/openmc-ace-to-hdf5 b/scripts/openmc-ace-to-hdf5 index 743cf4a50..90031f0fa 100755 --- a/scripts/openmc-ace-to-hdf5 +++ b/scripts/openmc-ace-to-hdf5 @@ -25,6 +25,13 @@ 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 optional --fission_energy_release argument will accept an HDF5 file +containing a library of fission energy release (ENDF MF=1 MT=458) data. A +library built from ENDF/B-VII.1 data is released with OpenMC and can be found at +openmc/data/fission_Q_data_endb71.h5. This data is necessary for +'fission-q-prompt' and 'fission-q-recoverable' tallies, but is not needed +otherwise. + """ class CustomFormatter(argparse.ArgumentDefaultsHelpFormatter, @@ -47,6 +54,8 @@ parser.add_argument('--xsdir', help='MCNP xsdir file that lists ' 'ACE libraries') parser.add_argument('--xsdata', help='Serpent xsdata file that lists ' 'ACE libraries') +parser.add_argument('--fission_energy_release', help='HDF5 file containing ' + 'fission energy release data') args = parser.parse_args() if not os.path.isdir(args.destination): @@ -111,6 +120,14 @@ for filename in ace_libraries: # Continuous-energy neutron data neutron = openmc.data.IncidentNeutron.from_ace( table, args.metastable) + + # Fission energy release data, if available + if args.fission_energy_release is not None: + fer = openmc.data.FissionEnergyRelease.from_compact_hdf5( + args.fission_energy_release, neutron) + if fer is not None: + neutron.fission_energy = fer + print(neutron.name) # Determine filename From 2fff5cd3d049152d30578befee9d259e70330376 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Tue, 2 Aug 2016 18:18:06 -0500 Subject: [PATCH 06/33] Improve fission energy release documentation --- docs/source/io_formats/fission_energy.rst | 53 +++++++++++++++ docs/source/io_formats/index.rst | 1 + docs/source/io_formats/nuclear_data.rst | 21 ++++++ docs/source/pythonapi/index.rst | 18 ++--- openmc/data/fission_energy.py | 81 ++++++++++++++++++++++- 5 files changed, 164 insertions(+), 10 deletions(-) create mode 100644 docs/source/io_formats/fission_energy.rst diff --git a/docs/source/io_formats/fission_energy.rst b/docs/source/io_formats/fission_energy.rst new file mode 100644 index 000000000..d73a95f60 --- /dev/null +++ b/docs/source/io_formats/fission_energy.rst @@ -0,0 +1,53 @@ +.. _usersguide_fission_energy: + +================================== +Fission Energy Release File Format +================================== + +This file is a compact HDF5 representation of the ENDF MT=1, MF=458 data (see +ENDF-102_ for details). It gives the information needed to compute the energy +carried away from fission reactions by each reaction product (e.g. fragment +nuclei, neutrons) which depends on the incident neutron energy. OpenMC is +distributed with one of these files under +openmc/data/fission_Q_data_endfb71.h5. More files of this format can be +created from ENDF files with the +``openmc.data.write_compact_458_library`` function. They can be read with the +``openmc.data.FissionEnergyRelease.from_compact_hdf5`` class method. + +:Attributes: - **comment** (*char[]*) -- An optional text comment + - **component order** (*char[][]*) -- An array of strings + specifying the order each reaction product occurs in the data + arrays. The components use the 2-3 letter abbreviations + specified in ENDF-102 e.g. EFR for fission fragments and ENP for + prompt neutrons. + +**//** + Nuclides are named by concatenating their Z and their A numbers. For + example, U235 is named 92235. Metastable nuclides are appended with an + '_m' and their metastable number. For example, the first excited isomer + of Am-242 is named 95242_m1. + +:Datasets: - **data** (*double[][][]*) -- The energy release coefficients. The + first axis indexes the component type. The second axis specifies + values or uncertainties. The third axis indexes the polynomial + order. If the data uses the Sher-Beck format, then the last axis + will have a length of one and ENDF-102 should be consulted for + energy dependence. Otherwise, the data uses the Madland format + which is a polynomial of incident energy. + + For example, if 'EFR' is given first in the **component order** + attribute and the data uses the Madland format, then the energy + released in the form of fission fragments at an incident energy + :math:`E` is given by + + .. math:: + \text{data}[0, 0, 0] + \text{data}[0, 0, 1] \cdot E + + \text{data}[0, 0, 2] \cdot E^2 + \ldots + + And its uncertainty is + + .. math:: + \text{data}[0, 1, 0] + \text{data}[0, 1, 1] \cdot E + + \text{data}[0, 1, 2] \cdot E^2 + \ldots + +.. _ENDF-102: http://www.nndc.bnl.gov/endfdocs/ENDF-102-2012.pdf diff --git a/docs/source/io_formats/index.rst b/docs/source/io_formats/index.rst index acab7e893..39b38fcf2 100644 --- a/docs/source/io_formats/index.rst +++ b/docs/source/io_formats/index.rst @@ -15,6 +15,7 @@ Data Files nuclear_data mgxs_library data_wmp + fission_energy ------------ Output Files diff --git a/docs/source/io_formats/nuclear_data.rst b/docs/source/io_formats/nuclear_data.rst index ba6a54eb1..7544ca1f5 100644 --- a/docs/source/io_formats/nuclear_data.rst +++ b/docs/source/io_formats/nuclear_data.rst @@ -55,6 +55,27 @@ Incident Neutron Data from fission. It is formatted as a reaction product, described in :ref:`product`. +**//fission_energy_release/** + +:Attributes: - **format** (*char[]*) -- The energy-dependence format. Either + 'Madland' or 'Sher-Beck' + +:Datasets: - **fragments** (*double[]*) -- Polynomial coefficients for energy + released in the form of fragments + - **prompt_neutrons** (*double[]* or :ref:`tabulated <1d_tabulated>`) + -- Energy released in the form of prompt neutrons. Polynomial if + the format is Madland or a table if Sher-Beck. + - **delayed_neutrons** (*double[]*) -- Polynomial coefficients for + energy released in the form of delayed neutrons + - **prompt_photons** (*double[]*) -- Polynomial coefficients for + energy released in the form of prompt photons + - **delayed_photons** (*double[]*) -- Polynomial coefficients for + energy released in the form of delayed photons + - **betas** (*double[]*) -- Polynomial coefficients for + energy released in the form of betas + - **neutrinos** (*double[]*) -- Polynomial coefficients for + energy released in the form of neutrinos + ------------------------------- Thermal Neutron Scattering Data ------------------------------- diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst index 36f161b3c..5d45c799d 100644 --- a/docs/source/pythonapi/index.rst +++ b/docs/source/pythonapi/index.rst @@ -335,6 +335,7 @@ Core Classes openmc.data.Tabulated1D openmc.data.ThermalScattering openmc.data.CoherentElastic + openmc.data.FissionEnergyRelease Angle-Energy Distributions -------------------------- @@ -368,21 +369,22 @@ Classes +++++++ .. autosummary:: - :toctree: generated - :nosignatures: - :template: myclass.rst + :toctree: generated + :nosignatures: + :template: myclass.rst - openmc.data.ace.Library - openmc.data.ace.Table + openmc.data.ace.Library + openmc.data.ace.Table Functions +++++++++ .. autosummary:: - :toctree: generated - :nosignatures: + :toctree: generated + :nosignatures: - openmc.data.ace.ascii_to_binary + openmc.data.ace.ascii_to_binary + openmc.data.write_compact_458_library .. _Jupyter: https://jupyter.org/ .. _NumPy: http://www.numpy.org/ diff --git a/openmc/data/fission_energy.py b/openmc/data/fission_energy.py index c22fef9a9..d7fe5fa72 100644 --- a/openmc/data/fission_energy.py +++ b/openmc/data/fission_energy.py @@ -1,7 +1,6 @@ from collections import Callable from copy import deepcopy import sys -#from warnings import warn import h5py import numpy as np @@ -199,6 +198,82 @@ def write_compact_458_library(endf_files, output_name=None, comment=None, class FissionEnergyRelease(object): + """Energy relased by fission reactions. + + Energy is carried away from fission reactions by many different particles. + The attributes of this class specify how much energy is released in the form + of fission fragments, neutrons, photons, etc. Each component is also (in + general) a function of the incident neutron energy. + + Following a fission reaction, most of the energy release is carried by the + daughter nuclei fragments. These fragments accelerate apart from the + Coulomb force on the time scale of ~10^-20 s [1]. Those fragments emit + prompt neutrons between ~10^-18 and ~10^-13 s after scission (although some + prompt neutrons may come directly from the scission point) [1]. Prompt + photons follow with a time scale of ~10^-14 to ~10^-7 s [1]. The fission + products then emit delayed neutrons with half lives between 0.1 and 100 s. + The remaining fission energy comes from beta decays of the fission products + which release beta particles, photons, and neutrinos (that escape the + reactor and do not produce usable heat). + + Use the class methods to instantiate this class from an HDF5 or ENDF + dataset. The :meth:`FissionEnergyRelease.from_hdf5` method builds this + class from the usual OpenMC HDF5 data files. + :meth:`FissionEnergyRelease.from_endf` uses ENDF-formatted data. + :meth:`FissionEnergyRelease.from_compact_hdf5` uses a different HDF5 format + that is meant to be compact and store the exact same data as the ENDF + format. Files with this format can be generated with the + :func:`openmc.data.write_compact_458_library` function. + + References + ---------- + [1] D. G. Madland, "Total prompt energy release in the neutron-induced + fission of ^235U, ^238U, and ^239Pu", Nuclear Physics A 772:113--137 (2006). + + + Attributes + ---------- + fragments : Callable + Function that accepts incident neutron energy value(s) and returns the + kinetic energy of the fission daughter nuclides (after prompt neutron + emission). + prompt_neutrons : Callable + Function of energy that returns the kinetic energy of prompt fission + neutrons. + delayed_neutrons : Callable + Function of energy that returns the kinetic energy of delayed neutrons + emitted from fission products. + prompt_photons : Callable + Function of energy that returns the kinetic energy of prompt fission + photons. + delayed_photons : Callable + Function of energy that returns the kinetic energy of delayed photons. + betas : Callable + Function of energy that returns the kinetic energy of delayed beta + particles. + neutrinos : Callable + Function of energy that returns the kinetic energy of neutrinos. + recoverable : Callable + Function of energy that returns the kinetic energy of all products that + can be absorbed in the reactor (all of the energy except for the + neutrinos). + total : Callable + Function of energy that returns the kinetic energy of all products. + q_prompt : Callable + Function of energy that returns the prompt fission Q-value (fragments + + prompt neutrons + prompt photons - incident neutron energy). + q_recoverable : Callable + Function of energy that returns the recoverable fission Q-value + (total release - neutrinos - incident neutron energy). This value is + sometimes referred to as the pseudo-Q-value. + q_total : Callable + Function of energy that returns the total fission Q-value (total release + - incident neutron energy). + form : str + Format used to compute the energy-dependence of the data. Either + 'Sher-Beck' or 'Madland'. + + """ def __init__(self): self._fragments = None self._prompt_neutrons = None @@ -453,8 +528,10 @@ class FissionEnergyRelease(object): obj.neutrinos = Polynomial(group['neutrinos'].value) if group.attrs['format'].decode() == 'Madland': + obj.form = 'Madland' obj.prompt_neutrons = Polynomial(group['prompt_neutrons'].value) elif group.attrs['format'].decode() == 'Sher-Beck': + obj.form = 'Sher-Beck' obj.prompt_neutrons = Tabulated1D.from_hdf5( group['prompt_neutrons']) else: @@ -471,7 +548,7 @@ class FissionEnergyRelease(object): fname : str Path to an HDF5 file containing fission energy release data. This file should have been generated form the - openmc.data.write_compact_458_library function. + :func:`openmc.data.write_compact_458_library` function. incident_neutron : openmc.data.IncidentNeutron Corresponding incident neutron dataset From ef89d40479015f673851020b2749993f771bbe33 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 3 Aug 2016 11:35:42 -0500 Subject: [PATCH 07/33] Add fission energy release test, score docs --- docs/source/usersguide/input.rst | 21 +++++++++++++++++++++ tests/test_tallies/inputs_true.dat | 2 +- tests/test_tallies/results_true.dat | 2 +- tests/test_tallies/test_tallies.py | 5 +++-- 4 files changed, 26 insertions(+), 4 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index e530097de..9ed30afeb 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1809,6 +1809,27 @@ The ```` element accepts the following sub-elements: | |:math:`\gamma`-rays are assumed to deposit their | | |energy locally. Units are MeV per source particle. | +----------------------+---------------------------------------------------+ + |fission-q-prompt |The prompt fission energy production rate. This | + | |energy comes in the form of fission fragment | + | |nuclei, prompt neutrons, and prompt | + | |:math:`\gamma`-rays. This value depends on the | + | |incident energy and it requires that the nuclear | + | |data library contains the optional fission energy | + | |release data. Energy is assumed to be deposited | + | |locally. Units are MeV per source particle. | + +----------------------+---------------------------------------------------+ + |fission-q-recoverable |The recoverable fission energy production rate. | + | |This energy comes in the form of fission fragment | + | |nuclei, prompt and delayed neutrons, prompt and | + | |delayed :math:`\gamma`-rays, and delayed | + | |:math:`\beta`-rays. This tally differs from the | + | |kappa-fission tally in that it is dependent on | + | |incident neutron energy and it requires that the | + | |nuclear data library contains the optional fission | + | |energy release data. Energy is assumed to be | + | |deposited locally. Units are MeV per source | + | |paticle. | + +----------------------+---------------------------------------------------+ .. note:: The ``analog`` estimator is actually identical to the ``collision`` diff --git a/tests/test_tallies/inputs_true.dat b/tests/test_tallies/inputs_true.dat index 9d67bc0d0..f5390eea1 100644 --- a/tests/test_tallies/inputs_true.dat +++ b/tests/test_tallies/inputs_true.dat @@ -1 +1 @@ -930af242a043f2676a000dbc5a2db6b148edcb31ed8c87dbaa35a8efb37a3be8cff30cdf4dc03f9c5c7eb4021f7e4c3327e64681cdd8fd8722c95c69db850227 \ No newline at end of file +1bef757d276362fdcd9405096b4cdcbd894f9215ed406493486a45193729be446c9a12242c887f89b6e209ec5beaaacb04dee2fd61e72b4f5c6a8712b776ed6e \ No newline at end of file diff --git a/tests/test_tallies/results_true.dat b/tests/test_tallies/results_true.dat index 7aa65e1c1..7d4763f98 100644 --- a/tests/test_tallies/results_true.dat +++ b/tests/test_tallies/results_true.dat @@ -1 +1 @@ -a51db2a4efc681805f85968e04411dc33beee0532c202f5179b9a82880ab60a75e53fa9141c81045ea1d2842372f2d8da900326f09382ea61dd80a3c9b43bba1 \ No newline at end of file +f8ad60a94994e6b126b64b8a3e10ac3647314fff37558f9e32f174e6baf944aa2e617ad15f3aa6d6b3fc4e4ead944afb3ca2ea606b40e50bf96e7ad86c86a12b \ No newline at end of file diff --git a/tests/test_tallies/test_tallies.py b/tests/test_tallies/test_tallies.py index ba0098513..387be8af5 100644 --- a/tests/test_tallies/test_tallies.py +++ b/tests/test_tallies/test_tallies.py @@ -123,8 +123,9 @@ class TalliesTestHarness(PyAPITestHarness): t.filters = [cell_filter] t.scores = ['absorption', 'delayed-nu-fission', 'events', 'fission', 'inverse-velocity', 'kappa-fission', '(n,2n)', '(n,n1)', - '(n,gamma)', 'nu-fission', 'scatter', 'elastic', 'total', - 'prompt-nu-fission'] + '(n,gamma)', 'nu-fission', 'scatter', 'elastic', + 'total', 'prompt-nu-fission', 'fission-q-prompt', + 'fission-q-recoverable'] score_tallies[0].estimator = 'tracklength' score_tallies[1].estimator = 'analog' score_tallies[2].estimator = 'collision' From 4ff005b439481439f74b1c5d8e2ae9dda50d07d2 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 3 Aug 2016 15:53:18 -0500 Subject: [PATCH 08/33] Fix incident neutron energy subtraction --- openmc/data/fission_energy.py | 4 ++-- src/endf.F90 | 4 ++++ tests/test_tallies/results_true.dat | 2 +- 3 files changed, 7 insertions(+), 3 deletions(-) diff --git a/openmc/data/fission_energy.py b/openmc/data/fission_energy.py index d7fe5fa72..334231a9f 100644 --- a/openmc/data/fission_energy.py +++ b/openmc/data/fission_energy.py @@ -603,12 +603,12 @@ class FissionEnergyRelease(object): data=self.prompt_neutrons.coef) q_prompt = (self.fragments + self.prompt_neutrons + - self.prompt_photons + Polynomial((-1.0, 0.0))) + self.prompt_photons + Polynomial((0.0, -1.0))) group.create_dataset('q_prompt', data=q_prompt.coef) q_recoverable = (self.fragments + self.prompt_neutrons + self.delayed_neutrons + self.prompt_photons + self.delayed_photons + self.betas + - Polynomial((-1.0, 0.0))) + Polynomial((0.0, -1.0))) group.create_dataset('q_recoverable', data=q_recoverable.coef) elif self.form == 'Sher-Beck': group.attrs['format'] = np.string_('Sher-Beck') diff --git a/src/endf.F90 b/src/endf.F90 index a836a5439..833082e1c 100644 --- a/src/endf.F90 +++ b/src/endf.F90 @@ -60,6 +60,10 @@ contains string = "events" case (SCORE_INVERSE_VELOCITY) string = "inverse-velocity" + case (SCORE_FISS_Q_PROMPT) + string = "fission-q-prompt" + case (SCORE_FISS_Q_RECOV) + string = "fission-q-recoverable" ! Normal ENDF-based reactions case (TOTAL_XS) diff --git a/tests/test_tallies/results_true.dat b/tests/test_tallies/results_true.dat index 7d4763f98..818d99a92 100644 --- a/tests/test_tallies/results_true.dat +++ b/tests/test_tallies/results_true.dat @@ -1 +1 @@ -f8ad60a94994e6b126b64b8a3e10ac3647314fff37558f9e32f174e6baf944aa2e617ad15f3aa6d6b3fc4e4ead944afb3ca2ea606b40e50bf96e7ad86c86a12b \ No newline at end of file +a6e5480c66e6510687bf281983b3f387da45005bb08d8cd0aa629e53c5a42a1b2fa8d76345ad2df497d85f2b273fbc5edc36105a71a73c1107479ca0af8329f6 \ No newline at end of file From 14b8ef6a1f6bd248a22e6c4cf806d27fcd812fca Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 5 Aug 2016 13:48:12 -0500 Subject: [PATCH 09/33] Make PyAPI Polynomial class, remove Constant1D --- openmc/data/function.py | 91 ++++++++++++++++++++++++++++++++++++++++- openmc/data/product.py | 28 +++---------- openmc/data/reaction.py | 6 +-- src/endf_header.F90 | 30 -------------- src/nuclide_header.F90 | 8 ++-- src/product_header.F90 | 8 ++-- src/tally.F90 | 44 +++++--------------- 7 files changed, 115 insertions(+), 100 deletions(-) diff --git a/openmc/data/function.py b/openmc/data/function.py index bea6f5e9a..e827e1283 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -1,3 +1,4 @@ +from abc import ABCMeta, abstractmethod from collections import Iterable, Callable from numbers import Real, Integral @@ -9,7 +10,53 @@ INTERPOLATION_SCHEME = {1: 'histogram', 2: 'linear-linear', 3: 'linear-log', 4: 'log-linear', 5: 'log-log'} -class Tabulated1D(object): +class Function1D(object): + """A function of one independent variable with HDF5 support.""" + + __meta__class = ABCMeta + + def __init__(self): pass + + @abstractmethod + def __call__(self): pass + + @abstractmethod + def to_hdf5(self, group, name='xy'): + """Write function to an HDF5 group + + Parameters + ---------- + group : h5py.Group + HDF5 group to write to + name : str + Name of the dataset to create + + """ + pass + + @classmethod + def from_hdf5(cls, dataset): + """Generate function from an HDF5 dataset + + Parameters + ---------- + dataset : h5py.Dataset + Dataset to read from + + Returns + ------- + openmc.data.Function1D + Function read from dataset + + """ + for subclass in cls.__subclasses__(): + if dataset.attrs['type'].decode() == subclass.__name__: + return subclass.from_hdf5(dataset) + raise ValueError("Unrecognized Function1D class: '" + + dataset.attrs['type'].decode() + "'") + + +class Tabulated1D(Function1D): """A one-dimensional tabulated function. This class mirrors the TAB1 type from the ENDF-6 format. A tabulated @@ -239,7 +286,7 @@ class Tabulated1D(object): """ dataset = group.create_dataset(name, data=np.vstack( [self.x, self.y])) - dataset.attrs['type'] = np.string_('tab1') + dataset.attrs['type'] = np.string_(type(self).__name__) dataset.attrs['breakpoints'] = self.breakpoints dataset.attrs['interpolation'] = self.interpolation @@ -258,6 +305,10 @@ class Tabulated1D(object): Function read from dataset """ + if dataset.attrs['type'].decode() != cls.__name__: + raise ValueError("Expected an HDF5 attribute 'type' equal to '" + + cls.__name__ + "'") + x = dataset.value[0, :] y = dataset.value[1, :] breakpoints = dataset.attrs['breakpoints'] @@ -304,6 +355,42 @@ class Tabulated1D(object): return Tabulated1D(x, y, breakpoints, interpolation) +class Polynomial(np.polynomial.Polynomial, Function1D): + def to_hdf5(self, group, name='xy'): + """Write polynomial function to an HDF5 group + + Parameters + ---------- + group : h5py.Group + HDF5 group to write to + name : str + Name of the dataset to create + + """ + dataset = group.create_dataset(name, data=self.coef) + dataset.attrs['type'] = np.string_(type(self).__name__) + + @classmethod + def from_hdf5(cls, dataset): + """Generate function from an HDF5 dataset + + Parameters + ---------- + dataset : h5py.Dataset + Dataset to read from + + Returns + ------- + openmc.data.Function1D + Function read from dataset + + """ + if dataset.attrs['type'].decode() != cls.__name__: + raise ValueError("Expected an HDF5 attribute 'type' equal to '" + + cls.__name__ + "'") + return cls(dataset.value) + + class Sum(object): """Sum of multiple functions. diff --git a/openmc/data/product.py b/openmc/data/product.py index dd276daac..116905f7a 100644 --- a/openmc/data/product.py +++ b/openmc/data/product.py @@ -3,10 +3,9 @@ from numbers import Real import sys import numpy as np -from numpy.polynomial.polynomial import Polynomial import openmc.checkvalue as cv -from .function import Tabulated1D +from .function import Tabulated1D, Polynomial, Function1D from .angle_energy import AngleEnergy if sys.version_info[0] >= 3: @@ -36,7 +35,7 @@ class Product(object): yield represents particles from prompt and delayed sources. particle : str What particle the reaction product is. - yield_ : float or openmc.data.Tabulated1D or numpy.polynomial.Polynomial + yield_ : float or openmc.data.Tabulated1D or openmc.data.Polynomial Yield of secondary particle in the reaction. """ @@ -47,7 +46,7 @@ class Product(object): self.emission_mode = 'prompt' self.distribution = [] self.applicability = [] - self.yield_ = 1 + self.yield_ = Polynomial((1,)) # 0-order polynomial i.e. a constant def __repr__(self): if isinstance(self.yield_, Real): @@ -120,7 +119,7 @@ class Product(object): @yield_.setter def yield_(self, yield_): cv.check_type('product yield', yield_, - (Real, Tabulated1D, Polynomial)) + (Tabulated1D, Polynomial)) self._yield = yield_ def to_hdf5(self, group): @@ -138,16 +137,7 @@ class Product(object): group.attrs['decay_rate'] = self.decay_rate # Write yield - if isinstance(self.yield_, Tabulated1D): - self.yield_.to_hdf5(group, 'yield') - dset = group['yield'] - dset.attrs['type'] = np.string_('tabulated') - elif isinstance(self.yield_, Polynomial): - dset = group.create_dataset('yield', data=self.yield_.coef) - dset.attrs['type'] = np.string_('polynomial') - else: - dset = group.create_dataset('yield', data=float(self.yield_)) - dset.attrs['type'] = np.string_('constant') + self.yield_.to_hdf5(group, 'yield') # Write applicability/distribution group.attrs['n_distribution'] = len(self.distribution) @@ -180,13 +170,7 @@ class Product(object): p.decay_rate = group.attrs['decay_rate'] # Read yield - yield_type = group['yield'].attrs['type'].decode() - if yield_type == 'constant': - p.yield_ = group['yield'].value - elif yield_type == 'polynomial': - p.yield_ = Polynomial(group['yield'].value) - elif yield_type == 'tabulated': - p.yield_ = Tabulated1D.from_hdf5(group['yield']) + p.yield_ = Function1D.from_hdf5(group['yield']) # Read applicability/distribution n_distribution = group.attrs['n_distribution'] diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index ad707d276..d35599579 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -5,13 +5,12 @@ from numbers import Real from warnings import warn 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 .function import Tabulated1D +from .function import Tabulated1D, Polynomial from .data import REACTION_NAME from .product import Product from .uncorrelated import UncorrelatedAngleEnergy @@ -465,7 +464,8 @@ class Reaction(object): idx = ace.jxs[11] + abs(ty) - 101 yield_ = Tabulated1D.from_ace(ace, idx) else: - yield_ = abs(ty) + # 0-order polynomial i.e. a constant + yield_ = Polynomial((abs(ty),)) neutron = Product('neutron') neutron.yield_ = yield_ diff --git a/src/endf_header.F90 b/src/endf_header.F90 index c4b4ef3ad..8d8aefaa3 100644 --- a/src/endf_header.F90 +++ b/src/endf_header.F90 @@ -30,17 +30,6 @@ module endf_header end subroutine function1d_from_hdf5_ end interface -!=============================================================================== -! CONSTANT1D represents a constant one-dimensional function -!=============================================================================== - - type, extends(Function1D) :: Constant1D - real(8) :: y - contains - procedure :: from_hdf5 => constant1d_from_hdf5 - procedure :: evaluate => constant1d_evaluate - end type Constant1D - !=============================================================================== ! POLYNOMIAL represents a one-dimensional function expressed as a polynomial !=============================================================================== @@ -72,25 +61,6 @@ module endf_header contains -!=============================================================================== -! Constant1D implementation -!=============================================================================== - - subroutine constant1d_from_hdf5(this, dset_id) - class(Constant1D), intent(inout) :: this - integer(HID_T), intent(in) :: dset_id - - call read_dataset(this % y, dset_id) - end subroutine constant1d_from_hdf5 - - pure function constant1d_evaluate(this, x) result(y) - class(Constant1D), intent(in) :: this - real(8), intent(in) :: x - real(8) :: y - - y = this % y - end function constant1d_evaluate - !=============================================================================== ! Polynomial implementation !=============================================================================== diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 8ff83482e..514958a07 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -10,7 +10,7 @@ module nuclide_header use constants use dict_header, only: DictIntInt use endf, only: reaction_name, is_fission, is_disappearance - use endf_header, only: Function1D, Constant1D, Polynomial, Tabulated1D + use endf_header, only: Function1D, Polynomial, Tabulated1D use error, only: fatal_error, warning use hdf5_interface, only: read_attribute, open_group, close_group, & open_dataset, read_dataset, close_dataset, get_shape @@ -283,11 +283,9 @@ module nuclide_header total_nu = open_dataset(nu_group, 'yield') call read_attribute(temp, total_nu, 'type') select case (temp) - case ('constant') - allocate(Constant1D :: this % total_nu) - case ('tabulated') + case ('Tabulated1D') allocate(Tabulated1D :: this % total_nu) - case ('polynomial') + case ('Polynomial') allocate(Polynomial :: this % total_nu) end select call this % total_nu % from_hdf5(total_nu) diff --git a/src/product_header.F90 b/src/product_header.F90 index f20adf0d4..a69929473 100644 --- a/src/product_header.F90 +++ b/src/product_header.F90 @@ -5,7 +5,7 @@ module product_header use angleenergy_header, only: AngleEnergyContainer use constants, only: ZERO, MAX_WORD_LEN, EMISSION_PROMPT, EMISSION_DELAYED, & EMISSION_TOTAL, NEUTRON, PHOTON - use endf_header, only: Tabulated1D, Function1D, Constant1D, Polynomial + use endf_header, only: Tabulated1D, Function1D, Polynomial use hdf5_interface, only: read_attribute, open_group, close_group, & open_dataset, close_dataset, read_dataset use random_lcg, only: prn @@ -109,11 +109,9 @@ contains yield = open_dataset(group_id, 'yield') call read_attribute(temp, yield, 'type') select case (temp) - case ('constant') - allocate(Constant1D :: this % yield) - case ('tabulated') + case ('Tabulated1D') allocate(Tabulated1D :: this % yield) - case ('polynomial') + case ('Polynomial') allocate(Polynomial :: this % yield) end select call this % yield % from_hdf5(yield) diff --git a/src/tally.F90 b/src/tally.F90 index ec46a6194..0a910a969 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1,7 +1,6 @@ module tally use constants - use endf_header, only: Constant1D use error, only: fatal_error use geometry_header use global @@ -247,24 +246,17 @@ contains ! reaction with neutrons in the exit channel if (p % event_MT == ELASTIC .or. p % event_MT == N_LEVEL .or. & (p % event_MT >= N_N1 .and. p % event_MT <= N_NC)) then - ! Don't waste time on very common reactions we know have multiplicities - ! of one. + ! Don't waste time on very common reactions we know have + ! multiplicities of one. score = p % last_wgt * flux else - m = nuclides(p%event_nuclide)%reaction_index% & + m = nuclides(p % event_nuclide) % reaction_index % & get_key(p % event_MT) ! Get yield and apply to score - associate (rxn => nuclides(p%event_nuclide)%reactions(m)) - select type (yield => rxn % products(1) % yield) - type is (Constant1D) - ! Grab the yield from the reaction - score = p % last_wgt * yield % y * flux - class default - ! the yield was already incorporated in to p % wgt per the - ! scattering routine - score = p % wgt * flux - end select + associate (rxn => nuclides(p % event_nuclide) % reactions(m)) + score = p % last_wgt * flux & + * rxn % products(1) % yield % evaluate(p % last_E) end associate end if @@ -289,16 +281,9 @@ contains get_key(p % event_MT) ! Get yield and apply to score - associate (rxn => nuclides(p%event_nuclide)%reactions(m)) - select type (yield => rxn % products(1) % yield) - type is (Constant1D) - ! Grab the yield from the reaction - score = p % last_wgt * yield % y * flux - class default - ! the yield was already incorporated in to p % wgt per the - ! scattering routine - score = p % wgt * flux - end select + associate (rxn => nuclides(p % event_nuclide) % reactions(m)) + score = p % last_wgt * flux & + * rxn % products(1) % yield % evaluate(p % last_E) end associate end if @@ -324,15 +309,8 @@ contains ! Get yield and apply to score associate (rxn => nuclides(p%event_nuclide)%reactions(m)) - select type (yield => rxn % products(1) % yield) - type is (Constant1D) - ! Grab the yield from the reaction - score = p % last_wgt * yield % y * flux - class default - ! the yield was already incorporated in to p % wgt per the - ! scattering routine - score = p % wgt * flux - end select + score = p % last_wgt * flux & + * rxn % products(1) % yield % evaluate(p % last_E) end associate end if From 264cb33b3635e54f8d9958b829c0d0660b25258c Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 7 Aug 2016 14:31:18 -0500 Subject: [PATCH 10/33] Address #690 comments --- data/fission_Q_data_endfb71.h5 | Bin 67543 -> 67543 bytes docs/source/io_formats/fission_energy.rst | 7 ++--- openmc/data/endf_utils.py | 5 ++-- openmc/data/fission_energy.py | 35 +++++++++------------- src/tally.F90 | 12 ++++---- 5 files changed, 25 insertions(+), 34 deletions(-) diff --git a/data/fission_Q_data_endfb71.h5 b/data/fission_Q_data_endfb71.h5 index 89e51b70783034b40ce4e207fdb3e34b5d5d33cc..7cc5a86b5249c1abacfbf23bf777639825164c9b 100644 GIT binary patch delta 8209 zcmeI1dst0bAII0JeJa;N(Hi2I(GcO>x-v}VQYl*Gami5bq>D=PzH9&1XaCl3{r1{x?QW$C zw^D@_uC1yUZ5{R}IY0khD%jWnFgY#I8NM_Z%z+j85dat$6%jEZau#SI1K?`f%85-B z3ZE&n!e&j0iu6<{^n&c7IZW`4JU7U<&tNkJ1)7|?Fb^lQwQ}9`rh^1w_nob}Npg2< zlk3!07v?lZ*Te$DlNqs{_>?)8l;ZbChve`24H?oGcsKxd$m9Ujnvj=jPK(U>l6g5klF4@3E& z8eMeLAR$m~7sLfok0x7Vq0P047l=_iz0$_Q>Wc`FVb(?3BjH-rrh{A*^=T5saQ=%i zi0hLs@GB%-F8==6f5h}_jp4I%r~D0=a*}i=Crvt?a~mY|ICR`im(-@M?3!+4n<$~P zk9xDNg@d!~sqTA+?lM1Jv_q@Jkinz8y@Lb#0eR`At1@tkWITBaYzv!#D3VLo;7a#e zTi{Fbs|Dz>>#QAcG@}MbU+h(Yz9hA5fc$`4H(lcDpq$`vO2HvwgNFItaM1TDzWEgFy`|zOyPC0Xi{4 zggUW$1PEv4a&@y1kj^aT$V)=NHipRA^E?!EV8wy)KaB%@nZcsBITOK9cITTZ8zzA} z#J~nDI8r_pq%c@ntAKEj&J6C{S`i6qSkddC>r9Z(iVc766$2gv5KP)VykPSq2x09i zcW(RzObO_jyqxa(aL!{ORo0y82T_La~Zs@neWBV zy)B&)oOaClM@ayi?xMP2=qo=wozEERjNrB|yNz>KJtr8WZC4EYq<=e(&k(vH_)^i> zk?%)gHw+!}1JCmr`tC^BIx5+g#I%!MsU**SeP8J8m+ih8KK{$9oX3fSO-Ce|@%Q#0}&kG@$5}$Kj@s_XCN8ew{h z2LvC@o$*N$z=E-JG=?W?dIs_t>0=Ol9#yhZnn8mxIuygAgkKFP)5gDSO^_=WTgJ`EY7ernTO!Hf)p3N zA2;V=Sh?`6@x@Xk%p(K zy~$N-$zE{UBn-P8Z5SlYs7i{1Z_Q$C3_jFvPsgNeKJj5tmmj3t(|(*=fMBcI3&uHA zim`Gbh8cMYQ~8Y8WCYD}?;6Kdis7{gL+Aug7#OceEKE)oi;jyCTxOYJ$e#8lSE;2w zZ83SNB<*>>L7%aN3~Ig03~+M=f;C=444bg}02a&*X&5fkRw(&SB(Fqp&W8uy%tDee zP>10^>Noj(hISQ#Z{t>V=QBR7!mx#p_gFroEFC3!TF?H*JJ@d^kT0t6UI;X70h(gm z@j1C?kX0W!I`~$`CI5d_{y(%TPc&?`9D*C90y$gppHS>g1X{um;H4UL!WxZaqk5waMp)Z2~u#?WJpuIUt!8cb`@50NJdV zS3$FIF7` z>salgnls0M3!`nZTW|u%<OD~w35t0FF-bH*Z8{r1(?A4;orJl1&C+Ze ztJem;1n#W1abTA!aF-R!y^5>BAFP-cH@^ntGgyzHyS?6l8?1cV;nxqqkrmeue46=R zZMDcZk@~X|oq+S!l8Nt1Lnr>TzA_WR&_?%~^XH3@by!$)zvFPd`FhUW=ao%fiq$cO zRfOcc)o=@jb4CVx^8I+V1;OLmQw7p~(7qhlis2Ex-FiM_Lk@zuwy}m~s-z#2f5))p zvsZdPV_+_V3(Hh(rTw5U_*UC6jO=7LfY11}4Z#_$zb)Z2PUm64=JDV``pE5^nC9gp zc)n@QG(IC{FNWQM1CH_i2-!y#)X`6~q!~1(9rly8?3+sCm#!-EVl+8*rTgRH$eq79 z8`_&(p;l}yTyy}#L+!e_)BLr|q?87$49!RU1y!{av|l<^rIP9XSaR?J3e22Ie46IhsX#HYJH{Uk}Q z`r)F3Qgji)vzlu|_>6Uh7~Z>7`WxSm`IpGf^?u|AK4ZXT40pU&+VdHzD+oI6SiXYK z_;3Zos92d{Kc|rD%Y&;3&U?3&lN>Z2zhA>JZ1!>ECP#{qT7=-0LjOv>r(a)3!u|Y1 zo5)v%w)*tzQgJ1BYH%OJL)C(D_xZ)0h=pS@hPBZ$LmEo+$>2ku9sjK zcPypQ;O%Sj(3+gOv7ySy~1#OXNL^FAAwZ}F8MmgR+>TctmSJA14??D z^BEssldltZRLl5`Q`K0o%o$*Bh}vHgQDe{VM+6W1DDwE6b)PUSy?@^LlByqh%9#H* zf`KiTua#!dm<;%be5APT=W6H={U?>lNjE*|Dj9;Cv)UdtFifseOZtLH{~*JOSyvY~ z@8{dc8zi+MiyeB9#ZFaYah-wBM77p+?gNo_@A6QjJ^oc@s$G=3DnwcpS}D?M>nf4f6udUk3bl(i z{bOB4-4HkFGzFe}&M4+>ijY%6my7tx8!S(rV9Sk=8t4B+}Z9 z6(a50e`THaiXi+qI+0doC#H)k)WMmbwH0lk9uuK5r)zR@5FcrkX7?76)`sPL)~tz^ zG_lG3Tt%}@t30wzq*al0WljzmNMr_Thkc(bnATj}FVfohgCgzH`H)Dv-Z?DNZcC1c zw7c6ek@k3U+@QhgEV;}-@h_X_5wv^b)M^*HsiKEbyY96RX}5H%IvqS&t#+H$P^8_5 zG!|(ObyHKVQ7Y|4T4nteFeCxFApvUCFI$TYG%4*wT06OOoem09t36sdn`({nbXSp9 z#dH^Gb%2LRYid2~l290`)+lp(e!fnnYD6znsaE=$N==EMsnjkXU@Bd{3pACk2jiNb&>2QZ7R2R$ueRSvP+Gz zp+_yJTw)g^gz91CwJ5p8!y@ap?890Sn@gVG?~I=M>3975{PT?0(5PUYlI z?rLfIQR1w)S&6ABgDj2al4>*fqjJsT5P#pW4vsvp%~v0%b?vamT)%UAh;)-OI!iyr z+`D_5-}AouIM-Ocl?8#5N46hlQ)YMMDWlz<*L^u^>{vgfwL|nHF-PdR8TTmZf>0wD zBa$3PNW)oC@%)f6ybp!GCsytKjs$w=g?O^7jp)o*YZmpln_kjUp_jW@w2NOdY@)uV z^G2x2TKhX2r0GJlP}sc7g(JyDNf{Nf`jRdL8Y;efpG~sqMoItLyd#oH;&ZyJ8-eRv zGF~!ClE?b#R-w}2n?GB?1o0syb|+Gw-#kunPL4*o%%W4f6mZ?FPw9oFzK;%cWTN<( z5`e%5L(gnr#-tDGjir*|Pomz68P$ityEmTAVPo9;Vrh2&L(Te(zV?zTePQ2J>C)xp z@9VAY+$ERv!|gmI{`zTlwz`~wGu0!29*qblC{tZ zg?Hc|hqx75p(6VDi(+(7RL2c1Zo>Q`w?*;XJdvR(H$HSkrIc{QJk}98)5rQ}q0Xp^ zI;;sObV2JV7v`OCL77q=Uv7+LxhqPe&hamwdm=yjcxv-iHEN`FcpN7x?G?zY>-OmU_x2VGfM~^YcnJe}V zs;SW7AXG|YThU+}ib^SXYW2lXWT4JZhRVWF5q0=bzcdW_QESZNjIqd>16!kvOlraAM`3U4ki(R^s7m31Yi%wg99fj(sH3RFSP$ohyT<{H{y`xbhH8(9w znuyetR^RtZ94ex;i!0R*`V5HG%X_t-ef-cn=76PO`{6O8y<@NTgOeY^#D&0+nmfbj zfO&H%Fvhu>N#HZ5s%7l)RNNnfUB5*qMFDgPOLZe~WlrrLHY3~}Lyut800Q0D z%|69uNC#qA(KV@*&8QhjV8ek+aVCZnw;;Fq@(Z%w=z*n2ZSB*!Ddu97%end7#)J3P z^GjNh_eSV(x4DVMNr%<@eP^2Dc(Xh|)yEDZP#ihEmQ~(WuN{n`sqU1Rt$=rx)^OhNfu0g*o?VB7(TtA7%I--JJ&aufWve75^;tkBpn*IgeQxa=iBipgh2m}o^Q#; z#ZU~3*T-}cXYes>z(gAV>*~La0bz_jSrUmME5WKen}MSU-1f+oiMNCASFdOc$G$ux zo=@Ws(cH>0VDLwuumQJYFkHTF`9M4fJ{dnvAaJN`$80vk5Q|}RaR-GsgHOhUiCjY# zPkp_ywwp^_Lww^fpzX)9*^F+Jh_pE&E7DjxiILDpUtoxN+Mx@ZaWIKMaQ^(S*lpP` z6+_(2{qM3FNy!A(+Bc=J88`)lwv+6HID^k=?^G@yXtHT4If~>XcNzv#9=)86I5v&I z9p?c$b{IcR$FRer$h2*=oAZ3P%$`ADRaSO?HerhVC2V64{JQ9Rax7^tq`m{H>wQmu@Uc3ow+f zSN&*WwD~RXNjFJV1_ATihTGzdhZ*Ux$?fB(;^p~ny|Iu$RPtI^@$&R6un2>Gc!zUr zM(#HRPE>WB#AeLS#PDSMsB|`C#9{*Zm2(!b8U3;_xYx>}L>b1MEN)Xd;@7a=n74%M z)OM-j=Ii;D7*+>rkBI_kG7EADyc+zXh#f>`E(Xn)+lI3l!FmFc{K9ZHL$wOSDD%q6 zY{tt~1WKJ`scgm#1J><6w)NS_aifEg`PAFd3CNEg?Xq9r;+(znIL&=rY@-G~`2R`! ze|*x`nXW|33pOVrZ+e0FBDw@q4y;!T&0BM$kkX)m zKF0cwEJ2ObVT60kE)-9#&$ANuAZKbd>{zo08FaL|aKFmE$cB2%wc1dIl4zKZ`*Ja!xzgkp9t-T8EPorDZY8g8J z9I~Os9*$41L)o+n?DwSL0y3x19q*pGjKXQTH4lTYqf#2?a-W>L$c8>TRXlA(Hq^Iq zW5^R!M8PVOD&x5_TMa6=Uch7q~{8!MCMs6M-{F!;0CAg3TC8j%7rO zGx%M0$`%5{HGlZC8NpjIWdBu2*bLheZZSP@!L)BUzL|A?cY0ZZ!O3;mYBuBOHX?0! z_0Cgc%63Lhx0PZz;j_oI55L#07OpL6d$?;;i0ilPwnXj4P)r6a5|4w=snb3JD}xKI z#2LTsONSYGk?!pBXUhoeKXWHSygZ+avi%q)#tbl>1=_WPk73~f?$k1;!jH|Eau7q; zZ4N%{I27ds3I>_Ie#O@=TU;Er6&Mm!$4&3<-i*VTQNf+?v**4vo!><}^=KuA$bezH z#e?7n*NT4_m&h0WY8zwV=M;tW0~fv*UR{z^TN$uKA+7E~@;iEWy?rpWF{OACIR-}0dv~rAApgo@tYBLCxN>BOLA*~LaR*d%%Xr/** - Nuclides are named by concatenating their Z and their A numbers. For - example, U235 is named 92235. Metastable nuclides are appended with an - '_m' and their metastable number. For example, the first excited isomer - of Am-242 is named 95242_m1. + Nuclides are named by concatenating their atomic symbol and mass number. For + example, 'U235' or 'Pu239'. Metastable nuclides are appended with an + '_m' and their metastable number. For example, 'Am242_m1' :Datasets: - **data** (*double[][][]*) -- The energy release coefficients. The first axis indexes the component type. The second axis specifies diff --git a/openmc/data/endf_utils.py b/openmc/data/endf_utils.py index bfd9ae5c8..1a77c60a5 100644 --- a/openmc/data/endf_utils.py +++ b/openmc/data/endf_utils.py @@ -12,9 +12,8 @@ import re def read_float(float_string): """Parse ENDF 6E11.0 formatted string into a float.""" assert len(float_string) == 11 - pattern = '([\s\\-]\d+\\.\d+)([\\+\\-]\d+)' - mantissa, exponent = re.match(pattern, float_string).groups() - return float(mantissa + 'e' + exponent) + pattern = r'([\s\-]\d+\.\d+)([\+\-]\d+)' + return float(re.sub(pattern, r'\1e\2', float_string)) def read_CONT_line(line): diff --git a/openmc/data/fission_energy.py b/openmc/data/fission_energy.py index 334231a9f..60cc43564 100644 --- a/openmc/data/fission_energy.py +++ b/openmc/data/fission_energy.py @@ -6,8 +6,9 @@ import h5py import numpy as np from numpy.polynomial.polynomial import Polynomial -from .function import Tabulated1D, Sum +from .data import ATOMIC_SYMBOL from .endf_utils import read_float, read_CONT_line, identify_nuclide +from .function import Tabulated1D, Sum import openmc.checkvalue as cv if sys.version_info[0] >= 3: @@ -70,11 +71,11 @@ def _extract_458_data(filename): labels = ('EFR', 'ENP', 'END', 'EGP', 'EGD', 'EB', 'ENU', 'ER', 'ET') # Associate each set of values and uncertainties with its label. - value = dict() - uncertainty = dict() - for i in range(len(labels)): - value[labels[i]] = data[2*i::18] - uncertainty[labels[i]] = data[2*i + 1::18] + value = {} + uncertainty = {} + for i, label in enumerate(labels): + value[label] = data[2*i::18] + uncertainty[label] = data[2*i + 1::18] # In ENDF/B-7.1, data for 2nd-order coefficients were mistakenly not # converted from MeV to eV. Check for this error and fix it if present. @@ -94,13 +95,6 @@ def _extract_458_data(filename): for coeffs in value.values(): coeffs[2] *= 1e-6 for coeffs in uncertainty.values(): coeffs[2] *= 1e-6 - # Perform the sanity check again... just in case. - for coeffs in value.values(): - second_order = coeffs[2] - if abs(second_order) * 1e12 > 1e8: - raise ValueError("Encountered a ludicrously large second-" - "order polynomial coefficient.") - # Convert eV to MeV. for coeffs in value.values(): for i in range(len(coeffs)): @@ -112,8 +106,8 @@ def _extract_458_data(filename): return value, uncertainty -def write_compact_458_library(endf_files, output_name=None, comment=None, - verbose=False): +def write_compact_458_library(endf_files, output_name='fission_Q_data.h5', + comment=None, verbose=False): """Read ENDF files, strip the MF=1 MT=458 data and write to small HDF5. Parameters @@ -130,7 +124,6 @@ def write_compact_458_library(endf_files, output_name=None, comment=None, """ # Open the output file. - if output_name is None: output_name = 'fission_Q_data.h5' out = h5py.File(output_name, 'w', libver='latest') # Write comments, if given. This commented out comment is the one used for @@ -179,7 +172,7 @@ def write_compact_458_library(endf_files, output_name=None, comment=None, value, uncertainty = data # Make a group for this isomer. - name = str(ident['Z']) + str(ident['A']) + name = ATOMIC_SYMBOL[ident['Z']] + str(ident['A']) if ident['LISO'] != 0: name += '_m' + str(ident['LISO']) nuclide_group = out.create_group(name) @@ -447,7 +440,7 @@ class FissionEnergyRelease(object): and p.emission_mode == 'prompt'] else: raise ValueError('IncidentNeutron data has no fission ' - 'reaction.') + 'reaction.') if len(nu_prompt) == 0: raise ValueError('Nu data is needed to compute fission energy ' 'release with the Sher-Beck format.') @@ -533,7 +526,7 @@ class FissionEnergyRelease(object): elif group.attrs['format'].decode() == 'Sher-Beck': obj.form = 'Sher-Beck' obj.prompt_neutrons = Tabulated1D.from_hdf5( - group['prompt_neutrons']) + group['prompt_neutrons']) else: raise ValueError('Unrecognized energy release format') @@ -565,14 +558,14 @@ class FissionEnergyRelease(object): components = [s.decode() for s in fin.attrs['component order']] - nuclide_name = str(incident_neutron.atomic_number) + nuclide_name = ATOMIC_SYMBOL[incident_neutron.atomic_number] nuclide_name += str(incident_neutron.mass_number) if incident_neutron.metastable != 0: nuclide_name += '_m' + str(incident_neutron.metastable) if nuclide_name not in fin: return None - data = {c : fin[nuclide_name + '/data'][i, 0, :] + data = {c: fin[nuclide_name + '/data'][i, 0, :] for i, c in enumerate(components)} return cls._from_dictionary(data, incident_neutron) diff --git a/src/tally.F90 b/src/tally.F90 index ab630ceb2..a949313bd 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -710,7 +710,7 @@ contains score = p % absorb_wgt * & nuc % reactions(nuc % index_fission(1)) % Q_value * & micro_xs(p % event_nuclide) % fission / & - micro_xs(p % event_nuclide) % absorption + micro_xs(p % event_nuclide) % absorption * flux end if end associate else @@ -724,7 +724,7 @@ contains score = p % last_wgt * & nuc % reactions(nuc % index_fission(1)) % Q_value * & micro_xs(p % event_nuclide) % fission / & - micro_xs(p % event_nuclide) % absorption + micro_xs(p % event_nuclide) % absorption * flux end if end associate end if @@ -785,7 +785,7 @@ contains score = p % absorb_wgt & * nuc % fission_q_prompt % evaluate(p % last_E) & * micro_xs(p % event_nuclide) % fission & - / micro_xs(p % event_nuclide) % absorption + / micro_xs(p % event_nuclide) % absorption * flux end if end associate else @@ -799,7 +799,7 @@ contains score = p % last_wgt & * nuc % fission_q_prompt % evaluate(p % last_E) & * micro_xs(p % event_nuclide) % fission & - / micro_xs(p % event_nuclide) % absorption + / micro_xs(p % event_nuclide) % absorption * flux end if end associate end if @@ -844,7 +844,7 @@ contains score = p % absorb_wgt & * nuc % fission_q_recov % evaluate(p % last_E) & * micro_xs(p % event_nuclide) % fission & - / micro_xs(p % event_nuclide) % absorption + / micro_xs(p % event_nuclide) % absorption * flux end if end associate else @@ -858,7 +858,7 @@ contains score = p % last_wgt & * nuc % fission_q_recov % evaluate(p % last_E) & * micro_xs(p % event_nuclide) % fission & - / micro_xs(p % event_nuclide) % absorption + / micro_xs(p % event_nuclide) % absorption * flux end if end associate end if From 75878d583d6b1e0ea158c62226d130238cd552d6 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 8 Aug 2016 10:06:02 -0500 Subject: [PATCH 11/33] Address #695 comments --- openmc/data/function.py | 2 +- openmc/data/product.py | 5 ++--- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/openmc/data/function.py b/openmc/data/function.py index e827e1283..798aa18cd 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -13,7 +13,7 @@ INTERPOLATION_SCHEME = {1: 'histogram', 2: 'linear-linear', 3: 'linear-log', class Function1D(object): """A function of one independent variable with HDF5 support.""" - __meta__class = ABCMeta + __metaclass__ = ABCMeta def __init__(self): pass diff --git a/openmc/data/product.py b/openmc/data/product.py index 116905f7a..eec47c956 100644 --- a/openmc/data/product.py +++ b/openmc/data/product.py @@ -35,7 +35,7 @@ class Product(object): yield represents particles from prompt and delayed sources. particle : str What particle the reaction product is. - yield_ : float or openmc.data.Tabulated1D or openmc.data.Polynomial + yield_ : openmc.data.Function1D Yield of secondary particle in the reaction. """ @@ -118,8 +118,7 @@ class Product(object): @yield_.setter def yield_(self, yield_): - cv.check_type('product yield', yield_, - (Tabulated1D, Polynomial)) + cv.check_type('product yield', yield_, Function1D) self._yield = yield_ def to_hdf5(self, group): From d86583c41683d87b9498a1a7095e49017f684f8c Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 8 Aug 2016 10:55:09 -0500 Subject: [PATCH 12/33] Infer Sher-Beck vs. Madland from data type --- docs/source/io_formats/nuclear_data.rst | 51 +++++++++++-------- openmc/data/fission_energy.py | 65 +++++++++---------------- src/nuclide_header.F90 | 14 ++---- tests/test_tallies/results_true.dat | 2 +- 4 files changed, 60 insertions(+), 72 deletions(-) diff --git a/docs/source/io_formats/nuclear_data.rst b/docs/source/io_formats/nuclear_data.rst index 7544ca1f5..2a16c02ce 100644 --- a/docs/source/io_formats/nuclear_data.rst +++ b/docs/source/io_formats/nuclear_data.rst @@ -57,24 +57,33 @@ Incident Neutron Data **//fission_energy_release/** -:Attributes: - **format** (*char[]*) -- The energy-dependence format. Either - 'Madland' or 'Sher-Beck' - -:Datasets: - **fragments** (*double[]*) -- Polynomial coefficients for energy - released in the form of fragments - - **prompt_neutrons** (*double[]* or :ref:`tabulated <1d_tabulated>`) - -- Energy released in the form of prompt neutrons. Polynomial if - the format is Madland or a table if Sher-Beck. - - **delayed_neutrons** (*double[]*) -- Polynomial coefficients for - energy released in the form of delayed neutrons - - **prompt_photons** (*double[]*) -- Polynomial coefficients for - energy released in the form of prompt photons - - **delayed_photons** (*double[]*) -- Polynomial coefficients for - energy released in the form of delayed photons - - **betas** (*double[]*) -- Polynomial coefficients for - energy released in the form of betas - - **neutrinos** (*double[]*) -- Polynomial coefficients for - energy released in the form of neutrinos +:Datasets: - **fragments** (:ref:`polynomial <1d_polynomial>`) -- Energy + released in the form of fragments as a function of incident + neutron energy. + - **prompt_neutrons** (:ref:`polynomial <1d_polynomial>` or + :ref:`tabulated <1d_tabulated>`) -- Energy released in the form of + prompt neutrons as a function of incident neutron energy. + - **delayed_neutrons** (:ref:`polynomial <1d_polynomial>`) -- Energy + released in the form of delayed neutrons as a function of incident + neutron energy. + - **prompt_photons** (:ref:`polynomial <1d_polynomial>`) -- Energy + released in the form of prompt photons as a function of incident + neutron energy. + - **delayed_photons** (:ref:`polynomial <1d_polynomial>`) -- Energy + released in the form of delayed photons as a function of incident + neutron energy. + - **betas** (:ref:`polynomial <1d_polynomial>`) -- Energy + released in the form of betas as a function of incident + neutron energy. + - **neutrinos** (:ref:`polynomial <1d_polynomial>`) -- Energy + released in the form of neutrinos as a function of incident + neutron energy. + - **q_prompt** (:ref:`polynomial <1d_polynomial>` or + :ref:`tabulated <1d_tabulated>`) -- The prompt fission Q-value + (fragments + prompt neutrons + prompt photons - incident energy) + - **q_recoverable** (:ref:`polynomial <1d_polynomial>` or + :ref:`tabulated <1d_tabulated>`) -- The recoverable fission Q-value + (Q_prompt + delayed neutrons + delayed photons + betas) ------------------------------- Thermal Neutron Scattering Data @@ -163,17 +172,19 @@ Tabulated :Object type: Dataset :Datatype: *double[2][]* :Description: x-values are listed first followed by corresponding y-values -:Attributes: - **type** (*char[]*) -- 'tabulated' +:Attributes: - **type** (*char[]*) -- 'Tabulated1D' - **breakpoints** (*int[]*) -- Region breakpoints - **interpolation** (*int[]*) -- Region interpolation codes +.. _1d_polynomial: + Polynomial ---------- :Object type: Dataset :Datatype: *double[]* :Description: Polynomial coefficients listed in order of increasing power -:Attributes: - **type** (*char[]*) -- 'polynomial' +:Attributes: - **type** (*char[]*) -- 'Polynomial' Coherent elastic scattering --------------------------- diff --git a/openmc/data/fission_energy.py b/openmc/data/fission_energy.py index 60cc43564..24743c19c 100644 --- a/openmc/data/fission_energy.py +++ b/openmc/data/fission_energy.py @@ -4,11 +4,10 @@ import sys import h5py import numpy as np -from numpy.polynomial.polynomial import Polynomial from .data import ATOMIC_SYMBOL from .endf_utils import read_float, read_CONT_line, identify_nuclide -from .function import Tabulated1D, Sum +from .function import Function1D, Tabulated1D, Polynomial, Sum import openmc.checkvalue as cv if sys.version_info[0] >= 3: @@ -262,9 +261,6 @@ class FissionEnergyRelease(object): q_total : Callable Function of energy that returns the total fission Q-value (total release - incident neutron energy). - form : str - Format used to compute the energy-dependence of the data. Either - 'Sher-Beck' or 'Madland'. """ def __init__(self): @@ -275,7 +271,6 @@ class FissionEnergyRelease(object): self._delayed_photons = None self._betas = None self._neutrinos = None - self._form = None @property def fragments(self): @@ -328,10 +323,6 @@ class FissionEnergyRelease(object): @property def q_total(self): return Sum([self.total, lambda E: -E]) - - @property - def form(self): - return self._form @fragments.setter def fragments(self, energy_release): @@ -368,11 +359,6 @@ class FissionEnergyRelease(object): cv.check_type('neutrinos', energy_release, Callable) self._neutrinos = energy_release - @form.setter - def form(self, form): - cv.check_value('format', form, ('Madland', 'Sher-Beck')) - self._form = form - @classmethod def _from_dictionary(cls, energy_release, incident_neutron): """Generate fission energy release data from a dictionary. @@ -401,7 +387,6 @@ class FissionEnergyRelease(object): # energy dependence. Otherwise, it is a polynomial. n_coeffs = len(energy_release['EFR']) if n_coeffs > 1: - out.form = 'Madland' out.fragments = Polynomial(energy_release['EFR']) out.prompt_neutrons = Polynomial(energy_release['ENP']) out.delayed_neutrons = Polynomial(energy_release['END']) @@ -410,12 +395,9 @@ class FissionEnergyRelease(object): out.betas = Polynomial(energy_release['EB']) out.neutrinos = Polynomial(energy_release['ENU']) else: - out.form = 'Sher-Beck' - - # EFR and ENP are energy independent. Polynomial is used because it - # has a __call__ attribute that handles Iterable inputs. The - # energy-dependence of END is unspecified in ENDF-102 so assume it - # is independent. + # EFR and ENP are energy independent. Use 0-order polynomials to + # make a constant function. The energy-dependence of END is + # unspecified in ENDF-102 so assume it is independent. out.fragments = Polynomial((energy_release['EFR'][0])) out.prompt_photons = Polynomial((energy_release['EGP'][0])) out.delayed_neutrons = Polynomial((energy_release['END'][0])) @@ -580,41 +562,40 @@ class FissionEnergyRelease(object): """ - group.create_dataset('fragments', data=self.fragments.coef) - group.create_dataset('delayed_neutrons', - data=self.delayed_neutrons.coef) - group.create_dataset('prompt_photons', - data=self.prompt_photons.coef) - group.create_dataset('delayed_photons', - data=self.delayed_photons.coef) - group.create_dataset('betas', data=self.betas.coef) - group.create_dataset('neutrinos', data=self.neutrinos.coef) - - if self.form == 'Madland': - group.attrs['format'] = np.string_('Madland') - group.create_dataset('prompt_neutrons', - data=self.prompt_neutrons.coef) - + self.fragments.to_hdf5(group, 'fragments') + self.prompt_neutrons.to_hdf5(group, 'prompt_neutrons') + self.delayed_neutrons.to_hdf5(group, 'delayed_neutrons') + self.prompt_photons.to_hdf5(group, 'prompt_photons') + self.delayed_photons.to_hdf5(group, 'delayed_photons') + self.betas.to_hdf5(group, 'betas') + self.neutrinos.to_hdf5(group, 'neutrinos') + + if isinstance(self.prompt_neutrons, Polynomial): + # Add the polynomials for the relevant components together. Use a + # Polynomial((0.0, -1.0)) to subtract incident energy. q_prompt = (self.fragments + self.prompt_neutrons + self.prompt_photons + Polynomial((0.0, -1.0))) - group.create_dataset('q_prompt', data=q_prompt.coef) + q_prompt.to_hdf5(group, 'q_prompt') q_recoverable = (self.fragments + self.prompt_neutrons + self.delayed_neutrons + self.prompt_photons + self.delayed_photons + self.betas + Polynomial((0.0, -1.0))) - group.create_dataset('q_recoverable', data=q_recoverable.coef) - elif self.form == 'Sher-Beck': - group.attrs['format'] = np.string_('Sher-Beck') - self.prompt_neutrons.to_hdf5(group, 'prompt_neutrons') + q_recoverable.to_hdf5(group, 'q_recoverable') + elif isinstance(self.prompt_neutrons, Tabulated1D): + # Make a Tabulated1D and evaluate the polynomial components at the + # table x points to get new y points. Subtract x from y to remove + # incident energy. q_prompt = deepcopy(self.prompt_neutrons) q_prompt.y += self.fragments(q_prompt.x) q_prompt.y += self.prompt_photons(q_prompt.x) + q_prompt.y -= q_prompt.x q_prompt.to_hdf5(group, 'q_prompt') q_recoverable = q_prompt q_recoverable.y += self.delayed_neutrons(q_recoverable.x) q_recoverable.y += self.delayed_photons(q_recoverable.x) q_recoverable.y += self.betas(q_recoverable.x) q_recoverable.to_hdf5(group, 'q_recoverable') + else: raise ValueError('Unrecognized energy release format') diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 4fcbf3af8..f32a93030 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -305,13 +305,13 @@ module nuclide_header hdf5_err) if (exists) then fer_group = open_group(group_id, 'fission_energy_release') - call read_attribute(temp, fer_group, 'format') - if (temp == 'Madland') then - ! The data uses the Madland format, i.e. polynomials + ! Check to see if this is polynomial or tabulated data + fer_dset = open_dataset(fer_group, 'q_prompt') + call read_attribute(temp, fer_dset, 'type') + if (temp == 'Polynomial') then ! Read the prompt Q-value allocate(Polynomial :: this % fission_q_prompt) - fer_dset = open_dataset(fer_group, 'q_prompt') call this % fission_q_prompt % from_hdf5(fer_dset) call close_dataset(fer_dset) @@ -320,13 +320,9 @@ module nuclide_header fer_dset = open_dataset(fer_group, 'q_recoverable') call this % fission_q_recov % from_hdf5(fer_dset) call close_dataset(fer_dset) - else if (temp == 'Sher-Beck') then - ! The data uses the Sher-Beck format. Python has handily converted this - ! format to Tabulated1Ds. - + else if (temp == 'Tabulated1D') then ! Read the prompt Q-value allocate(Tabulated1D :: this % fission_q_prompt) - fer_dset = open_dataset(fer_group, 'q_prompt') call this % fission_q_prompt % from_hdf5(fer_dset) call close_dataset(fer_dset) diff --git a/tests/test_tallies/results_true.dat b/tests/test_tallies/results_true.dat index 818d99a92..d018a65da 100644 --- a/tests/test_tallies/results_true.dat +++ b/tests/test_tallies/results_true.dat @@ -1 +1 @@ -a6e5480c66e6510687bf281983b3f387da45005bb08d8cd0aa629e53c5a42a1b2fa8d76345ad2df497d85f2b273fbc5edc36105a71a73c1107479ca0af8329f6 \ No newline at end of file +5a0f3f1ae244ada7d8c9f444d7a98c2589a9720e78174a276cf96162df6728bab6d534f136f65cb0386c3235eb7d5db47a0dd6504636ca77e7eb375e6978a84d \ No newline at end of file From b3c134e9bbef99e48261a18c8d665fdc58228d5a Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 8 Aug 2016 11:09:50 -0500 Subject: [PATCH 13/33] Update FissionEnergyRelease.from_hdf5 --- openmc/data/fission_energy.py | 23 +++++++---------------- 1 file changed, 7 insertions(+), 16 deletions(-) diff --git a/openmc/data/fission_energy.py b/openmc/data/fission_energy.py index 24743c19c..03d3ed84e 100644 --- a/openmc/data/fission_energy.py +++ b/openmc/data/fission_energy.py @@ -495,22 +495,13 @@ class FissionEnergyRelease(object): obj = cls() - obj.fragments = Polynomial(group['fragments'].value) - obj.delayed_neutrons = Polynomial(group['delayed_neutrons'].value) - obj.prompt_photons = Polynomial(group['prompt_photons'].value) - obj.delayed_photons = Polynomial(group['delayed_photons'].value) - obj.betas = Polynomial(group['betas'].value) - obj.neutrinos = Polynomial(group['neutrinos'].value) - - if group.attrs['format'].decode() == 'Madland': - obj.form = 'Madland' - obj.prompt_neutrons = Polynomial(group['prompt_neutrons'].value) - elif group.attrs['format'].decode() == 'Sher-Beck': - obj.form = 'Sher-Beck' - obj.prompt_neutrons = Tabulated1D.from_hdf5( - group['prompt_neutrons']) - else: - raise ValueError('Unrecognized energy release format') + obj.fragments = Function1D.from_hdf5(group['fragments']) + obj.prompt_neutrons = Function1D.from_hdf5(group['prompt_neutrons']) + obj.delayed_neutrons = Function1D.from_hdf5(group['delayed_neutrons']) + obj.prompt_photons = Function1D.from_hdf5(group['prompt_photons']) + obj.delayed_photons = Function1D.from_hdf5(group['delayed_photons']) + obj.betas = Function1D.from_hdf5(group['betas']) + obj.neutrinos = Function1D.from_hdf5(group['neutrinos']) return obj From 5b0ba57792aca84b8831d5f93fc9edb476a9b229 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 8 Aug 2016 11:23:00 -0500 Subject: [PATCH 14/33] Allow longer score names in statepoint files --- src/endf.F90 | 2 +- src/state_point.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/endf.F90 b/src/endf.F90 index 833082e1c..07094d5d9 100644 --- a/src/endf.F90 +++ b/src/endf.F90 @@ -14,7 +14,7 @@ contains pure function reaction_name(MT) result(string) integer, intent(in) :: MT - character(20) :: string + character(MAX_WORD_LEN) :: string select case (MT) ! Special reactions for tallies diff --git a/src/state_point.F90 b/src/state_point.F90 index abe034897..39532652f 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -51,7 +51,7 @@ contains integer(HID_T) :: file_id integer(HID_T) :: cmfd_group, tallies_group, tally_group, meshes_group, & mesh_group, filter_group, runtime_group - character(20), allocatable :: str_array(:) + character(MAX_WORD_LEN), allocatable :: str_array(:) character(MAX_FILE_LEN) :: filename type(RegularMesh), pointer :: meshp type(TallyObject), pointer :: tally From 7d23e09710387b0479cccf62eb1117f5c7e11b9f Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 8 Aug 2016 12:51:30 -0500 Subject: [PATCH 15/33] Update Travis data --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 30708de36..173af7602 100644 --- a/.travis.yml +++ b/.travis.yml @@ -42,7 +42,7 @@ install: true before_script: - if [[ ! -e $HOME/nndc_hdf5/cross_sections.xml ]]; then - wget https://anl.box.com/shared/static/6pwyfjnufam0sb96kqwwrve6vdn8m7u4.xz -O - | tar -C $HOME -xvJ; + wget https://anl.box.com/shared/static/dqkwdl7o4lauo91h3mgrn9qno6a3c8mp.xz -O - | tar -C $HOME -xvJ; fi - export OPENMC_CROSS_SECTIONS=$HOME/nndc_hdf5/cross_sections.xml From 651ee4240b24fa476bd8ba47606f9d89f9d317f7 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 10 Aug 2016 11:29:32 -0500 Subject: [PATCH 16/33] Address #698 comments --- docs/source/io_formats/fission_energy.rst | 4 ++-- openmc/data/function.py | 10 +++++----- scripts/openmc-ace-to-hdf5 | 2 +- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/docs/source/io_formats/fission_energy.rst b/docs/source/io_formats/fission_energy.rst index 1d3231ee3..768db56eb 100644 --- a/docs/source/io_formats/fission_energy.rst +++ b/docs/source/io_formats/fission_energy.rst @@ -9,8 +9,8 @@ ENDF-102_ for details). It gives the information needed to compute the energy carried away from fission reactions by each reaction product (e.g. fragment nuclei, neutrons) which depends on the incident neutron energy. OpenMC is distributed with one of these files under -openmc/data/fission_Q_data_endfb71.h5. More files of this format can be -created from ENDF files with the +data/fission_Q_data_endfb71.h5. More files of this format can be created from +ENDF files with the ``openmc.data.write_compact_458_library`` function. They can be read with the ``openmc.data.FissionEnergyRelease.from_compact_hdf5`` class method. diff --git a/openmc/data/function.py b/openmc/data/function.py index 798aa18cd..0d4fd3b1d 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -13,12 +13,11 @@ INTERPOLATION_SCHEME = {1: 'histogram', 2: 'linear-linear', 3: 'linear-log', class Function1D(object): """A function of one independent variable with HDF5 support.""" - __metaclass__ = ABCMeta - def __init__(self): pass - @abstractmethod - def __call__(self): pass + def __call__(self): + raise NotImplemented('Subclasses of Function1D should overwrite the ' + '__call__ and to_hdf5 methods') @abstractmethod def to_hdf5(self, group, name='xy'): @@ -32,7 +31,8 @@ class Function1D(object): Name of the dataset to create """ - pass + raise NotImplemented('Subclasses of Function1D should overwrite the ' + '__call__ and to_hdf5 methods') @classmethod def from_hdf5(cls, dataset): diff --git a/scripts/openmc-ace-to-hdf5 b/scripts/openmc-ace-to-hdf5 index 90031f0fa..2051de00b 100755 --- a/scripts/openmc-ace-to-hdf5 +++ b/scripts/openmc-ace-to-hdf5 @@ -124,7 +124,7 @@ for filename in ace_libraries: # Fission energy release data, if available if args.fission_energy_release is not None: fer = openmc.data.FissionEnergyRelease.from_compact_hdf5( - args.fission_energy_release, neutron) + args.fission_energy_release, neutron) if fer is not None: neutron.fission_energy = fer From 043c56f3384e738be9eba41ac1d8e5aa4dd17766 Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Thu, 11 Aug 2016 11:31:29 -0400 Subject: [PATCH 17/33] Now track OpenCG lattice names. Now strip hyphens on nuclide names --- openmc/nuclide.py | 10 +++++++++- openmc/opencg_compatible.py | 3 ++- 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/openmc/nuclide.py b/openmc/nuclide.py index 14609161f..961ab201f 100644 --- a/openmc/nuclide.py +++ b/openmc/nuclide.py @@ -1,5 +1,6 @@ from numbers import Integral import sys +import warnings from openmc.checkvalue import check_type @@ -98,7 +99,14 @@ class Nuclide(object): @name.setter def name(self, name): check_type('name', name, basestring) - self._name = name + + if '-' in name: + new_name = name.strip('-', '') + msg = 'OpenMC nuclide names follow the GND standard. Nuclide ' \ + '"{}" is being transformed to "{}".'.format(name, new_name) + warnings.warn(msg, DeprecationWarning) + + self._name = name.replace('-', '') @xs.setter def xs(self, xs): diff --git a/openmc/opencg_compatible.py b/openmc/opencg_compatible.py index 93a257f46..a3b35e055 100644 --- a/openmc/opencg_compatible.py +++ b/openmc/opencg_compatible.py @@ -911,6 +911,7 @@ def get_openmc_lattice(opencg_lattice): if lattice_id in OPENMC_LATTICES: return OPENMC_LATTICES[lattice_id] + name = opencg_lattice.name dimension = opencg_lattice.dimension width = opencg_lattice.width offset = opencg_lattice.offset @@ -941,7 +942,7 @@ def get_openmc_lattice(opencg_lattice): ((np.array(width, dtype=np.float64) * np.array(dimension, dtype=np.float64))) / -2.0 - openmc_lattice = openmc.RectLattice(lattice_id=lattice_id) + openmc_lattice = openmc.RectLattice(lattice_id=lattice_id, name=name) openmc_lattice.pitch = width openmc_lattice.universes = universe_array openmc_lattice.lower_left = lower_left From 8bc7afe91fe396bb79670221d6a235bf086dfc64 Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Thu, 11 Aug 2016 11:48:21 -0400 Subject: [PATCH 18/33] Fixed bug in Nuclide hyphen strip --- openmc/nuclide.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/nuclide.py b/openmc/nuclide.py index 961ab201f..6f7bd79b8 100644 --- a/openmc/nuclide.py +++ b/openmc/nuclide.py @@ -101,7 +101,7 @@ class Nuclide(object): check_type('name', name, basestring) if '-' in name: - new_name = name.strip('-', '') + new_name = name.replace('-', '') msg = 'OpenMC nuclide names follow the GND standard. Nuclide ' \ '"{}" is being transformed to "{}".'.format(name, new_name) warnings.warn(msg, DeprecationWarning) From 12c8672804b04391cd63d59ff75b3d45d10fb670 Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Thu, 11 Aug 2016 12:20:02 -0400 Subject: [PATCH 19/33] Reworded nuclide de-hyphen deprecation warning message --- openmc/nuclide.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openmc/nuclide.py b/openmc/nuclide.py index 6f7bd79b8..124c86a35 100644 --- a/openmc/nuclide.py +++ b/openmc/nuclide.py @@ -102,8 +102,8 @@ class Nuclide(object): if '-' in name: new_name = name.replace('-', '') - msg = 'OpenMC nuclide names follow the GND standard. Nuclide ' \ - '"{}" is being transformed to "{}".'.format(name, new_name) + msg = 'OpenMC nuclides follow the GND naming convention. Nuclide ' \ + '"{}" is being renamed as "{}".'.format(name, new_name) warnings.warn(msg, DeprecationWarning) self._name = name.replace('-', '') From fe63a59c5b192ec3962586c6173f2b0d4d3e34dc Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Thu, 11 Aug 2016 13:38:56 -0400 Subject: [PATCH 20/33] Moved get_thermal_name(...) from scripts to openmc.data.thermal --- openmc/data/thermal.py | 22 +++++++++++++++++++++- openmc/material.py | 3 ++- scripts/openmc-update-inputs | 21 ++------------------- 3 files changed, 25 insertions(+), 21 deletions(-) diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index 9de39f518..275836be7 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -40,6 +40,26 @@ _THERMAL_NAMES = {'al': 'c_Al27', 'al27': 'c_Al27', 'zrzrh': 'c_Zr_in_ZrH', 'zr-h': 'c_Zr_in_ZrH', 'zr/h': 'c_Zr_in_ZrH'} +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 + + class CoherentElastic(object): r"""Coherent elastic scattering data from a crystalline material @@ -394,4 +414,4 @@ class ThermalScattering(object): pairs = np.fromiter(map(lambda p: p[0], ace.pairs), int) table.zaids = pairs[np.nonzero(pairs)] - return table + return table \ No newline at end of file diff --git a/openmc/material.py b/openmc/material.py index d7ffd04e9..e1f3da97a 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -6,6 +6,7 @@ from xml.etree import ElementTree as ET import sys import openmc +import openmc.data import openmc.checkvalue as cv from openmc.clean_xml import sort_xml_elements, clean_xml_indentation @@ -484,7 +485,7 @@ class Material(object): 'non-string cross-section identifier "{1}"'.format(self._id, xs) raise ValueError(msg) - self._sab.append((name, xs)) + self._sab.append((openmc.data.get_thermal_name(name), xs)) def make_isotropic_in_lab(self): for nuclide, percent, percent_type in self._nuclides: diff --git a/scripts/openmc-update-inputs b/scripts/openmc-update-inputs index 5b5bf0978..613e4a224 100755 --- a/scripts/openmc-update-inputs +++ b/scripts/openmc-update-inputs @@ -253,23 +253,6 @@ 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 @@ -298,13 +281,13 @@ def update_materials(root): for sab in material.findall('sab'): if 'name' in sab.attrib: sabname = sab.attrib['name'] - sab.set('name', get_thermal_name(sabname)) + sab.set('name', openmc.data.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) + name_elem.text = openmc.data.get_thermal_name(sabname) was_updated = True return was_updated From b791450b751a294cd42904b0638cba6b4fe72b9c Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Thu, 11 Aug 2016 13:46:07 -0400 Subject: [PATCH 21/33] Added conditional to pass qualifying S(a,b) names through converter method --- openmc/data/thermal.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index 275836be7..8b6e35b61 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -43,7 +43,9 @@ _THERMAL_NAMES = {'al': 'c_Al27', 'al27': 'c_Al27', 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: + if name in _THERMAL_NAMES.values(): + return name + elif name.lower() in _THERMAL_NAMES: return _THERMAL_NAMES[name.lower()] else: # Make an educated guess?? This actually works well for @@ -52,7 +54,7 @@ def get_thermal_name(name): matches = get_close_matches( name.lower(), _THERMAL_NAMES.keys(), cutoff=0.5) if len(matches) > 0: - return _THERMAL_NAMES[matches[0]] + '.' + xs + return _THERMAL_NAMES[matches[0]] else: # OK, we give up. Just use the ACE name. return 'c_' + name From 11a2c31031eeabe29595ab335716739bc2c24cb1 Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Thu, 11 Aug 2016 22:12:47 -0400 Subject: [PATCH 22/33] Removed unreachable return in get_thermal_name(...) --- openmc/data/thermal.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index 8b6e35b61..cd44f21f9 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -59,8 +59,6 @@ def get_thermal_name(name): # OK, we give up. Just use the ACE name. return 'c_' + name - return name - class CoherentElastic(object): r"""Coherent elastic scattering data from a crystalline material @@ -416,4 +414,4 @@ class ThermalScattering(object): pairs = np.fromiter(map(lambda p: p[0], ace.pairs), int) table.zaids = pairs[np.nonzero(pairs)] - return table \ No newline at end of file + return table From 531a1726078a4b2d7abacf1752a04187a38805ed Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Thu, 11 Aug 2016 22:21:40 -0400 Subject: [PATCH 23/33] Now appending metatable suffix to nuclide names per suggestion by @paulromano --- openmc/nuclide.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/openmc/nuclide.py b/openmc/nuclide.py index 124c86a35..00ae25ac6 100644 --- a/openmc/nuclide.py +++ b/openmc/nuclide.py @@ -102,6 +102,10 @@ class Nuclide(object): if '-' in name: new_name = name.replace('-', '') + new_name = new_name.replace('Nat', '0') + if new_name.endswith('m'): + new_name = new_name[:-1] + '_m1' + msg = 'OpenMC nuclides follow the GND naming convention. Nuclide ' \ '"{}" is being renamed as "{}".'.format(name, new_name) warnings.warn(msg, DeprecationWarning) From cdb882648aff932644870f6bd96c45c3c0c032c7 Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Thu, 11 Aug 2016 22:27:53 -0400 Subject: [PATCH 24/33] Now printing a UserWarning in place of a DeprecationWarning for nuclide renaming --- openmc/nuclide.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/nuclide.py b/openmc/nuclide.py index 00ae25ac6..2c6998cd5 100644 --- a/openmc/nuclide.py +++ b/openmc/nuclide.py @@ -108,7 +108,7 @@ class Nuclide(object): msg = 'OpenMC nuclides follow the GND naming convention. Nuclide ' \ '"{}" is being renamed as "{}".'.format(name, new_name) - warnings.warn(msg, DeprecationWarning) + warnings.warn(msg) self._name = name.replace('-', '') From 284058319097c4777d0e3d0c7814b7ed6c4ef665 Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Thu, 11 Aug 2016 23:00:59 -0400 Subject: [PATCH 25/33] Fixed nuclide renaming bug --- openmc/nuclide.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/openmc/nuclide.py b/openmc/nuclide.py index 2c6998cd5..68f95beb9 100644 --- a/openmc/nuclide.py +++ b/openmc/nuclide.py @@ -99,19 +99,18 @@ class Nuclide(object): @name.setter def name(self, name): check_type('name', name, basestring) + self._name = name if '-' in name: - new_name = name.replace('-', '') - new_name = new_name.replace('Nat', '0') - if new_name.endswith('m'): - new_name = new_name[:-1] + '_m1' + self._name = name.replace('-', '') + self._name = self._name.replace('Nat', '0') + if self._name.endswith('m'): + self._name = self._name[:-1] + '_m1' msg = 'OpenMC nuclides follow the GND naming convention. Nuclide ' \ - '"{}" is being renamed as "{}".'.format(name, new_name) + '"{}" is being renamed as "{}".'.format(name, self._name) warnings.warn(msg) - self._name = name.replace('-', '') - @xs.setter def xs(self, xs): check_type('cross-section identifier', xs, basestring) From 6ecf89424e62d65f60f029d83b93645c58762753 Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Fri, 12 Aug 2016 10:02:39 -0400 Subject: [PATCH 26/33] Added warning message for S(a,b) table renaming --- openmc/material.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/openmc/material.py b/openmc/material.py index e1f3da97a..da173c92e 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -485,7 +485,13 @@ class Material(object): 'non-string cross-section identifier "{1}"'.format(self._id, xs) raise ValueError(msg) - self._sab.append((openmc.data.get_thermal_name(name), xs)) + new_name = openmc.data.get_thermal_name(name) + if new_name != name: + msg = 'OpenMC S(a,b) tables follow the GND naming convention. ' \ + 'Table "{}" is being renamed as "{}".'.format(name, new_name) + warnings.warn(msg) + + self._sab.append((new_name, xs)) def make_isotropic_in_lab(self): for nuclide, percent, percent_type in self._nuclides: From 5b219fd6aa36f58cc5313b259b6f098ccf35a2d1 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 12 Aug 2016 11:32:31 -0500 Subject: [PATCH 27/33] Fix doc build for metaclassed Polynomials --- docs/source/conf.py | 3 +++ openmc/data/function.py | 13 +++++-------- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/docs/source/conf.py b/docs/source/conf.py index 4aa000f38..1baea2b03 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -28,6 +28,9 @@ MOCK_MODULES = ['numpy', 'numpy.polynomial', 'numpy.polynomial.polynomial', 'h5py', 'pandas', 'opencg'] sys.modules.update((mod_name, MagicMock()) for mod_name in MOCK_MODULES) +import numpy as np +np.polynomial.Polynomial = MagicMock + # If extensions (or modules to document with autodoc) are in another directory, # add these directories to sys.path here. If the directory is relative to the diff --git a/openmc/data/function.py b/openmc/data/function.py index 768250259..88d27ce93 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -1,6 +1,7 @@ from abc import ABCMeta, abstractmethod from collections import Iterable, Callable from numbers import Real, Integral +from six import with_metaclass import numpy as np @@ -10,14 +11,11 @@ INTERPOLATION_SCHEME = {1: 'histogram', 2: 'linear-linear', 3: 'linear-log', 4: 'log-linear', 5: 'log-log'} -class Function1D(object): +class Function1D(with_metaclass(ABCMeta, object)): """A function of one independent variable with HDF5 support.""" - def __init__(self): pass - - def __call__(self): - raise NotImplemented('Subclasses of Function1D should overwrite the ' - '__call__ and to_hdf5 methods') + @abstractmethod + def __call__(self): pass @abstractmethod def to_hdf5(self, group, name='xy'): @@ -31,8 +29,7 @@ class Function1D(object): Name of the dataset to create """ - raise NotImplemented('Subclasses of Function1D should overwrite the ' - '__call__ and to_hdf5 methods') + pass @classmethod def from_hdf5(cls, dataset): From f4e0525a053741fda2e507aa2ca7954e7d1d73d9 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 12 Aug 2016 14:16:46 -0500 Subject: [PATCH 28/33] Remove six dependency --- openmc/data/function.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/openmc/data/function.py b/openmc/data/function.py index 88d27ce93..a7397e785 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -1,7 +1,6 @@ from abc import ABCMeta, abstractmethod from collections import Iterable, Callable from numbers import Real, Integral -from six import with_metaclass import numpy as np @@ -11,9 +10,11 @@ INTERPOLATION_SCHEME = {1: 'histogram', 2: 'linear-linear', 3: 'linear-log', 4: 'log-linear', 5: 'log-log'} -class Function1D(with_metaclass(ABCMeta, object)): +class Function1D(object): """A function of one independent variable with HDF5 support.""" + __metaclass__ = ABCMeta + @abstractmethod def __call__(self): pass From c5df6ce146abeee0d83447aa7a1deccf354b9ade Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 12 Aug 2016 16:40:13 -0500 Subject: [PATCH 29/33] Fix mesh filter max iterator check --- src/tally_filter.F90 | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index c0e1c8853..67d045284 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -295,8 +295,12 @@ contains search_iter = 0 do while (any(ijk0(:m % n_dimension) < 1) & .or. any(ijk0(:m % n_dimension) > m % dimension)) - if (search_iter == MAX_SEARCH_ITER) call fatal_error("Failed to & - &find a mesh intersection on a tally mesh filter.") + if (search_iter == MAX_SEARCH_ITER) then + call warning("Failed to find a mesh intersection on a tally mesh & + &filter.") + next_bin = NO_BIN_FOUND + return + end if do j = 1, m % n_dimension if (abs(uvw(j)) < FP_PRECISION) then @@ -315,6 +319,8 @@ contains else ijk0(j) = ijk0(j) - 1 end if + + search_iter = search_iter + 1 end do distance = d(j) xyz0 = xyz0 + distance * uvw From aca4578f46327e3978f95b9dcfd99d2610b2fff6 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 12 Aug 2016 17:57:54 -0500 Subject: [PATCH 30/33] Fix Sphinx bullet list error --- docs/source/io_formats/fission_energy.rst | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/source/io_formats/fission_energy.rst b/docs/source/io_formats/fission_energy.rst index 768db56eb..f80db4569 100644 --- a/docs/source/io_formats/fission_energy.rst +++ b/docs/source/io_formats/fission_energy.rst @@ -26,7 +26,8 @@ ENDF files with the example, 'U235' or 'Pu239'. Metastable nuclides are appended with an '_m' and their metastable number. For example, 'Am242_m1' -:Datasets: - **data** (*double[][][]*) -- The energy release coefficients. The +:Datasets: + - **data** (*double[][][]*) -- The energy release coefficients. The first axis indexes the component type. The second axis specifies values or uncertainties. The third axis indexes the polynomial order. If the data uses the Sher-Beck format, then the last axis From 328c6a0e886d0011f3e88b03a53028e7119d0250 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 13 Aug 2016 17:40:28 -0500 Subject: [PATCH 31/33] Allow estimator to be set for MGXS --- openmc/mgxs/mgxs.py | 74 ++++++++++++++++++++++----------------------- 1 file changed, 37 insertions(+), 37 deletions(-) diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index 7aa12d041..b94235155 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -47,6 +47,10 @@ DOMAIN_TYPES = ['cell', 'material', 'mesh'] +ESTIMATOR_TYPES = ['tracklength', + 'collision', + 'analog'] + # Supported domain classes _DOMAINS = (openmc.Cell, openmc.Universe, @@ -102,7 +106,7 @@ class MGXS(object): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'analog'} + estimator : {'tracklength', 'collision', 'analog'} The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section @@ -149,6 +153,7 @@ class MGXS(object): self._rxn_type = None self._by_nuclide = None self._nuclides = None + self._estimator = 'tracklength' self._domain = None self._domain_type = None self._energy_groups = None @@ -250,7 +255,7 @@ class MGXS(object): @property def estimator(self): - return 'tracklength' + return self._estimator @property def tallies(self): @@ -368,6 +373,11 @@ class MGXS(object): cv.check_iterable_type('nuclides', nuclides, basestring) self._nuclides = nuclides + @estimator.setter + def estimator(self, estimator): + cv.check_value('estimator', estimator, ESTIMATOR_TYPES) + self._estimator = estimator + @domain.setter def domain(self, domain): cv.check_type('domain', domain, _DOMAINS) @@ -1643,7 +1653,7 @@ class MatrixMGXS(MGXS): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'analog'} + estimator : {'tracklength', 'collision', 'analog'} The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section @@ -1692,10 +1702,6 @@ class MatrixMGXS(MGXS): return [[energy], [energy, energyout]] - @property - def estimator(self): - return 'analog' - def get_xs(self, in_groups='all', out_groups='all', subdomains='all', nuclides='all', xs_type='macro', order_groups='increasing', @@ -2072,7 +2078,7 @@ class TotalXS(MGXS): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'analog'} + estimator : {'tracklength', 'collision', 'analog'} The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys @@ -2190,7 +2196,7 @@ class TransportXS(MGXS): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'analog'} + estimator : {'tracklength', 'collision', 'analog'} The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys @@ -2233,6 +2239,7 @@ class TransportXS(MGXS): super(TransportXS, self).__init__(domain, domain_type, groups, by_nuclide, name) self._rxn_type = 'transport' + self._estimator = 'analog' @property def scores(self): @@ -2245,10 +2252,6 @@ class TransportXS(MGXS): energyout_filter = openmc.Filter('energyout', group_edges) return [[energy_filter], [energy_filter], [energyout_filter]] - @property - def estimator(self): - return 'analog' - @property def rxn_rate_tally(self): if self._rxn_rate_tally is None: @@ -2320,7 +2323,7 @@ class NuTransportXS(TransportXS): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'analog'} + estimator : {'tracklength', 'collision', 'analog'} The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys @@ -2441,7 +2444,7 @@ class AbsorptionXS(MGXS): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'analog'} + estimator : {'tracklength', 'collision', 'analog'} The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys @@ -2557,7 +2560,7 @@ class CaptureXS(MGXS): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'analog'} + estimator : {'tracklength', 'collision', 'analog'} The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys @@ -2679,7 +2682,7 @@ class FissionXS(MGXS): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'analog'} + estimator : {'tracklength', 'collision', 'analog'} The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys @@ -2790,7 +2793,7 @@ class NuFissionXS(MGXS): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'analog'} + estimator : {'tracklength', 'collision', 'analog'} The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys @@ -2906,7 +2909,7 @@ class KappaFissionXS(MGXS): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'analog'} + estimator : {'tracklength', 'collision', 'analog'} The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys @@ -3019,7 +3022,7 @@ class ScatterXS(MGXS): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'analog'} + estimator : {'tracklength', 'collision', 'analog'} The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys @@ -3134,7 +3137,7 @@ class NuScatterXS(MGXS): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'analog'} + estimator : {'tracklength', 'collision', 'analog'} The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys @@ -3177,10 +3180,7 @@ class NuScatterXS(MGXS): super(NuScatterXS, self).__init__(domain, domain_type, groups, by_nuclide, name) self._rxn_type = 'nu-scatter' - - @property - def estimator(self): - return 'analog' + self._estimator = 'analog' class ScatterMatrixXS(MatrixMGXS): @@ -3268,7 +3268,7 @@ class ScatterMatrixXS(MatrixMGXS): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'analog'} + estimator : {'tracklength', 'collision', 'analog'} The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys @@ -3314,6 +3314,7 @@ class ScatterMatrixXS(MatrixMGXS): self._correction = 'P0' self._legendre_order = 0 self._hdf5_key = 'scatter matrix' + self._estimator = 'analog' def __deepcopy__(self, memo): clone = super(ScatterMatrixXS, self).__deepcopy__(memo) @@ -3929,7 +3930,7 @@ class NuScatterMatrixXS(ScatterMatrixXS): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'analog'} + estimator : {'tracklength', 'collision', 'analog'} The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys @@ -4051,7 +4052,7 @@ class MultiplicityMatrixXS(MatrixMGXS): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'analog'} + estimator : {'tracklength', 'collision', 'analog'} The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys @@ -4094,6 +4095,7 @@ class MultiplicityMatrixXS(MatrixMGXS): super(MultiplicityMatrixXS, self).__init__(domain, domain_type, groups, by_nuclide, name) self._rxn_type = 'multiplicity matrix' + self._estimator = 'analog' @property def scores(self): @@ -4198,7 +4200,7 @@ class NuFissionMatrixXS(MatrixMGXS): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'analog'} + estimator : {'tracklength', 'collision', 'analog'} The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys @@ -4242,6 +4244,7 @@ class NuFissionMatrixXS(MatrixMGXS): groups, by_nuclide, name) self._rxn_type = 'nu-fission' self._hdf5_key = 'nu-fission matrix' + self._estimator = 'analog' class Chi(MGXS): @@ -4313,7 +4316,7 @@ class Chi(MGXS): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'analog'} + estimator : {'tracklength', 'collision', 'analog'} The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys @@ -4355,6 +4358,7 @@ class Chi(MGXS): groups=None, by_nuclide=False, name=''): super(Chi, self).__init__(domain, domain_type, groups, by_nuclide, name) self._rxn_type = 'chi' + self._estimator = 'analog' @property def scores(self): @@ -4372,10 +4376,6 @@ class Chi(MGXS): def tally_keys(self): return ['nu-fission-in', 'nu-fission-out'] - @property - def estimator(self): - return 'analog' - @property def rxn_rate_tally(self): if self._rxn_rate_tally is None: @@ -4803,7 +4803,7 @@ class ChiPrompt(Chi): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'analog'} + estimator : {'tracklength', 'collision', 'analog'} The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys @@ -4918,7 +4918,7 @@ class InverseVelocity(MGXS): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'analog'} + estimator : {'tracklength', 'collision', 'analog'} The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys @@ -5052,7 +5052,7 @@ class PromptNuFissionXS(MGXS): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'analog'} + estimator : {'tracklength', 'collision', 'analog'} The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys From 8e0f3f4f79fedcbd3892a0164d336f6315f97f45 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 15 Aug 2016 06:47:19 -0500 Subject: [PATCH 32/33] Respond to @wbinventor comments on #704 --- openmc/__init__.py | 2 +- openmc/mgxs/mgxs.py | 32 +++++++++++++++++--------------- openmc/tallies.py | 6 ++++-- 3 files changed, 22 insertions(+), 18 deletions(-) diff --git a/openmc/__init__.py b/openmc/__init__.py index 026ccce11..a0492ee40 100644 --- a/openmc/__init__.py +++ b/openmc/__init__.py @@ -13,10 +13,10 @@ from openmc.settings import * from openmc.surface import * from openmc.universe import * from openmc.mesh import * -from openmc.mgxs_library import * from openmc.filter import * from openmc.trigger import * from openmc.tallies import * +from openmc.mgxs_library import * from openmc.cmfd import * from openmc.executor import * from openmc.statepoint import * diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index b94235155..d68dc508f 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -13,6 +13,7 @@ import numpy as np import openmc import openmc.checkvalue as cv +from openmc.tallies import ESTIMATOR_TYPES from openmc.mgxs import EnergyGroups if sys.version_info[0] >= 3: @@ -39,7 +40,6 @@ MGXS_TYPES = ['total', 'inverse-velocity', 'prompt-nu-fission'] - # Supported domain types DOMAIN_TYPES = ['cell', 'distribcell', @@ -47,10 +47,6 @@ DOMAIN_TYPES = ['cell', 'material', 'mesh'] -ESTIMATOR_TYPES = ['tracklength', - 'collision', - 'analog'] - # Supported domain classes _DOMAINS = (openmc.Cell, openmc.Universe, @@ -148,7 +144,6 @@ class MGXS(object): def __init__(self, domain=None, domain_type=None, energy_groups=None, by_nuclide=False, name=''): - self._name = '' self._rxn_type = None self._by_nuclide = None @@ -165,6 +160,7 @@ class MGXS(object): self._loaded_sp = False self._derived = False self._hdf5_key = None + self._valid_estimators = ESTIMATOR_TYPES self.name = name self.by_nuclide = by_nuclide @@ -375,7 +371,7 @@ class MGXS(object): @estimator.setter def estimator(self, estimator): - cv.check_value('estimator', estimator, ESTIMATOR_TYPES) + cv.check_value('estimator', estimator, self._valid_estimators) self._estimator = estimator @domain.setter @@ -2196,7 +2192,7 @@ class TransportXS(MGXS): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'collision', 'analog'} + estimator : 'analog' The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys @@ -2240,6 +2236,7 @@ class TransportXS(MGXS): groups, by_nuclide, name) self._rxn_type = 'transport' self._estimator = 'analog' + self._valid_estimators = ['analog'] @property def scores(self): @@ -2323,7 +2320,7 @@ class NuTransportXS(TransportXS): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'collision', 'analog'} + estimator : 'analog' The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys @@ -3181,6 +3178,7 @@ class NuScatterXS(MGXS): groups, by_nuclide, name) self._rxn_type = 'nu-scatter' self._estimator = 'analog' + self._valid_estimators = ['analog'] class ScatterMatrixXS(MatrixMGXS): @@ -3268,7 +3266,7 @@ class ScatterMatrixXS(MatrixMGXS): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'collision', 'analog'} + estimator : 'analog' The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys @@ -3315,6 +3313,7 @@ class ScatterMatrixXS(MatrixMGXS): self._legendre_order = 0 self._hdf5_key = 'scatter matrix' self._estimator = 'analog' + self._valid_estimators = ['analog'] def __deepcopy__(self, memo): clone = super(ScatterMatrixXS, self).__deepcopy__(memo) @@ -3930,7 +3929,7 @@ class NuScatterMatrixXS(ScatterMatrixXS): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'collision', 'analog'} + estimator : 'analog' The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys @@ -4052,7 +4051,7 @@ class MultiplicityMatrixXS(MatrixMGXS): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'collision', 'analog'} + estimator : 'analog' The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys @@ -4096,6 +4095,7 @@ class MultiplicityMatrixXS(MatrixMGXS): by_nuclide, name) self._rxn_type = 'multiplicity matrix' self._estimator = 'analog' + self._valid_estimators = ['analog'] @property def scores(self): @@ -4200,7 +4200,7 @@ class NuFissionMatrixXS(MatrixMGXS): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'collision', 'analog'} + estimator : 'analog' The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys @@ -4245,6 +4245,7 @@ class NuFissionMatrixXS(MatrixMGXS): self._rxn_type = 'nu-fission' self._hdf5_key = 'nu-fission matrix' self._estimator = 'analog' + self._valid_estimators = ['analog'] class Chi(MGXS): @@ -4316,7 +4317,7 @@ class Chi(MGXS): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'collision', 'analog'} + estimator : 'analog' The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys @@ -4359,6 +4360,7 @@ class Chi(MGXS): super(Chi, self).__init__(domain, domain_type, groups, by_nuclide, name) self._rxn_type = 'chi' self._estimator = 'analog' + self._valid_estimators = ['analog'] @property def scores(self): @@ -4803,7 +4805,7 @@ class ChiPrompt(Chi): tally_keys : list of str The keys into the tallies dictionary for each tally used to compute the multi-group cross section - estimator : {'tracklength', 'collision', 'analog'} + estimator : 'analog' The tally estimator used to compute the multi-group cross section tallies : collections.OrderedDict OpenMC tallies needed to compute the multi-group cross section. The keys diff --git a/openmc/tallies.py b/openmc/tallies.py index c68b0faac..f645fa162 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -40,6 +40,9 @@ _SCORE_CLASSES = (basestring, CrossScore, AggregateScore) _NUCLIDE_CLASSES = (basestring, Nuclide, CrossNuclide, AggregateNuclide) _FILTER_CLASSES = (Filter, CrossFilter, AggregateFilter) +# Valid types of estimators +ESTIMATOR_TYPES = ['tracklength', 'collision', 'analog'] + def reset_auto_tally_id(): """Reset counter for auto-generated tally IDs.""" @@ -387,8 +390,7 @@ class Tally(object): @estimator.setter def estimator(self, estimator): - cv.check_value('estimator', estimator, - ['analog', 'tracklength', 'collision']) + cv.check_value('estimator', estimator, ESTIMATOR_TYPES) self._estimator = estimator @triggers.setter From 6be762858e14b9e1db69cf1366fbde004dc4cd51 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 16 Aug 2016 09:33:26 -0500 Subject: [PATCH 33/33] Allow estimator to be set for mgxs.Library --- openmc/mgxs/library.py | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index ee11d0ef6..5685681e4 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -11,6 +11,7 @@ import numpy as np import openmc import openmc.mgxs import openmc.checkvalue as cv +from openmc.tallies import ESTIMATOR_TYPES if sys.version_info[0] >= 3: @@ -39,7 +40,7 @@ class Library(object): mgxs_types : Iterable of str The types of cross sections in the library (e.g., ['total', 'scatter']) name : str, optional - Name of the multi-group cross section. library Used as a label to + Name of the multi-group cross section library. Used as a label to identify tallies in OpenMC 'tallies.xml' file. Attributes @@ -64,6 +65,9 @@ class Library(object): The highest legendre moment in the scattering matrices (default is 0) energy_groups : openmc.mgxs.EnergyGroups Energy group structure for energy condensation + estimator : str or None + The tally estimator used to compute multi-group cross sections. If None, + the default for each MGXS type is used. tally_trigger : openmc.Trigger An (optional) tally precision trigger given to each tally used to compute the cross section @@ -102,6 +106,7 @@ class Library(object): self._sp_filename = None self._keff = None self._sparse = False + self._estimator = None self.name = name self.openmc_geometry = openmc_geometry @@ -206,6 +211,10 @@ class Library(object): def tally_trigger(self): return self._tally_trigger + @property + def estimator(self): + return self._estimator + @property def num_groups(self): return self.energy_groups.num_groups @@ -327,6 +336,11 @@ class Library(object): cv.check_type('tally trigger', tally_trigger, openmc.Trigger) self._tally_trigger = tally_trigger + @estimator.setter + def estimator(self, estimator): + cv.check_value('estimator', estimator, ESTIMATOR_TYPES) + self._estimator = estimator + @sparse.setter def sparse(self, sparse): """Convert tally data from NumPy arrays to SciPy list of lists (LIL) @@ -368,9 +382,11 @@ class Library(object): mgxs.domain_type = self.domain_type mgxs.energy_groups = self.energy_groups mgxs.by_nuclide = self.by_nuclide + if self.estimator is not None: + mgxs.estimator = self.estimator # If a tally trigger was specified, add it to the MGXS - if self.tally_trigger: + if self.tally_trigger is not None: mgxs.tally_trigger = self.tally_trigger # Specify whether to use a transport ('P0') correction