Filter weight implementation (#3345)

Co-authored-by: Grego01-biot <grego01@pc-neutronic-06.psfc.mit.edu>
Co-authored-by: Patrick Shriwise <pshriwise@gmail.com>
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

CMakeLists.txt

@@ -425,6 +425,7 @@ list(APPEND libopenmc_SOURCES
   src/tallies/filter_surface.cpp
   src/tallies/filter_time.cpp
   src/tallies/filter_universe.cpp
+  src/tallies/filter_weight.cpp
   src/tallies/filter_zernike.cpp
   src/tallies/tally.cpp
   src/tallies/tally_scoring.cpp

docs/source/pythonapi/base.rst

@@ -143,6 +143,7 @@ Constructing Tallies
    openmc.SpatialLegendreFilter
    openmc.SphericalHarmonicsFilter
    openmc.TimeFilter
+   openmc.WeightFilter
    openmc.ZernikeFilter
    openmc.ZernikeRadialFilter
    openmc.ParentNuclideFilter

docs/source/pythonapi/capi.rst

@@ -91,6 +91,7 @@ Classes
    Tally
    UniverseFilter
    UnstructuredMesh
+   WeightFilter
    WeightWindows
    ZernikeFilter
    ZernikeRadialFilter

include/openmc/tallies/filter.h

@@ -44,6 +44,7 @@ enum class FilterType {
   SURFACE,
   TIME,
   UNIVERSE,
+  WEIGHT,
   ZERNIKE,
   ZERNIKE_RADIAL
 };

include/openmc/tallies/filter_weight.h

@@ -0,0 +1,51 @@
+#ifndef OPENMC_TALLIES_FILTER_WEIGHT_H
+#define OPENMC_TALLIES_FILTER_WEIGHT_H
+
+#include <string>
+
+#include "openmc/span.h"
+#include "openmc/tallies/filter.h"
+#include "openmc/vector.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Bins the weights of the particles.
+//==============================================================================
+
+class WeightFilter : public Filter {
+public:
+  //----------------------------------------------------------------------------
+  // Constructors, destructors
+
+  ~WeightFilter() = default;
+
+  //----------------------------------------------------------------------------
+  // Methods
+
+  std::string type_str() const override { return "weight"; }
+  FilterType type() const override { return FilterType::WEIGHT; }
+
+  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 vector<double>& bins() const { return bins_; }
+  void set_bins(span<const double> bins);
+
+protected:
+  //----------------------------------------------------------------------------
+  // Data members
+  vector<double> bins_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_WEIGHT_H

openmc/filter.py

@@ -25,7 +25,7 @@ _FILTER_TYPES = (
     'energyout', 'mu', 'musurface', 'polar', 'azimuthal', 'distribcell', 'delayedgroup',
     'energyfunction', 'cellfrom', 'materialfrom', 'legendre', 'spatiallegendre',
     'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'cellinstance',
-    'collision', 'time', 'parentnuclide'
+    'collision', 'time', 'parentnuclide', 'weight'
 )
 
 _CURRENT_NAMES = (
@@ -2328,3 +2328,26 @@ class EnergyFunctionFilter(Filter):
             {self.short_name.lower(): filter_bins})])
 
         return df
+
+
+class WeightFilter(RealFilter):
+    """Bins tally events based on the incoming particle weight.
+
+    Parameters
+    ----------
+    Values : Iterable of float
+        A list or iterable of the weight boundaries, as float values.
+    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 weights by which to filter
+    num_bins : int
+        The number of filter bins
+    values : numpy.ndarray
+        Array of weight boundaries
+    """

openmc/lib/filter.py

@@ -610,6 +610,10 @@ class UniverseFilter(Filter):
     filter_type = 'universe'
 
 
+class WeightFilter(Filter):
+    filter_type = 'weight'
+
+
 class ZernikeFilter(Filter):
     filter_type = 'zernike'
 

src/hdf5_interface.cpp

@@ -92,7 +92,7 @@ void get_shape_attr(hid_t obj_id, const char* name, hsize_t* dims)
   H5Aclose(attr);
 }
 
