diff --git a/CMakeLists.txt b/CMakeLists.txt index 9517fe2ec..bb3575c9d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -420,6 +420,7 @@ set(LIBOPENMC_FORTRAN_SRC src/tallies/tally_filter_distribcell.F90 src/tallies/tally_filter_energy.F90 src/tallies/tally_filter_energyfunc.F90 + src/tallies/tally_filter_legendre.F90 src/tallies/tally_filter_material.F90 src/tallies/tally_filter_mesh.F90 src/tallies/tally_filter_meshsurface.F90 diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index 070b5bd20..84d2a39da 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -118,6 +118,7 @@ Constructing Tallies openmc.DistribcellFilter openmc.DelayedGroupFilter openmc.EnergyFunctionFilter + openmc.LegendreFilter openmc.Mesh openmc.Trigger openmc.TallyDerivative diff --git a/openmc/__init__.py b/openmc/__init__.py index 9b13fa692..243796a13 100644 --- a/openmc/__init__.py +++ b/openmc/__init__.py @@ -15,6 +15,7 @@ from openmc.surface import * from openmc.universe import * from openmc.lattice import * from openmc.filter import * +from openmc.filter_legendre import * from openmc.trigger import * from openmc.tally_derivative import * from openmc.tallies import * diff --git a/openmc/filter_legendre.py b/openmc/filter_legendre.py new file mode 100644 index 000000000..075fa02db --- /dev/null +++ b/openmc/filter_legendre.py @@ -0,0 +1,131 @@ +from numbers import Integral +from xml.etree import ElementTree as ET + +import numpy as np +import pandas as pd + +import openmc.checkvalue as cv +from . import Filter + + +class LegendreFilter(Filter): + r"""Score Legendre expansion moments up to specified order. + + This filter allows scores to be multiplied by Legendre polynomials of the + change in particle angle ($\mu$) up to a user-specified order. + + Parameters + ---------- + order : int + Maximum Legendre polynomial order + filter_id : int or None + Unique identifier for the filter + + Attributes + ---------- + order : int + Maximum Legendre polynomial order + id : int + Unique identifier for the filter + num_bins : int + The number of filter bins + + """ + + def __init__(self, order, filter_id=None): + self.order = order + self.id = filter_id + + def __hash__(self): + string = type(self).__name__ + '\n' + string += '{: <16}=\t{}\n'.format('\tOrder', self.order) + return hash(string) + + def __repr__(self): + string = type(self).__name__ + '\n' + string += '{: <16}=\t{}\n'.format('\tOrder', self.order) + string += '{: <16}=\t{}\n'.format('\tID', self.id) + return string + + @property + def order(self): + return self._order + + @order.setter + def order(self, order): + cv.check_type('Legendre order', order, Integral) + cv.check_greater_than('Legendre order', order, 0, equality=True) + self._order = order + + @property + def num_bins(self): + return self._order + 1 + + @classmethod + def from_hdf5(cls, group, **kwargs): + if group['type'].value.decode() != cls.short_name.lower(): + raise ValueError("Expected HDF5 data for filter type '" + + cls.short_name.lower() + "' but got '" + + group['type'].value.decode() + " instead") + + filter_id = int(group.name.split('/')[-1].lstrip('filter ')) + + out = cls(group['order'].value, filter_id) + + return out + + def get_pandas_dataframe(self, data_size, stride, **kwargs): + """Builds a Pandas DataFrame for the Filter's bins. + + This method constructs a Pandas DataFrame object for the filter with + columns annotated by filter bin information. This is a helper method for + :meth:`Tally.get_pandas_dataframe`. + + Parameters + ---------- + data_size : Integral + The total number of bins in the tally corresponding to this filter + stride : int + Stride in memory for the filter + + Returns + ------- + pandas.DataFrame + A Pandas DataFrame with a column that is filled with strings + indicating Legendre orders. The number of rows in the DataFrame is + the same as the total number of bins in the corresponding tally. + + See also + -------- + Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe() + + """ + # Initialize Pandas DataFrame + df = pd.DataFrame() + + bins = np.array(['P{}'.format(i) for i in range(self.order + 1)]) + filter_bins = np.repeat(bins, stride) + tile_factor = data_size // len(filter_bins) + filter_bins = np.tile(filter_bins, tile_factor) + df = pd.concat([df, pd.DataFrame( + {self.short_name.lower(): filter_bins})]) + + return df + + def to_xml_element(self): + """Return XML Element representing the filter. + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing Legendre filter data + + """ + element = ET.Element('filter') + element.set('id', str(self.id)) + element.set('type', self.short_name.lower()) + + subelement = ET.SubElement(element, 'order') + subelement.text = str(self.order) + + return element diff --git a/src/constants.F90 b/src/constants.F90 index 0e138f0ee..29b68ae1d 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -358,7 +358,7 @@ module constants integer, parameter :: NO_BIN_FOUND = -1 ! Tally filter and map types - integer, parameter :: N_FILTER_TYPES = 16 + integer, parameter :: N_FILTER_TYPES = 17 integer, parameter :: & FILTER_UNIVERSE = 1, & FILTER_MATERIAL = 2, & @@ -375,7 +375,8 @@ module constants FILTER_DELAYEDGROUP = 13, & FILTER_ENERGYFUNCTION = 14, & FILTER_CELLFROM = 15, & - FILTER_MESHSURFACE = 16 + FILTER_MESHSURFACE = 16, & + FILTER_LEGENDRE = 17 ! Mesh types integer, parameter :: & diff --git a/src/tallies/tally_filter.F90 b/src/tallies/tally_filter.F90 index 3d6eaaef2..056fc198a 100644 --- a/src/tallies/tally_filter.F90 +++ b/src/tallies/tally_filter.F90 @@ -17,6 +17,7 @@ module tally_filter use tally_filter_distribcell use tally_filter_energy use tally_filter_energyfunc + use tally_filter_legendre use tally_filter_material use tally_filter_mesh use tally_filter_meshsurface @@ -64,6 +65,8 @@ contains type_ = 'energyout' type is (EnergyFunctionFilter) type_ = 'energyfunction' + type is (LegendreFilter) + type_ = 'legendre' type is (MaterialFilter) type_ = 'material' type is (MeshFilter) @@ -135,6 +138,8 @@ contains allocate(EnergyoutFilter :: filters(index) % obj) case ('energyfunction') allocate(EnergyFunctionFilter :: filters(index) % obj) + case ('legendre') + allocate(LegendreFilter :: filters(index) % obj) case ('material') allocate(MaterialFilter :: filters(index) % obj) case ('mesh') diff --git a/src/tallies/tally_filter_legendre.F90 b/src/tallies/tally_filter_legendre.F90 new file mode 100644 index 000000000..59fb57a72 --- /dev/null +++ b/src/tallies/tally_filter_legendre.F90 @@ -0,0 +1,86 @@ +module tally_filter_legendre + + use, intrinsic :: ISO_C_BINDING + + use hdf5, only: HID_T + + use constants + use error + use hdf5_interface + use math, only: calc_pn + use particle_header, only: Particle + use string, only: to_str + use tally_filter_header + use xml_interface + + implicit none + private + +!=============================================================================== +! LEGENDREFILTER gives Legendre moments of the change in scattering angle +!=============================================================================== + + type, public, extends(TallyFilter) :: LegendreFilter + integer :: order + contains + procedure :: from_xml + procedure :: get_all_bins + procedure :: to_statepoint + procedure :: text_label + end type LegendreFilter + +contains + +!=============================================================================== +! LegendreFilter methods +!=============================================================================== + + subroutine from_xml(this, node) + class(LegendreFilter), intent(inout) :: this + type(XMLNode), intent(in) :: node + + ! Get specified order + call get_node_value(node, "order", this % order) + this % n_bins = this % order + 1 + end subroutine from_xml + + subroutine get_all_bins(this, p, estimator, match) + class(LegendreFilter), intent(in) :: this + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + type(TallyFilterMatch), intent(inout) :: match + + integer :: i + real(8) :: wgt + + ! TODO: Use recursive formula to calculate higher orders + do i = 0, this % order + wgt = calc_pn(i, p % mu) + call match % bins % push_back(i + 1) + call match % weights % push_back(wgt) + end do + end subroutine get_all_bins + + subroutine to_statepoint(this, filter_group) + class(LegendreFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + + call write_dataset(filter_group, "type", "legendre") + call write_dataset(filter_group, "n_bins", this % n_bins) + call write_dataset(filter_group, "order", this % order) + end subroutine to_statepoint + + function text_label(this, bin) result(label) + class(LegendreFilter), intent(in) :: this + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label + + label = "Legendre expansion order " // trim(to_str(bin - 1)) + end function text_label + +!=============================================================================== +! C API FUNCTIONS +!=============================================================================== + + +end module tally_filter_legendre diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index 0495ce14d..1f221c83b 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -354,6 +354,9 @@ contains j = FILTER_AZIMUTHAL type is (EnergyFunctionFilter) j = FILTER_ENERGYFUNCTION + type is (LegendreFilter) + j = FILTER_LEGENDRE + this % estimator = ESTIMATOR_ANALOG end select this % find_filter(j) = i end do