Implement filter for cosine of angle of surface crossing (#2768)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

CMakeLists.txt

@@ -416,6 +416,7 @@ list(APPEND libopenmc_SOURCES
   src/tallies/filter_meshborn.cpp
   src/tallies/filter_meshsurface.cpp
   src/tallies/filter_mu.cpp
+  src/tallies/filter_musurface.cpp
   src/tallies/filter_particle.cpp
   src/tallies/filter_polar.cpp
   src/tallies/filter_sph_harm.cpp

docs/source/pythonapi/base.rst

@@ -133,6 +133,7 @@ Constructing Tallies
    openmc.EnergyFilter
    openmc.EnergyoutFilter
    openmc.MuFilter
+   openmc.MuSurfaceFilter
    openmc.PolarFilter
    openmc.AzimuthalFilter
    openmc.DistribcellFilter

include/openmc/tallies/filter.h

@@ -36,6 +36,7 @@ enum class FilterType {
   MESHBORN,
   MESH_SURFACE,
   MU,
+  MUSURFACE,
   PARTICLE,
   POLAR,
   SPHERICAL_HARMONICS,

include/openmc/tallies/filter_mu.h

@@ -40,7 +40,7 @@ public:
 
   void set_bins(gsl::span<double> bins);
 
-private:
+protected:
   //----------------------------------------------------------------------------
   // Data members
 

include/openmc/tallies/filter_musurface.h

@@ -0,0 +1,34 @@
+#ifndef OPENMC_TALLIES_FILTER_MU_SURFACE_H
+#define OPENMC_TALLIES_FILTER_MU_SURFACE_H
+
+#include <gsl/gsl-lite.hpp>
+
+#include "openmc/tallies/filter_mu.h"
+#include "openmc/vector.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Bins the incoming-outgoing direction cosine.  This is only used for surface
+//! crossings.
+//==============================================================================
+
+class MuSurfaceFilter : public MuFilter {
+public:
+  //----------------------------------------------------------------------------
+  // Constructors, destructors
+
+  ~MuSurfaceFilter() = default;
+
+  //----------------------------------------------------------------------------
+  // Methods
+
+  std::string type_str() const override { return "musurface"; }
+  FilterType type() const override { return FilterType::MUSURFACE; }
+
+  void get_all_bins(const Particle& p, TallyEstimator estimator,
+    FilterMatch& match) const override;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_MU_SURFACE_H

openmc/filter.py

@@ -22,7 +22,7 @@ from ._xml import get_text
 
 _FILTER_TYPES = (
     'universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy',
-    'energyout', 'mu', 'polar', 'azimuthal', 'distribcell', 'delayedgroup',
+    'energyout', 'mu', 'musurface', 'polar', 'azimuthal', 'distribcell', 'delayedgroup',
     'energyfunction', 'cellfrom', 'materialfrom', 'legendre', 'spatiallegendre',
     'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'cellinstance',
     'collision', 'time'
@@ -765,7 +765,7 @@ class ParticleFilter(Filter):
     @Filter.bins.setter
     def bins(self, bins):
         cv.check_type('bins', bins, Sequence, str)
-        bins = np.atleast_1d(bins)    
+        bins = np.atleast_1d(bins)
         for edge in bins:
             cv.check_value('filter bin', edge, _PARTICLES)
         self._bins = bins
@@ -1850,6 +1850,44 @@ class MuFilter(RealFilter):
                 cv.check_less_than('filter value', x, 1., equality=True)
 
 
+class MuSurfaceFilter(MuFilter):
+    """Bins tally events based on the angle of surface crossing.
+
+    This filter bins events based on the cosine of the angle between the
+    direction of the particle and the normal to the surface at the point it
+    crosses. Only used in conjunction with a SurfaceFilter and current score.
+
+    .. versionadded:: 0.15.1
+
+    Parameters
+    ----------
+    values : int or Iterable of Real
+        A grid of surface crossing angles which the events will be divided into.
+        Values represent the cosine of the angle between the direction of the
+        particle and the normal to the surface at the point it crosses. If an
+        iterable is given, the values will be used explicitly as grid points. If
+        a single int is given, the range [-1, 1] will be divided equally into
+        that number of bins.
+    filter_id : int
+        Unique identifier for the filter
+
+    Attributes
+    ----------
+    values : numpy.ndarray
+        An array of values for which each successive pair constitutes a range of
+        surface crossing angle cosines for a single bin.
+    id : int
+        Unique identifier for the filter
+    bins : numpy.ndarray
+        An array of shape (N, 2) where each row is a pair of cosines of surface
+        crossing angle for a single filter
+    num_bins : Integral
+        The number of filter bins
+
+    """
+    # Note: inherits implementation from MuFilter
+
+
 class PolarFilter(RealFilter):
     """Bins tally events based on the incident particle's direction.
 

openmc/lib/filter.py

@@ -21,9 +21,9 @@ __all__ = [
     'CellInstanceFilter', 'CollisionFilter', 'DistribcellFilter', 'DelayedGroupFilter',
     'EnergyFilter', 'EnergyoutFilter', 'EnergyFunctionFilter', 'LegendreFilter',
     'MaterialFilter', 'MaterialFromFilter', 'MeshFilter', 'MeshBornFilter',
-    'MeshSurfaceFilter', 'MuFilter', 'ParticleFilter',
-    'PolarFilter', 'SphericalHarmonicsFilter', 'SpatialLegendreFilter', 'SurfaceFilter',
-    'UniverseFilter', 'ZernikeFilter', 'ZernikeRadialFilter', 'filters'
+    'MeshSurfaceFilter', 'MuFilter', 'MuSurfaceFilter', 'ParticleFilter', 'PolarFilter',
+    'SphericalHarmonicsFilter', 'SpatialLegendreFilter', 'SurfaceFilter', 'UniverseFilter',
+    'ZernikeFilter', 'ZernikeRadialFilter', 'filters'
 ]
 
 # Tally functions
@@ -540,6 +540,10 @@ class MuFilter(Filter):
     filter_type = 'mu'
 
 
+class MuSurfaceFilter(Filter):
+    filter_type = 'musurface'
+
+
 class ParticleFilter(Filter):
     filter_type = 'particle'
 
@@ -642,6 +646,7 @@ _FILTER_TYPE_MAP = {
     'meshborn': MeshBornFilter,
     'meshsurface': MeshSurfaceFilter,
     'mu': MuFilter,
+    'musurface': MuSurfaceFilter,
     'particle': ParticleFilter,
     'polar': PolarFilter,
     'sphericalharmonics': SphericalHarmonicsFilter,

src/tallies/filter.cpp

@@ -26,6 +26,7 @@
 #include "openmc/tallies/filter_meshborn.h"
 #include "openmc/tallies/filter_meshsurface.h"
 #include "openmc/tallies/filter_mu.h"
+#include "openmc/tallies/filter_musurface.h"
 #include "openmc/tallies/filter_particle.h"
 #include "openmc/tallies/filter_polar.h"
 #include "openmc/tallies/filter_sph_harm.h"
@@ -133,6 +134,8 @@ Filter* Filter::create(const std::string& type, int32_t id)
     return Filter::create<MeshSurfaceFilter>(id);
   } else if (type == "mu") {
     return Filter::create<MuFilter>(id);
+  } else if (type == "musurface") {
+    return Filter::create<MuSurfaceFilter>(id);
   } else if (type == "particle") {
     return Filter::create<ParticleFilter>(id);
   } else if (type == "polar") {

src/tallies/filter_musurface.cpp

@@ -0,0 +1,36 @@
+#include "openmc/tallies/filter_musurface.h"
+
+#include <cmath> // for abs, copysign
+
+#include "openmc/search.h"
+#include "openmc/surface.h"
+#include "openmc/tallies/tally_scoring.h"
+
+namespace openmc {
+
+void MuSurfaceFilter::get_all_bins(
+  const Particle& p, TallyEstimator estimator, FilterMatch& match) const
+{
+  // Get surface normal (and make sure it is a unit vector)
+  const auto surf {model::surfaces[std::abs(p.surface()) - 1].get()};
+  auto n = surf->normal(p.r());
+  n /= n.norm();
+
+  // Determine whether normal should be pointing in or out
+  if (p.surface() < 0)
+    n *= -1;
+
+  // Determine cosine of angle between normal and particle direction
+  double mu = p.u().dot(n);
+  if (std::abs(mu) > 1.0)
+    mu = std::copysign(1.0, mu);
+
+  // Find matching bin
+  if (mu >= bins_.front() && mu <= bins_.back()) {
+    auto bin = lower_bound_index(bins_.begin(), bins_.end(), mu);
+    match.bins_.push_back(bin);
+    match.weights_.push_back(1.0);
+  }
+}
+
+} // namespace openmc

tests/regression_tests/filter_musurface/inputs_true.dat

@@ -0,0 +1,37 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material depletable="true" id="1">
+      <density units="g/cm3" value="10.0"/>
+      <nuclide ao="1.0" name="U235"/>
+    </material>
+    <material id="2">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="1.0" name="Zr90"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1" universe="1"/>
+    <cell id="2" material="2" region="1 -2" universe="1"/>
+    <surface coeffs="0.0 0.0 1.0" id="1" type="z-cylinder"/>
+    <surface boundary="vacuum" coeffs="0.0 0.0 3.0" id="2" type="z-cylinder"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>1000</particles>
+    <batches>5</batches>
+    <inactive>0</inactive>
+  </settings>
+  <tallies>
+    <filter id="1" type="surface">
+      <bins>1</bins>
+    </filter>
+    <filter id="2" type="musurface">
+      <bins>-1.0 -0.5 0.0 0.5 1.0</bins>
+    </filter>
+    <tally id="1">
+      <filters>1 2</filters>
+      <scores>current</scores>
+    </tally>
+  </tallies>
+</model>

tests/regression_tests/filter_musurface/results_true.dat

@@ -0,0 +1,11 @@
+k-combined:
+1.157005E-01 7.587090E-03
+tally 1:
+0.000000E+00
+0.000000E+00
+0.000000E+00
+0.000000E+00
+8.770000E-01
+1.608710E-01
+3.909000E+00
+3.063035E+00

tests/regression_tests/filter_musurface/test.py

@@ -0,0 +1,43 @@
+import numpy as np
+from math import pi
+
+import openmc
+import pytest
+
+from tests.testing_harness import PyAPITestHarness
+
+
+@pytest.fixture
+def model():
+    model = openmc.Model()
+    fuel = openmc.Material()
+    fuel.set_density('g/cm3', 10.0)
+    fuel.add_nuclide('U235', 1.0)
+    zr = openmc.Material()
+    zr.set_density('g/cm3', 1.0)
+    zr.add_nuclide('Zr90', 1.0)
+
+    cyl1 = openmc.ZCylinder(r=1.0)
+    cyl2 = openmc.ZCylinder(r=3.0, boundary_type='vacuum')
+    cell1 = openmc.Cell(fill=fuel, region=-cyl1)
+    cell2 = openmc.Cell(fill=zr, region=+cyl1 & -cyl2)
+    model.geometry = openmc.Geometry([cell1, cell2])
+
+    model.settings.batches = 5
+    model.settings.inactive = 0
+    model.settings.particles = 1000
+
+    # Create a tally for current through the first surface binned by mu
+    surf_filter = openmc.SurfaceFilter([cyl1])
+    mu_filter = openmc.MuSurfaceFilter([-1.0, -0.5, 0.0, 0.5, 1.0])
+    tally = openmc.Tally()
+    tally.filters = [surf_filter, mu_filter]
+    tally.scores = ['current']
+    model.tallies.append(tally)
+
+    return model
+
+
+def test_filter_musurface(model):
+    harness = PyAPITestHarness('statepoint.5.h5', model)
+    harness.main()

tests/unit_tests/test_filter_musurface.py

@@ -0,0 +1,38 @@
+import openmc
+
+
+def test_musurface(run_in_tmpdir):
+    sphere = openmc.Sphere(r=1.0, boundary_type='vacuum')
+    cell = openmc.Cell(region=-sphere)
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([cell])
+    model.settings.particles = 1000
+    model.settings.batches = 10
+    E = 1.0
+    model.settings.source = openmc.IndependentSource(
+        space=openmc.stats.Point(),
+        angle=openmc.stats.Isotropic(),
+        energy=openmc.stats.delta_function(E),
+    )
+    model.settings.run_mode = "fixed source"
+
+    filter1 = openmc.MuSurfaceFilter(200)
+    filter2 = openmc.SurfaceFilter(sphere)
+    tally = openmc.Tally()
+    tally.filters = [filter1, filter2]
+    tally.scores = ['current']
+    model.tallies = openmc.Tallies([tally])
+
+    # Run OpenMC
+    sp_filename = model.run()
+
+    # Get current binned by mu
+    with openmc.StatePoint(sp_filename) as sp:
+        current_mu = sp.tallies[tally.id].mean.ravel()
+
+    # All contributions should show up in last bin
+    assert current_mu[-1] == 1.0
+    for element in current_mu[:-1]:
+        assert element == 0.0
+
+