From 9d2f804ca09ca6591c378c85a66315c099c75101 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 19 Jul 2018 10:57:19 -0500 Subject: [PATCH] Refactor how recoverable/prompt/total are calculated using sum_functions --- openmc/data/fission_energy.py | 122 ++++++++-------------------------- openmc/data/function.py | 42 ++++++++++++ 2 files changed, 69 insertions(+), 95 deletions(-) diff --git a/openmc/data/fission_energy.py b/openmc/data/fission_energy.py index 1a3ddbf90..0f30ecfea 100644 --- a/openmc/data/fission_energy.py +++ b/openmc/data/fission_energy.py @@ -6,25 +6,18 @@ import sys import h5py import numpy as np -from .data import ATOMIC_SYMBOL, EV_PER_MEV +from .data import EV_PER_MEV from .endf import get_cont_record, get_list_record, get_tab1_record, Evaluation -from .function import Function1D, Tabulated1D, Polynomial, Sum +from .function import Function1D, Tabulated1D, Polynomial, sum_functions import openmc.checkvalue as cv from openmc.mixin import EqualityMixin -_LABELS = ('EFR', 'ENP', 'END', 'EGP', 'EGD', 'EB', 'ENU', 'ER', 'ET') -_NAMES = { - 'EFR': 'fragments', - 'ENP': 'prompt_neutrons', - 'END': 'delayed_neutrons', - 'EGP': 'prompt_photons', - 'EGD': 'delayed_photons', - 'EB': 'betas', - 'ENU': 'neutrinos', - 'ER': 'recoverable', - 'ET': 'total' -} +_NAMES = ( + 'fragments', 'prompt_neutrons', 'delayed_neutrons', + 'prompt_photons', 'delayed_photons', 'betas', + 'neutrinos', 'recoverable', 'total' +) class FissionEnergyRelease(EqualityMixin): @@ -141,27 +134,33 @@ class FissionEnergyRelease(EqualityMixin): @property def recoverable(self): - return Sum([self.fragments, self.prompt_neutrons, self.delayed_neutrons, - self.prompt_photons, self.delayed_photons, self.betas]) + components = ['fragments', 'prompt_neutrons', 'delayed_neutrons', + 'prompt_photons', 'delayed_photons', 'betas'] + return sum_functions(getattr(self, c) for c in components) @property def total(self): - return Sum([self.fragments, self.prompt_neutrons, self.delayed_neutrons, - self.prompt_photons, self.delayed_photons, self.betas, - self.neutrinos]) + components = ['fragments', 'prompt_neutrons', 'delayed_neutrons', + 'prompt_photons', 'delayed_photons', 'betas', + 'neutrinos'] + return sum_functions(getattr(self, c) for c in components) @property def q_prompt(self): - return Sum([self.fragments, self.prompt_neutrons, self.prompt_photons, - lambda E: -E]) + # Use a polynomial to subtract incident energy. + funcs = [self.fragments, self.prompt_neutrons, self.prompt_photons, + Polynomial((0.0, -1.0))] + return sum_functions(funcs) @property def q_recoverable(self): - return Sum([self.recoverable, lambda E: -E]) + # Use a polynomial to subtract incident energy. + return sum_functions([self.recoverable, Polynomial((0.0, -1.0))]) @property def q_total(self): - return Sum([self.total, lambda E: -E]) + # Use a polynomial to subtract incident energy. + return sum_functions([self.total, Polynomial((0.0, -1.0))]) @fragments.setter def fragments(self, energy_release): @@ -247,9 +246,8 @@ class FissionEnergyRelease(EqualityMixin): # Associate each set of values and uncertainties with its label. functions = {} - for i, label in enumerate(_LABELS): + for i, name in enumerate(_NAMES): coeffs = data[2*i::18] - name = _NAMES[label] # Ignore recoverable and total since we recalculate those directly if name in ('recoverable', 'total'): @@ -327,13 +325,12 @@ class FissionEnergyRelease(EqualityMixin): # Determine which component it is ifc = items[3] - name = _NAMES[_LABELS[ifc - 1]] + name = _NAMES[ifc - 1] # Replace value in dictionary functions[name] = eifc # Build the object - print(functions) return cls(**functions) @classmethod @@ -361,44 +358,7 @@ class FissionEnergyRelease(EqualityMixin): neutrinos = Function1D.from_hdf5(group['neutrinos']) return cls(fragments, prompt_neutrons, delayed_neutrons, prompt_photons, - prompt_photons, delayed_photons, betas, neutrinos) - - @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 - :func:`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 = 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, :] - for i, c in enumerate(components)} - - return cls._from_dictionary(data, incident_neutron) + delayed_photons, betas, neutrinos) def to_hdf5(self, group): """Write energy release data to an HDF5 group @@ -417,33 +377,5 @@ class FissionEnergyRelease(EqualityMixin): 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))) - 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))) - 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') + self.q_prompt.to_hdf5(group, 'q_prompt') + self.q_recoverable.to_hdf5(group, 'q_recoverable') diff --git a/openmc/data/function.py b/openmc/data/function.py index 3d09e44fc..f7636f25c 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -1,5 +1,7 @@ from abc import ABCMeta, abstractmethod from collections.abc import Iterable, Callable +from functools import reduce +from itertools import zip_longest from numbers import Real, Integral import numpy as np @@ -13,6 +15,46 @@ INTERPOLATION_SCHEME = {1: 'histogram', 2: 'linear-linear', 3: 'linear-log', 4: 'log-linear', 5: 'log-log'} +def sum_functions(funcs): + """Add tabulated/polynomials functions together + + Parameters + ---------- + funcs : list of Function1D + Functions to add + + Returns + ------- + Function1D + Sum of polynomial/tabulated functions + + """ + # Copy so we can iterate multiple times + funcs = list(funcs) + + # Get x values for all tabulated components + xs = [] + for f in funcs: + if isinstance(f, Tabulated1D): + xs.append(f.x) + if not np.all(f.interpolation == 2): + raise ValueError('Only linear-linear tabulated functions ' + 'can be combined') + + if xs: + # Take the union of all energies (sorted) + x = reduce(np.union1d, xs) + + # Evaluate each function and add together + y = sum(f(x) for f in funcs) + return Tabulated1D(x, y) + else: + # If no tabulated functions are present, we need to combine the + # polynomials by adding their coefficients + coeffs = [sum(x) for x in zip_longest(*funcs, fillvalue=0.0)] + return Polynomial(coeffs) + + class Function1D(EqualityMixin, metaclass=ABCMeta): """A function of one independent variable with HDF5 support.""" @abstractmethod