Implement Legendre expansion filter for change in scattering angle

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

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

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 *

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

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 :: &

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')

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

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