Merge pull request #1798 from ameliajo/collisionFilter

Added Collision filter

CMakeLists.txt

@@ -330,6 +330,7 @@ list(APPEND libopenmc_SOURCES
   src/tallies/filter_distribcell.cpp
   src/tallies/filter_energyfunc.cpp
   src/tallies/filter_energy.cpp
+  src/tallies/filter_collision.cpp
   src/tallies/filter_legendre.cpp
   src/tallies/filter_material.cpp
   src/tallies/filter_mesh.cpp

docs/source/pythonapi/base.rst

@@ -118,6 +118,7 @@ Constructing Tallies
    openmc.CellFromFilter
    openmc.CellbornFilter
    openmc.CellInstanceFilter
+   openmc.CollisionFilter
    openmc.SurfaceFilter
    openmc.MeshFilter
    openmc.MeshSurfaceFilter

include/openmc/tallies/filter_collision.h

@@ -0,0 +1,54 @@
+#ifndef OPENMC_TALLIES_FILTER_COLLISIONS_H
+#define OPENMC_TALLIES_FILTER_COLLISIONS_H
+
+#include <vector>
+#include <unordered_map>
+#include <gsl/gsl>
+
+#include "openmc/tallies/filter.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Bins the incident neutron energy.
+//==============================================================================
+
+class CollisionFilter : public Filter {
+public:
+  //----------------------------------------------------------------------------
+  // Constructors, destructors
+
+  ~CollisionFilter() = default;
+
+  //----------------------------------------------------------------------------
+  // Methods
+
+  std::string type() const override { return "collision"; }
+
+  void from_xml(pugi::xml_node node) override;
+
+  void get_all_bins(const Particle& p, TallyEstimator estimator, 
+    FilterMatch& match) const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  //----------------------------------------------------------------------------
+  // Accessors
+
+  const std::vector<int>& bins() const { return bins_; }
+  void set_bins(gsl::span<const int> bins);
+
+protected:
+  //----------------------------------------------------------------------------
+  // Data members
+
+  std::vector<int> bins_;
+
+  std::unordered_map<int,int> map_;
+
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_COLLISIONS_H

openmc/filter.py

@@ -22,7 +22,8 @@ _FILTER_TYPES = (
     'universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy',
     'energyout', 'mu', 'polar', 'azimuthal', 'distribcell', 'delayedgroup',
     'energyfunction', 'cellfrom', 'legendre', 'spatiallegendre',
-    'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'cellinstance'
+    'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'cellinstance',
+    'collision'
 )
 
 _CURRENT_NAMES = (
@@ -964,6 +965,64 @@ class MeshSurfaceFilter(MeshFilter):
         return pd.concat([df, pd.DataFrame(filter_dict)])
 
 
+class CollisionFilter(Filter):
+    """Bins tally events based on the number of collisions.
+
+    Parameters
+    ----------
+    bins : Iterable of int
+        A list or iterable of the number of collisions, as integer values. 
+        The events whose post-scattering collision number equals one of 
+        the provided values will be counted. 
+    filter_id : int
+        Unique identifier for the filter
+
+    Attributes
+    ----------
+    id : int
+        Unique identifier for the filter
+    bins : numpy.ndarray
+        An array of integer values representing the number of collisions events
+        by which to filter
+    num_bins : int
+        The number of filter bins
+
+    """
+
+    def __init__(self, bins, filter_id=None):
+        self.bins = np.asarray(bins)
+        self.id = filter_id
+
+    def __repr__(self):
+        string = type(self).__name__ + '\n'
+        string += '{: <16}=\t{}\n'.format('\tValues', self.bins)
+        string += '{: <16}=\t{}\n'.format('\tID', self.id)
+        return string
+
+    @Filter.bins.setter
+    def bins(self, bins):
+        Filter.bins.__set__(self, np.asarray(bins))
+
+    def check_bins(self, bins):
+        for x in bins:
+            # Values should be integers
+            cv.check_type('filter value', x, Integral)
+            cv.check_greater_than('filter value', x, 0, equality=True)
+
+    def to_xml_element(self):
+        """Return XML Element representing the Filter.
+
+        Returns
+        -------
+        element : xml.etree.ElementTree.Element
+            XML element containing filter data
+
+        """
+        element = super().to_xml_element()
+        element[0].text = ' '.join(str(x) for x in self.bins)
+        return element
+
+
 class RealFilter(Filter):
     """Tally modifier that describes phase-space and other characteristics
 

openmc/lib/filter.py

@@ -16,10 +16,10 @@ from .mesh import _get_mesh
 
 __all__ = [
     'Filter', 'AzimuthalFilter', 'CellFilter', 'CellbornFilter', 'CellfromFilter',
-    'CellInstanceFilter', 'DistribcellFilter', 'DelayedGroupFilter', 'EnergyFilter',
-    'EnergyoutFilter', 'EnergyFunctionFilter', 'LegendreFilter', 'MaterialFilter',
-    'MeshFilter', 'MeshSurfaceFilter', 'MuFilter', 'ParticleFilter', 'PolarFilter',
-    'SphericalHarmonicsFilter', 'SpatialLegendreFilter', 'SurfaceFilter',
+    'CellInstanceFilter', 'CollisionFilter', 'DistribcellFilter', 'DelayedGroupFilter', 
+    'EnergyFilter', 'EnergyoutFilter', 'EnergyFunctionFilter', 'LegendreFilter', 
+    'MaterialFilter', 'MeshFilter', 'MeshSurfaceFilter', 'MuFilter', 'ParticleFilter', 
+    'PolarFilter', 'SphericalHarmonicsFilter', 'SpatialLegendreFilter', 'SurfaceFilter',
     'UniverseFilter', 'ZernikeFilter', 'ZernikeRadialFilter', 'filters'
 ]
 
@@ -176,6 +176,10 @@ class EnergyFilter(Filter):
             self._index, len(energies), energies_p)
 
 
+class CollisionFilter(Filter):
+    filter_type = 'collision'
+
+
 class EnergyoutFilter(EnergyFilter):
     filter_type = 'energyout'
 

src/tallies/filter.cpp

@@ -15,6 +15,7 @@
 #include "openmc/tallies/filter_cellborn.h"
 #include "openmc/tallies/filter_cellfrom.h"
 #include "openmc/tallies/filter_cell_instance.h"
+#include "openmc/tallies/filter_collision.h"
 #include "openmc/tallies/filter_delayedgroup.h"
 #include "openmc/tallies/filter_distribcell.h"
 #include "openmc/tallies/filter_energyfunc.h"
@@ -125,6 +126,8 @@ Filter* Filter::create(const std::string& type, int32_t id)
     return Filter::create<EnergyFunctionFilter>(id);
   } else if (type == "energy") {
     return Filter::create<EnergyFilter>(id);
+  } else if (type == "collision") {
+    return Filter::create<CollisionFilter>(id);
   } else if (type == "energyout") {
     return Filter::create<EnergyoutFilter>(id);
   } else if (type == "legendre") {

src/tallies/filter_collision.cpp

@@ -0,0 +1,63 @@
+#include "openmc/tallies/filter_collision.h"
+
+#include <fmt/core.h>
+
+#include "openmc/capi.h"
+#include "openmc/search.h"
+#include "openmc/settings.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+//==============================================================================
+// CollisionFilter implementation
+//==============================================================================
+
+void CollisionFilter::from_xml(pugi::xml_node node)
+{
+  auto bins = get_node_array<int>(node, "bins");
+  this->set_bins(bins);
+}
+
+void CollisionFilter::set_bins(gsl::span<const int> bins)
+{
+  // Clear existing bins
+  bins_.clear();
+  bins_.reserve(bins.size());
+  map_.clear();
+
+  // Copy bins
+  for (gsl::index i = 0; i < bins.size(); ++i) {
+    bins_.push_back(bins[i]);
+    map_[bins[i]] = i;
+  }
+
+  n_bins_ = bins_.size();
+}
+
+void CollisionFilter::get_all_bins(
+  const Particle& p, TallyEstimator estimator, FilterMatch& match) const
+{
+  // Get the number of collisions for the particle
+  auto n = p.n_collision_;
+
+  // Bin the collision number. Must fit exactly the desired collision number.
+  auto search = map_.find(n);
+  if (search != map_.end()){
+    match.bins_.push_back(search->second);
+    match.weights_.push_back(1.0);
+  }
+}
+
+void CollisionFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  write_dataset(filter_group, "bins", bins_);
+}
+
+std::string CollisionFilter::text_label(int bin) const
+{
+  return fmt::format("Collision Number {}", bins_[bin]);
+}
+
+} // namespace openmc

src/tallies/tally.cpp

@@ -20,6 +20,7 @@
 #include "openmc/tallies/filter_cellfrom.h"
 #include "openmc/tallies/filter_delayedgroup.h"
 #include "openmc/tallies/filter_energy.h"
+#include "openmc/tallies/filter_collision.h"
 #include "openmc/tallies/filter_legendre.h"
 #include "openmc/tallies/filter_mesh.h"
 #include "openmc/tallies/filter_meshsurface.h"

tests/regression_tests/tallies/inputs_true.dat

@@ -350,6 +350,9 @@
   <filter id="12" type="universe">
     <bins>1 2 3 4 6 8</bins>
   </filter>
+  <filter id="15" type="collision">
+    <bins>1 2 5 3 6</bins>
+  </filter>
   <filter id="13" type="cell">
     <bins>10 21 22 23 60</bins>
   </filter>
@@ -422,23 +425,31 @@
   <tally id="15">
     <filters>10</filters>
     <scores>scatter nu-scatter</scores>
+    <estimator>analog</estimator>
   </tally>
   <tally id="16">
     <filters>11</filters>
     <scores>scatter nu-scatter flux total</scores>
+    <estimator>analog</estimator>
   </tally>
   <tally id="17">
     <filters>11</filters>
     <scores>flux total</scores>
+    <estimator>collision</estimator>
   </tally>
   <tally id="18">
     <filters>11</filters>
     <scores>flux total</scores>
+    <estimator>tracklength</estimator>
   </tally>
   <tally id="19">
     <filters>12</filters>
     <scores>total</scores>
   </tally>
+  <tally id="34">
+    <filters>15</filters>
+    <scores>scatter</scores>
+  </tally>
   <tally id="20">
     <filters>13</filters>
     <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 fission-q-prompt fission-q-recoverable decay-rate</scores>

tests/regression_tests/tallies/results_true.dat

@@ -1 +1 @@
-0f036e2c34f3e9cf350d4f6bac6b0c5933d77d0665a5d461a3911c29d8f1d49b74100e680920da286c63538fabb86cbcd17532c6aad6048ee78e6cef75dd76ae
+bde3149c47242d232bf0c138609d081f9cdb11851aac7486b485c764f0d954c79bad35ffd3d78f7c3485dbdc5e93120a14ff2a21da2856a5bad6bf598551d1ee

tests/regression_tests/tallies/test.py

@@ -103,23 +103,23 @@ def test_tallies():
     legendre_tally = Tally()
     legendre_tally.filters = [legendre_filter]
     legendre_tally.scores = ['scatter', 'nu-scatter']
-    legendre_tally.estimatir = 'analog'
+    legendre_tally.estimator = 'analog'
 
     harmonics_filter = SphericalHarmonicsFilter(order=4)
     harmonics_tally = Tally()
     harmonics_tally.filters = [harmonics_filter]
     harmonics_tally.scores = ['scatter', 'nu-scatter', 'flux', 'total']
-    harmonics_tally.estimatir = 'analog'
+    harmonics_tally.estimator = 'analog'
 
     harmonics_tally2 = Tally()
     harmonics_tally2.filters = [harmonics_filter]
     harmonics_tally2.scores = ['flux', 'total']
-    harmonics_tally2.estimatir = 'collision'
+    harmonics_tally2.estimator = 'collision'
 
     harmonics_tally3 = Tally()
     harmonics_tally3.filters = [harmonics_filter]
     harmonics_tally3.scores = ['flux', 'total']
-    harmonics_tally3.estimatir = 'tracklength'
+    harmonics_tally3.estimator = 'tracklength'
 
     universe_tally = Tally()
     universe_tally.filters = [
@@ -169,12 +169,19 @@ def test_tallies():
     fusion_tally.scores = ['H1-production', 'H2-production', 'H3-production',
         'He3-production', 'He4-production', 'heating', 'damage-energy']
 
+    n_collision = (1, 2, 5, 3, 6)
+    collision_filter = CollisionFilter(n_collision)
+    collision_tally = Tally()
+    collision_tally.filters = [collision_filter]
+    collision_tally.scores = ['scatter']
+
     model.tallies += [
         azimuthal_tally1, azimuthal_tally2, azimuthal_tally3,
         cellborn_tally, dg_tally, energy_tally, energyout_tally,
         transfer_tally, material_tally, mu_tally1, mu_tally2,
         polar_tally1, polar_tally2, polar_tally3, legendre_tally,
-        harmonics_tally, harmonics_tally2, harmonics_tally3, universe_tally]
+        harmonics_tally, harmonics_tally2, harmonics_tally3, 
+        universe_tally, collision_tally]
     model.tallies += score_tallies
     model.tallies += flux_tallies
     model.tallies += all_nuclide_tallies

tests/unit_tests/test_filters.py

@@ -46,6 +46,22 @@ def test_cell_instance():
     assert instances.apply(lambda x: x in (0, 1, 2)).all()
 
 
+def test_collision():
+    f = openmc.CollisionFilter([1, 5, 3, 2, 8])
+    assert f.bins[0] == 1
+    assert f.bins[1] == 5
+    assert f.bins[-1] == 8
+    assert len(f.bins) == 5
+
+    # Make sure __repr__ works
+    repr(f)
+
+    # to_xml_element()
+    elem = f.to_xml_element()
+    assert elem.tag == 'filter'
+    assert elem.attrib['type'] == 'collision'
+
+
 def test_legendre():
     n = 5
     f = openmc.LegendreFilter(n)