-hid_t create_group(hid_t parent_id, char const* name)
+hid_t create_group(hid_t parent_id, const char* name)
 {
   hid_t out = H5Gcreate(parent_id, name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
   if (out < 0) {

src/tallies/filter.cpp

@@ -36,6 +36,7 @@
 #include "openmc/tallies/filter_surface.h"
 #include "openmc/tallies/filter_time.h"
 #include "openmc/tallies/filter_universe.h"
+#include "openmc/tallies/filter_weight.h"
 #include "openmc/tallies/filter_zernike.h"
 #include "openmc/xml_interface.h"
 
@@ -154,6 +155,8 @@ Filter* Filter::create(const std::string& type, int32_t id)
     return Filter::create<TimeFilter>(id);
   } else if (type == "universe") {
     return Filter::create<UniverseFilter>(id);
+  } else if (type == "weight") {
+    return Filter::create<WeightFilter>(id);
   } else if (type == "zernike") {
     return Filter::create<ZernikeFilter>(id);
   } else if (type == "zernikeradial") {

src/tallies/filter_weight.cpp

@@ -0,0 +1,63 @@
+#include "openmc/tallies/filter_weight.h"
+
+#include <algorithm> // for is_sorted
+#include <stdexcept> // for runtime_error
+
+#include <fmt/core.h>
+
+#include "openmc/search.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+//==============================================================================
+// WeightFilter implementation
+//==============================================================================
+
+void WeightFilter::from_xml(pugi::xml_node node)
+{
+  auto bins = get_node_array<double>(node, "bins");
+  this->set_bins(bins);
+}
+
+void WeightFilter::set_bins(span<const double> bins)
+{
+  if (!std::is_sorted(bins.begin(), bins.end())) {
+    throw std::runtime_error {"Weight bins must be monotonically increasing."};
+  }
+
+  // Clear existing bins
+  bins_.clear();
+  bins_.reserve(bins.size());
+
+  // Copy bins
+  bins_.insert(bins_.end(), bins.begin(), bins.end());
+  n_bins_ = bins_.size() - 1;
+}
+
+void WeightFilter::get_all_bins(
+  const Particle& p, TallyEstimator estimator, FilterMatch& match) const
+{
+  // Get particle weight
+  double wgt = p.wgt_last();
+
+  // Bin the weight
+  if (wgt >= bins_.front() && wgt <= bins_.back()) {
+    auto bin = lower_bound_index(bins_.begin(), bins_.end(), wgt);
+    match.bins_.push_back(bin);
+    match.weights_.push_back(1.0);
+  }
+}
+
+void WeightFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  write_dataset(filter_group, "bins", bins_);
+}
+
+std::string WeightFilter::text_label(int bin) const
+{
+  return fmt::format("Weight [{}, {}]", bins_[bin], bins_[bin + 1]);
+}
+
+} // namespace openmc

tests/unit_tests/test_filter_weight.py

@@ -0,0 +1,44 @@
+import openmc
+import numpy as np
+
+
+def test_weightfilter(run_in_tmpdir):
+    steel = openmc.Material(name='Stainless Steel')
+    steel.set_density('g/cm3', 8.00)
+    steel.add_nuclide('Fe56', 1.0)
+
+    sphere = openmc.Sphere(r=50.0, boundary_type='vacuum')
+    cell = openmc.Cell(region=-sphere, fill=steel)
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([cell])
+    model.settings.particles = 100
+    model.settings.batches = 10
+
+    model.settings.source = openmc.IndependentSource(
+        energy=openmc.stats.delta_function(14e6),
+    )
+    model.settings.run_mode = "fixed source"
+
+    radius = list(range(1, 50))
+    sphere_mesh = openmc.SphericalMesh(radius)
+    mesh_filter = openmc.MeshFilter(sphere_mesh)
+    weight_filter = openmc.WeightFilter(
+        [0.999, 0.9999, 0.99999, 0.999999, 1.0, 1.000001 ,1.00001, 1.0001, 1.001]
+    )
+
+    tally = openmc.Tally()
+    tally.filters = [mesh_filter, weight_filter]
+    tally.estimator = 'analog'
+    tally.scores = ['flux']
+    model.tallies = openmc.Tallies([tally])
+
+    # Run OpenMC
+    model.run(apply_tally_results=True)
+
+    # Get current binned by mu
+    neutron_flux = tally.mean.reshape(48, 8)
+
+    # All contributions should show up in the fourth bin
+    assert np.all(neutron_flux[:, 3] != 0.0)
+    neutron_flux[:, 3] = 0.0
+    assert np.all(neutron_flux == 0.0)

tests/unit_tests/test_filters.py

@@ -312,3 +312,24 @@ def test_energy_filter():
     msg = 'Unable to set "filter value" to "-1.2" since it is less than "0.0"'
     with raises(ValueError, match=msg):
         openmc.EnergyFilter([-1.2, 0.25, 0.5])
+
+
+def test_weight():
+    f = openmc.WeightFilter([0.01, 0.1, 1.0, 10.0])
+    expected_bins = [[0.01, 0.1], [0.1, 1.0], [1.0, 10.0]]
+
+    assert np.allclose(f.bins, expected_bins)
+    assert len(f.bins) == 3
+
+    # Make sure __repr__ works
+    repr(f)
+
+    # to_xml_element()
+    elem = f.to_xml_element()
+    assert elem.tag == 'filter'
+    assert elem.attrib['type'] == 'weight'
+
+    # from_xml_element()
+    new_f = openmc.Filter.from_xml_element(elem)
+    assert new_f.id == f.id
+    assert np.allclose(new_f.bins, f.bins)