diff --git a/openmc/filter.py b/openmc/filter.py index 88a88f99d..b9e4c4286 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1436,3 +1436,79 @@ class DelayedGroupFilter(IntegralFilter): filter's bins. """ + + +class LinLinEnergyFilter(Filter): + """ + + Parameters + ---------- + energy : Iterable of Real + A grid of energy values in eV. + y : iterable of Real + A grid of interpolant values in eV. + + Attributes + ---------- + energy : Iterable of Real + A grid of energy values in eV. + y : iterable of Real + A grid of interpolant values in eV. + num_bins : Integral + The number of filter bins (always 1 for this filter) + stride : Integral + The number of filter, nuclide and score bins within each of this + filter's bins. + + """ + + def __init__(self, energy, y): + self.energy = energy + self.y = y + self._stride = None + + @property + def energy(self): + return self._energy + + @property + def y(self): + return self._y + + @energy.setter + def energy(self, energy): + # Make sure the energy grid is a numpy array. + energy = np.array(energy) + + # If the grid is 0D numpy array, promote to 1D. + if energy.shape == (): + energy.shape = (1,) + + # Make sure the values are Real and positive. + cv.check_type('filter energy grid', energy, Iterable, Real) + for E in energy: + cv.check_greater_than('filter energy grid', E, 0, equality=True) + + self._energy = energy + + @y.setter + def y(self, y): + # Make sure the values are in a numpy array. + y = np.array(y) + + # If the array is 0D, promote to 1D. + if y.shape == (): + y.shape = (1,) + + # Make sure the values are Real. + cv.check_type('filter interpolant values', y, Iterable, Real) + + self._y = y + + def to_xml(self): + """Return XML Element representing the Filter.""" + element = ET.Element('filter') + element.set('type', self.short_name.lower()) + element.set('energy', ' '.join(str(e) for e in self.energy)) + element.set('y', ' '.join(str(y) for y in self.y)) + return element diff --git a/src/constants.F90 b/src/constants.F90 index 044ea1c6b..bd9a0508b 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -342,7 +342,7 @@ module constants integer, parameter :: NO_BIN_FOUND = -1 ! Tally filter and map types - integer, parameter :: N_FILTER_TYPES = 13 + integer, parameter :: N_FILTER_TYPES = 14 integer, parameter :: & FILTER_UNIVERSE = 1, & FILTER_MATERIAL = 2, & @@ -356,7 +356,8 @@ module constants FILTER_MU = 10, & FILTER_POLAR = 11, & FILTER_AZIMUTHAL = 12, & - FILTER_DELAYEDGROUP = 13 + FILTER_DELAYEDGROUP = 13, & + FILTER_LINLINENERGY = 14 ! Mesh types integer, parameter :: & diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 1bbf68f48..9279a48b8 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2973,18 +2973,21 @@ contains temp_str = to_lower(temp_str) ! Determine number of bins - if (check_for_node(node_filt, "bins")) then - if (temp_str == 'energy' .or. temp_str == 'energyout' .or. & - temp_str == 'mu' .or. temp_str == 'polar' .or. & - temp_str == 'azimuthal') then - n_words = get_arraysize_double(node_filt, "bins") - else - n_words = get_arraysize_integer(node_filt, "bins") + select case(temp_str) + case ("energy", "energyout", "mu", "polar", "azimuthal") + if (.not. check_for_node(node_filt, "bins")) then + call fatal_error("Bins not set in filter on tally " & + // trim(to_str(t % id))) end if - else - call fatal_error("Bins not set in filter on tally " & - // trim(to_str(t % id))) - end if + n_words = get_arraysize_double(node_filt, "bins") + case ("mesh", "universe", "material", "cell", "distribcell", & + "cellborn", "surface", "delayedgroup") + if (.not. check_for_node(node_filt, "bins")) then + call fatal_error("Bins not set in filter on tally " & + // trim(to_str(t % id))) + end if + n_words = get_arraysize_integer(node_filt, "bins") + end select ! Determine type of filter select case (temp_str) @@ -3282,6 +3285,39 @@ contains ! Set the filter index in the tally find_filter array t % find_filter(FILTER_AZIMUTHAL) = j + case ('linlinenergy') + ! Allocate and declare the filter type. + allocate(LinLinEnergyFilter::t % filters(j) % obj) + select type (filt => t % filters(j) % obj) + type is (LinLinEnergyFilter) + filt % n_bins = 1 + ! Make sure this is continuous-energy mode. + if (.not. run_CE) then + call fatal_error("LinLinEnergy filters are only supported for & + &continuous-energy transport calculations") + end if + + ! Allocate and store energy grid. + if (.not. check_for_node(node_filt, "energy")) then + call fatal_error("Energy grid not specified for LinLinEnergy & + &filter on tally " // trim(to_str(t % id))) + end if + n_words = get_arraysize_double(node_filt, "energy") + allocate(filt % energy(n_words)) + call get_node_array(node_filt, "energy", filt % energy) + + ! Allocate and store interpolant values. + if (.not. check_for_node(node_filt, "y")) then + call fatal_error("y values not specified for LinLinEnergy & + &filter on tally " // trim(to_str(t % id))) + end if + n_words = get_arraysize_double(node_filt, "y") + allocate(filt % y(n_words)) + call get_node_array(node_filt, "y", filt % y) + end select + ! Set the filter index in the tally find_filter array + t % find_filter(FILTER_LINLINENERGY) = j + case default ! Specified tally filter is invalid, raise error call fatal_error("Unknown filter type '" & diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index 37265cc5f..064e9d646 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -191,6 +191,19 @@ module tally_filter procedure :: text_label => text_label_azimuthal end type AzimuthalFilter +!=============================================================================== +! LinLinEnergyFilter +!=============================================================================== + type, extends(TallyFilter) :: LinLinEnergyFilter + real(8), allocatable :: energy(:) + real(8), allocatable :: y(:) + + contains + procedure :: get_next_bin => get_next_bin_linlinenergy + procedure :: to_statepoint => to_statepoint_linlinenergy + procedure :: text_label => text_label_linlinenergy + end type LinLinEnergyFilter + contains !=============================================================================== @@ -1209,6 +1222,83 @@ contains // trim(to_str(E1)) // ")" end function text_label_azimuthal +!=============================================================================== +! LinLinEnergyFilter methods +!=============================================================================== + subroutine get_next_bin_linlinenergy(this, p, estimator, current_bin, & + next_bin, weight) + class(LinLinEnergyFilter), intent(in) :: this + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + integer, value, intent(in) :: current_bin + integer, intent(out) :: next_bin + real(8), intent(out) :: weight + + integer :: n, indx + real(8) :: E, f + + select type(this) + type is (LinLinEnergyFilter) + if (current_bin == NO_BIN_FOUND) then + n = size(this % energy) + + ! Make sure the correct energy is used. + if (estimator == ESTIMATOR_TRACKLENGTH) then + E = p % E + else + E = p % last_E + end if + + ! Check if energy of the particle is within energy bins. + if (E < this % energy(1) .or. E > this % energy(n)) then + next_bin = NO_BIN_FOUND + weight = ONE + + else + ! Search to find incoming energy bin. + indx = binary_search(this % energy, n, E) + + ! Compute an interpolation factor between nearest bins. + f = (E - this % energy(indx)) & + / (this % energy(indx+1) - this % energy(indx)) + + ! Interpolate on the lin-lin grid. + next_bin = 1 + weight = (ONE - f) * this % y(indx) + f * this % y(indx+1) + end if + + else + next_bin = NO_BIN_FOUND + weight = ONE + end if + end select + end subroutine get_next_bin_linlinenergy + + subroutine to_statepoint_linlinenergy(this, filter_group) + class(LinLinEnergyFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + + select type(this) + type is (LinLinEnergyFilter) + call write_dataset(filter_group, "type", "linlinenergy") + call write_dataset(filter_group, "energy", this % energy) + call write_dataset(filter_group, "y", this % y) + end select + end subroutine to_statepoint_linlinenergy + + function text_label_linlinenergy(this, bin) result(label) + class(LinLinEnergyFilter), intent(in) :: this + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label + + select type(this) + type is (LinLinEnergyFilter) + label = "Lin-Lin Energy Interpolation [" & + // trim(to_str(this % energy(1))) // ", ..., " & + // trim(to_str(this % energy(size(this % energy)))) // "]" + end select + end function text_label_linlinenergy + !=============================================================================== ! FIND_OFFSET (for distribcell) uses a given map number, a target cell ID, and ! a target offset to build a string which is the path from the base universe to