Meshborn filter (#2925)

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

CMakeLists.txt

@@ -411,6 +411,7 @@ list(APPEND libopenmc_SOURCES
   src/tallies/filter_material.cpp
   src/tallies/filter_materialfrom.cpp
   src/tallies/filter_mesh.cpp
+  src/tallies/filter_meshborn.cpp
   src/tallies/filter_meshsurface.cpp
   src/tallies/filter_mu.cpp
   src/tallies/filter_particle.cpp

docs/source/capi/index.rst

@@ -420,6 +420,16 @@ Functions
    :return: Return status (negative if an error occurred)
    :rtype: int
 
+.. c:function:: int openmc_mesh_filter_get_mesh(int32_t index, int32_t* index_mesh)
+
+   Get the mesh for a mesh filter
+
+   :param int32_t index: Index in the filters array
+   :param index_mesh: Index in the meshes array
+   :type index_mesh: int32_t*
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
 .. c:function:: int openmc_mesh_filter_set_mesh(int32_t index, int32_t index_mesh)
 
    Set the mesh for a mesh filter
@@ -429,6 +439,98 @@ Functions
    :return: Return status (negative if an error occurred)
    :rtype: int
 
+.. c:function:: int openmc_mesh_filter_get_translation(int32_t index, double translation[3])
+
+   Get the 3-D translation coordinates for a mesh filter
+
+   :param int32_t index: Index in the filters array
+   :param double[3] translation: 3-D translation coordinates
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_mesh_filter_set_translation(int32_t index, double translation[3])
+
+   Set the 3-D translation coordinates for a mesh filter
+
+   :param int32_t index: Index in the filters array
+   :param double[3] translation: 3-D translation coordinates
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_meshborn_filter_get_mesh(int32_t index, int32_t* index_mesh)
+
+   Get the mesh for a meshborn filter
+
+   :param int32_t index: Index in the filters array
+   :param index_mesh: Index in the meshes array
+   :type index_mesh: int32_t*
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_meshborn_filter_set_mesh(int32_t index, int32_t index_mesh)
+
+   Set the mesh for a meshborn filter
+
+   :param int32_t index: Index in the filters array
+   :param int32_t index_mesh: Index in the meshes array
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_meshborn_filter_get_translation(int32_t index, double translation[3])
+
+   Get the 3-D translation coordinates for a meshborn filter
+
+   :param int32_t index: Index in the filters array
+   :param double[3] translation: 3-D translation coordinates
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_meshborn_filter_set_translation(int32_t index, double translation[3])
+
+   Set the 3-D translation coordinates for a meshborn filter
+
+   :param int32_t index: Index in the filters array
+   :param double[3] translation: 3-D translation coordinates
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh)
+
+   Get the mesh for a mesh surface filter
+
+   :param int32_t index: Index in the filters array
+   :param index_mesh: Index in the meshes array
+   :type index_mesh: int32_t*
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh)
+
+   Set the mesh for a mesh surface filter
+
+   :param int32_t index: Index in the filters array
+   :param int32_t index_mesh: Index in the meshes array
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_meshsurface_filter_get_translation(int32_t index, double translation[3])
+
+   Get the 3-D translation coordinates for a mesh surface filter
+
+   :param int32_t index: Index in the filters array
+   :param double[3] translation: 3-D translation coordinates
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
+.. c:function:: int openmc_meshsurface_filter_set_translation(int32_t index, double translation[3])
+
+   Set the 3-D translation coordinates for a mesh surface filter
+
+   :param int32_t index: Index in the filters array
+   :param double[3] translation: 3-D translation coordinates
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
 .. c:function:: int openmc_next_batch()
 
    Simulate next batch of particles. Must be called after openmc_simulation_init().

docs/source/pythonapi/base.rst

@@ -128,6 +128,7 @@ Constructing Tallies
    openmc.CollisionFilter
    openmc.SurfaceFilter
    openmc.MeshFilter
+   openmc.MeshBornFilter
    openmc.MeshSurfaceFilter
    openmc.EnergyFilter
    openmc.EnergyoutFilter

docs/source/pythonapi/capi.rst

@@ -59,6 +59,7 @@ Classes
    MaterialFilter
    Material
    MeshFilter
+   MeshBornFilter
    MeshSurfaceFilter
    Nuclide
    RectilinearMesh

include/openmc/particle_data.h

@@ -262,6 +262,10 @@ public:
   int& cell_last(int i) { return cell_last_[i]; }
   const int& cell_last(int i) const { return cell_last_[i]; }
 
+  // Coordinates at birth
+  Position& r_born() { return r_born_; }
+  const Position& r_born() const { return r_born_; }
+
   // Coordinates of last collision or reflective/periodic surface
   // crossing for current tallies
   Position& r_last_current() { return r_last_current_; }
@@ -323,6 +327,7 @@ private:
   int n_coord_last_ {1};  //!< number of current coordinates
   vector<int> cell_last_; //!< coordinates for all levels
 
+  Position r_born_;         //!< coordinates at birth
   Position r_last_current_; //!< coordinates of the last collision or
                             //!< reflective/periodic surface crossing for
                             //!< current tallies

include/openmc/tallies/filter.h

@@ -33,6 +33,7 @@ enum class FilterType {
   MATERIAL,
   MATERIALFROM,
   MESH,
+  MESHBORN,
   MESH_SURFACE,
   MU,
   PARTICLE,

include/openmc/tallies/filter_meshborn.h

@@ -0,0 +1,26 @@
+#ifndef OPENMC_TALLIES_FILTER_MESHBORN_H
+#define OPENMC_TALLIES_FILTER_MESHBORN_H
+
+#include <cstdint>
+
+#include "openmc/position.h"
+#include "openmc/tallies/filter_mesh.h"
+
+namespace openmc {
+
+class MeshBornFilter : public MeshFilter {
+public:
+  //----------------------------------------------------------------------------
+  // Methods
+
+  std::string type_str() const override { return "meshborn"; }
+  FilterType type() const override { return FilterType::MESHBORN; }
+
+  void get_all_bins(const Particle& p, TallyEstimator estimator,
+    FilterMatch& match) const override;
+
+  std::string text_label(int bin) const override;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_MESHBORN_H

openmc/filter.py

@@ -982,6 +982,34 @@ class MeshFilter(Filter):
         if translation:
             out.translation = [float(x) for x in translation.split()]
         return out
+    
+
+class MeshBornFilter(MeshFilter):
+    """Filter events by the mesh cell a particle originated from.
+
+    Parameters
+    ----------
+    mesh : openmc.MeshBase
+        The mesh object that events will be tallied onto
+    filter_id : int
+        Unique identifier for the filter
+
+    Attributes
+    ----------
+    mesh : openmc.MeshBase
+        The mesh object that events will be tallied onto
+    id : int
+        Unique identifier for the filter
+    translation : Iterable of float
+        This array specifies a vector that is used to translate (shift)
+        the mesh for this filter
+    bins : list of tuple
+        A list of mesh indices for each filter bin, e.g. [(1, 1, 1), (2, 1, 1),
+        ...]
+    num_bins : Integral
+        The number of filter bins
+
+    """
 
 
 class MeshSurfaceFilter(MeshFilter):

openmc/lib/filter.py

@@ -20,7 +20,8 @@ __all__ = [
     'Filter', 'AzimuthalFilter', 'CellFilter', 'CellbornFilter', 'CellfromFilter',
     'CellInstanceFilter', 'CollisionFilter', 'DistribcellFilter', 'DelayedGroupFilter',
     'EnergyFilter', 'EnergyoutFilter', 'EnergyFunctionFilter', 'LegendreFilter',
-    'MaterialFilter', 'MaterialFromFilter', 'MeshFilter', 'MeshSurfaceFilter', 'MuFilter', 'ParticleFilter',
+    'MaterialFilter', 'MaterialFromFilter', 'MeshFilter', 'MeshBornFilter',
+    'MeshSurfaceFilter', 'MuFilter', 'ParticleFilter',
     'PolarFilter', 'SphericalHarmonicsFilter', 'SpatialLegendreFilter', 'SurfaceFilter',
     'UniverseFilter', 'ZernikeFilter', 'ZernikeRadialFilter', 'filters'
 ]
@@ -89,18 +90,36 @@ _dll.openmc_mesh_filter_get_mesh.errcheck = _error_handler
 _dll.openmc_mesh_filter_set_mesh.argtypes = [c_int32, c_int32]
 _dll.openmc_mesh_filter_set_mesh.restype = c_int
 _dll.openmc_mesh_filter_set_mesh.errcheck = _error_handler
-_dll.openmc_meshsurface_filter_get_mesh.argtypes = [c_int32, POINTER(c_int32)]
-_dll.openmc_meshsurface_filter_get_mesh.restype = c_int
-_dll.openmc_meshsurface_filter_get_mesh.errcheck = _error_handler
-_dll.openmc_meshsurface_filter_set_mesh.argtypes = [c_int32, c_int32]
-_dll.openmc_meshsurface_filter_set_mesh.restype = c_int
-_dll.openmc_meshsurface_filter_set_mesh.errcheck = _error_handler
 _dll.openmc_mesh_filter_get_translation.argtypes = [c_int32, POINTER(c_double*3)]
 _dll.openmc_mesh_filter_get_translation.restype = c_int
 _dll.openmc_mesh_filter_get_translation.errcheck = _error_handler
 _dll.openmc_mesh_filter_set_translation.argtypes = [c_int32, POINTER(c_double*3)]
 _dll.openmc_mesh_filter_set_translation.restype = c_int
 _dll.openmc_mesh_filter_set_translation.errcheck = _error_handler
+_dll.openmc_meshborn_filter_get_mesh.argtypes = [c_int32, POINTER(c_int32)]
+_dll.openmc_meshborn_filter_get_mesh.restype = c_int
+_dll.openmc_meshborn_filter_get_mesh.errcheck = _error_handler
+_dll.openmc_meshborn_filter_set_mesh.argtypes = [c_int32, c_int32]
+_dll.openmc_meshborn_filter_set_mesh.restype = c_int
+_dll.openmc_meshborn_filter_set_mesh.errcheck = _error_handler
+_dll.openmc_meshborn_filter_get_translation.argtypes = [c_int32, POINTER(c_double*3)]
+_dll.openmc_meshborn_filter_get_translation.restype = c_int
+_dll.openmc_meshborn_filter_get_translation.errcheck = _error_handler
+_dll.openmc_meshborn_filter_set_translation.argtypes = [c_int32, POINTER(c_double*3)]
+_dll.openmc_meshborn_filter_set_translation.restype = c_int
+_dll.openmc_meshborn_filter_set_translation.errcheck = _error_handler
+_dll.openmc_meshsurface_filter_get_mesh.argtypes = [c_int32, POINTER(c_int32)]
+_dll.openmc_meshsurface_filter_get_mesh.restype = c_int
+_dll.openmc_meshsurface_filter_get_mesh.errcheck = _error_handler
+_dll.openmc_meshsurface_filter_set_mesh.argtypes = [c_int32, c_int32]
+_dll.openmc_meshsurface_filter_set_mesh.restype = c_int
+_dll.openmc_meshsurface_filter_set_mesh.errcheck = _error_handler
+_dll.openmc_meshsurface_filter_get_translation.argtypes = [c_int32, POINTER(c_double*3)]
+_dll.openmc_meshsurface_filter_get_translation.restype = c_int
+_dll.openmc_meshsurface_filter_get_translation.errcheck = _error_handler
+_dll.openmc_meshsurface_filter_set_translation.argtypes = [c_int32, POINTER(c_double*3)]
+_dll.openmc_meshsurface_filter_set_translation.restype = c_int
+_dll.openmc_meshsurface_filter_set_translation.errcheck = _error_handler
 _dll.openmc_new_filter.argtypes = [c_char_p, POINTER(c_int32)]
 _dll.openmc_new_filter.restype = c_int
 _dll.openmc_new_filter.errcheck = _error_handler
@@ -347,6 +366,34 @@ class MaterialFromFilter(Filter):
 
 
 class MeshFilter(Filter):
+    """Mesh filter stored internally.
+
+    This class exposes a Mesh filter that is stored internally in the OpenMC
+    library. To obtain a view of a Mesh filter with a given ID, use the
+    :data:`openmc.lib.filters` mapping.
+
+    Parameters
+    ----------
+    mesh : openmc.lib.Mesh
+        Mesh to use for the filter
+    uid : int or None
+        Unique ID of the Mesh filter
+    new : bool
+        When `index` is None, this argument controls whether a new object is
+        created or a view of an existing object is returned.
+    index : int
+        Index in the `filters` array.
+
+    Attributes
+    ----------
+    filter_type : str
+        Type of filter
+    mesh : openmc.lib.Mesh
+        Mesh used for the filter
+    translation : Iterable of float
+        3-D coordinates of the translation vector
+
+    """
     filter_type = 'mesh'
 
     def __init__(self, mesh=None, uid=None, new=True, index=None):
@@ -375,7 +422,92 @@ class MeshFilter(Filter):
         _dll.openmc_mesh_filter_set_translation(self._index, (c_double*3)(*translation))
 
 
+class MeshBornFilter(Filter):
+    """MeshBorn filter stored internally.
+
+    This class exposes a MeshBorn filter that is stored internally in the OpenMC
+    library. To obtain a view of a MeshBorn filter with a given ID, use the
+    :data:`openmc.lib.filters` mapping.
+
+    Parameters
+    ----------
+    mesh : openmc.lib.Mesh
+        Mesh to use for the filter
+    uid : int or None
+        Unique ID of the MeshBorn filter
+    new : bool
+        When `index` is None, this argument controls whether a new object is
+        created or a view of an existing object is returned.
+    index : int
+        Index in the `filters` array.
+
+    Attributes
+    ----------
+    filter_type : str
+        Type of filter
+    mesh : openmc.lib.Mesh
+        Mesh used for the filter
+    translation : Iterable of float
+        3-D coordinates of the translation vector
+
+    """
+    filter_type = 'meshborn'
+
+    def __init__(self, mesh=None, uid=None, new=True, index=None):
+        super().__init__(uid, new, index)
+        if mesh is not None:
+            self.mesh = mesh
+
+    @property
+    def mesh(self):
+        index_mesh = c_int32()
+        _dll.openmc_meshborn_filter_get_mesh(self._index, index_mesh)
+        return _get_mesh(index_mesh.value)
+
+    @mesh.setter
+    def mesh(self, mesh):
+        _dll.openmc_meshborn_filter_set_mesh(self._index, mesh._index)
+
+    @property
+    def translation(self):
+        translation = (c_double*3)()
+        _dll.openmc_meshborn_filter_get_translation(self._index, translation)
+        return tuple(translation)
+
+    @translation.setter
+    def translation(self, translation):
+        _dll.openmc_meshborn_filter_set_translation(self._index, (c_double*3)(*translation))
+
+
 class MeshSurfaceFilter(Filter):
+    """MeshSurface filter stored internally.
+
+    This class exposes a MeshSurface filter that is stored internally in the
+    OpenMC library. To obtain a view of a MeshSurface filter with a given ID,
+    use the :data:`openmc.lib.filters` mapping.
+
+    Parameters
+    ----------
+    mesh : openmc.lib.Mesh
+        Mesh to use for the filter
+    uid : int or None
+        Unique ID of the MeshSurface filter
+    new : bool
+        When `index` is None, this argument controls whether a new object is
+        created or a view of an existing object is returned.
+    index : int
+        Index in the `filters` array.
+
+    Attributes
+    ----------
+    filter_type : str
+        Type of filter
+    mesh : openmc.lib.Mesh
+        Mesh used for the filter
+    translation : Iterable of float
+        3-D coordinates of the translation vector
+
+    """
     filter_type = 'meshsurface'
 
     def __init__(self, mesh=None, uid=None, new=True, index=None):
@@ -396,12 +528,12 @@ class MeshSurfaceFilter(Filter):
     @property
     def translation(self):
         translation = (c_double*3)()
-        _dll.openmc_mesh_filter_get_translation(self._index, translation)
+        _dll.openmc_meshsurface_filter_get_translation(self._index, translation)
         return tuple(translation)
 
     @translation.setter
     def translation(self, translation):
-        _dll.openmc_mesh_filter_set_translation(self._index, (c_double*3)(*translation))
+        _dll.openmc_meshsurface_filter_set_translation(self._index, (c_double*3)(*translation))
 
 
 class MuFilter(Filter):
@@ -507,6 +639,7 @@ _FILTER_TYPE_MAP = {
     'material': MaterialFilter,
     'materialfrom': MaterialFromFilter,
     'mesh': MeshFilter,
+    'meshborn': MeshBornFilter,
     'meshsurface': MeshSurfaceFilter,
     'mu': MuFilter,
     'particle': ParticleFilter,

src/particle.cpp

@@ -121,6 +121,7 @@ void Particle::from_source(const SourceSite* src)
   wgt_last() = src->wgt;
   r() = src->r;
   u() = src->u;
+  r_born() = src->r;
   r_last_current() = src->r;
   r_last() = src->r;
   u_last() = src->u;

src/tallies/filter.cpp

@@ -23,6 +23,7 @@
 #include "openmc/tallies/filter_material.h"
 #include "openmc/tallies/filter_materialfrom.h"
 #include "openmc/tallies/filter_mesh.h"
+#include "openmc/tallies/filter_meshborn.h"
 #include "openmc/tallies/filter_meshsurface.h"
 #include "openmc/tallies/filter_mu.h"
 #include "openmc/tallies/filter_particle.h"
@@ -126,6 +127,8 @@ Filter* Filter::create(const std::string& type, int32_t id)
     return Filter::create<MaterialFromFilter>(id);
   } else if (type == "mesh") {
     return Filter::create<MeshFilter>(id);
+  } else if (type == "meshborn") {
+    return Filter::create<MeshBornFilter>(id);
   } else if (type == "meshsurface") {
     return Filter::create<MeshSurfaceFilter>(id);
   } else if (type == "mu") {

src/tallies/filter_mesh.cpp

@@ -161,6 +161,7 @@ extern "C" int openmc_mesh_filter_get_translation(
   // Check the filter type
   const auto& filter = model::tally_filters[index];
   if (filter->type() != FilterType::MESH &&
+      filter->type() != FilterType::MESHBORN &&
       filter->type() != FilterType::MESH_SURFACE) {
     set_errmsg("Tried to get a translation from a non-mesh-based filter.");
     return OPENMC_E_INVALID_TYPE;
@@ -186,6 +187,7 @@ extern "C" int openmc_mesh_filter_set_translation(
   const auto& filter = model::tally_filters[index];
   // Check the filter type
   if (filter->type() != FilterType::MESH &&
+      filter->type() != FilterType::MESHBORN &&
       filter->type() != FilterType::MESH_SURFACE) {
     set_errmsg("Tried to set mesh on a non-mesh-based filter.");
     return OPENMC_E_INVALID_TYPE;

src/tallies/filter_meshborn.cpp

@@ -0,0 +1,61 @@
+#include "openmc/tallies/filter_meshborn.h"
+
+#include "openmc/capi.h"
+#include "openmc/constants.h"
+#include "openmc/error.h"
+#include "openmc/mesh.h"
+
+namespace openmc {
+
+void MeshBornFilter::get_all_bins(
+  const Particle& p, TallyEstimator estimator, FilterMatch& match) const
+{
+  Position r_born = p.r_born();
+
+  // apply translation if present
+  if (translated_) {
+    r_born -= translation();
+  }
+
+  auto bin = model::meshes[mesh_]->get_bin(r_born);
+  if (bin >= 0) {
+    match.bins_.push_back(bin);
+    match.weights_.push_back(1.0);
+  }
+}
+
+std::string MeshBornFilter::text_label(int bin) const
+{
+  auto& mesh = *model::meshes.at(mesh_);
+  return mesh.bin_label(bin) + " (born)";
+}
+
+//==============================================================================
+// C-API functions
+//==============================================================================
+
+extern "C" int openmc_meshborn_filter_get_mesh(
+  int32_t index, int32_t* index_mesh)
+{
+  return openmc_mesh_filter_get_mesh(index, index_mesh);
+}
+
+extern "C" int openmc_meshborn_filter_set_mesh(
+  int32_t index, int32_t index_mesh)
+{
+  return openmc_mesh_filter_set_mesh(index, index_mesh);
+}
+
+extern "C" int openmc_meshborn_filter_get_translation(
+  int32_t index, double translation[3])
+{
+  return openmc_mesh_filter_get_translation(index, translation);
+}
+
+extern "C" int openmc_meshborn_filter_set_translation(
+  int32_t index, double translation[3])
+{
+  return openmc_mesh_filter_set_translation(index, translation);
+}
+
+} // namespace openmc

src/tallies/filter_meshsurface.cpp

@@ -100,4 +100,16 @@ extern "C" int openmc_meshsurface_filter_set_mesh(
   return openmc_mesh_filter_set_mesh(index, index_mesh);
 }
 
+extern "C" int openmc_meshsurface_filter_get_translation(
+  int32_t index, double translation[3])
+{
+  return openmc_mesh_filter_get_translation(index, translation);
+}
+
+extern "C" int openmc_meshsurface_filter_set_translation(
+  int32_t index, double translation[3])
+{
+  return openmc_mesh_filter_set_translation(index, translation);
+}
+
 } // namespace openmc

src/tallies/tally.cpp

@@ -26,6 +26,7 @@
 #include "openmc/tallies/filter_energy.h"
 #include "openmc/tallies/filter_legendre.h"
 #include "openmc/tallies/filter_mesh.h"
+#include "openmc/tallies/filter_meshborn.h"
 #include "openmc/tallies/filter_meshsurface.h"
 #include "openmc/tallies/filter_particle.h"
 #include "openmc/tallies/filter_sph_harm.h"

tests/regression_tests/filter_meshborn/inputs_true.dat

@@ -0,0 +1,51 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1">
+      <density units="g/cm3" value="1.0"/>
+      <nuclide ao="1.0" name="H1"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1" universe="1"/>
+    <surface boundary="reflective" coeffs="0.0 0.0 0.0 10.0" id="1" type="sphere"/>
+  </geometry>
+  <settings>
+    <run_mode>fixed source</run_mode>
+    <particles>2000</particles>
+    <batches>8</batches>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="box">
+        <parameters>0.0 -10.0 -10.0 10.0 10.0 10.0</parameters>
+      </space>
+    </source>
+  </settings>
+  <tallies>
+    <mesh id="1">
+      <dimension>2 2 1</dimension>
+      <lower_left>-10.0 -10.0 -10.0</lower_left>
+      <upper_right>10.0 10.0 10.0</upper_right>
+    </mesh>
+    <filter id="1" type="mesh">
+      <bins>1</bins>
+    </filter>
+    <filter id="2" type="meshborn">
+      <bins>1</bins>
+    </filter>
+    <tally id="1" name="scatter">
+      <filters>1 2</filters>
+      <scores>scatter</scores>
+    </tally>
+    <tally id="2" name="scatter-mesh">
+      <filters>1</filters>
+      <scores>scatter</scores>
+    </tally>
+    <tally id="3" name="scatter-meshborn">
+      <filters>2</filters>
+      <scores>scatter</scores>
+    </tally>
+    <tally id="4" name="scatter-total">
+      <scores>scatter</scores>
+    </tally>
+  </tallies>
+</model>

tests/regression_tests/filter_meshborn/results_true.dat

@@ -0,0 +1,54 @@
+tally 1:
+0.000000E+00
+0.000000E+00
+2.631246E+01
+8.845079E+01
+0.000000E+00
+0.000000E+00
+2.450265E+00
+9.462266E-01
+0.000000E+00
+0.000000E+00
+3.878380E+02
+1.881752E+04
+0.000000E+00
+0.000000E+00
+2.932956E+01
+1.091674E+02
+0.000000E+00
+0.000000E+00
+1.837753E+00
+5.195343E-01
+0.000000E+00
+0.000000E+00
+2.944919E+01
+1.095819E+02
+0.000000E+00
+0.000000E+00
+2.921731E+01
+1.097387E+02
+0.000000E+00
+0.000000E+00
+4.019442E+02
+2.021184E+04
+tally 2:
+2.876273E+01
+1.060244E+02
+4.171676E+02
+2.176683E+04
+3.128695E+01
+1.238772E+02
+4.311615E+02
+2.325871E+04
+tally 3:
+0.000000E+00
+0.000000E+00
+4.452055E+02
+2.478148E+04
+0.000000E+00
+0.000000E+00
+4.631732E+02
+2.683862E+04
+tally 4:
+9.083787E+02
+1.031695E+05

tests/regression_tests/filter_meshborn/test.py

@@ -0,0 +1,107 @@
+"""Test the meshborn filter using a fixed source calculation on a H1 sphere.
+
+"""
+
+from numpy.testing import assert_allclose
+import numpy as np
+import openmc
+import pytest
+
+from tests.testing_harness import PyAPITestHarness
+
+
+RTOL = 1.0e-7
+ATOL = 0.0
+
+
+@pytest.fixture
+def model():
+    """Sphere of H1 with one hemisphere containing the source (x>0) and one
+    hemisphere with no source (x<0).
+
+    """
+    openmc.reset_auto_ids()
+    model = openmc.Model()
+
+    # Materials
+    h1 = openmc.Material()
+    h1.add_nuclide("H1", 1.0)
+    h1.set_density("g/cm3", 1.0)
+    model.materials = openmc.Materials([h1])
+
+    # Core geometry
+    r = 10.0
+    sphere = openmc.Sphere(r=r, boundary_type="reflective")
+    core = openmc.Cell(fill=h1, region=-sphere)
+    model.geometry = openmc.Geometry([core])
+
+    # Settings
+    model.settings.run_mode = 'fixed source'
+    model.settings.particles = 2000
+    model.settings.batches = 8
+    distribution = openmc.stats.Box((0., -r, -r), (r, r, r))
+    model.settings.source = openmc.IndependentSource(space=distribution)
+
+    # Tallies
+    mesh = openmc.RegularMesh()
+    mesh.dimension = (2, 2, 1)
+    mesh.lower_left = (-r, -r, -r)
+    mesh.upper_right = (r, r, r)
+
+    f_1 = openmc.MeshFilter(mesh)
+    f_2 = openmc.MeshBornFilter(mesh)
+
+    t_1 = openmc.Tally(name="scatter")
+    t_1.filters = [f_1, f_2]
+    t_1.scores = ["scatter"]
+
+    t_2 = openmc.Tally(name="scatter-mesh")
+    t_2.filters = [f_1]
+    t_2.scores = ["scatter"]
+
+    t_3 = openmc.Tally(name="scatter-meshborn")
+    t_3.filters = [f_2]
+    t_3.scores = ["scatter"]
+
+    t_4 = openmc.Tally(name="scatter-total")
+    t_4.scores = ["scatter"]
+
+    model.tallies = [t_1, t_2, t_3, t_4]
+
+    return model
+
+
+class MeshBornFilterTest(PyAPITestHarness):
+
+    def _compare_results(self):
+        """Additional unit tests on the tally results to check consistency."""
+        with openmc.StatePoint(self.statepoint_name) as sp:
+
+            t1 = sp.get_tally(name="scatter").mean.reshape(4, 4)
+            t2 = sp.get_tally(name="scatter-mesh").mean.reshape(4)
+            t3 = sp.get_tally(name="scatter-meshborn").mean.reshape(4)
+            t4 = sp.get_tally(name="scatter-total").mean.reshape(1)
+
+            # Consistency between mesh+meshborn matrix tally and meshborn tally
+            for i in range(4):
+                assert_allclose(t1[:, i].sum(), t3[i], rtol=RTOL, atol=ATOL)
+
+            # Consistency between mesh+meshborn matrix tally and mesh tally
+            for i in range(4):
+                assert_allclose(t1[i, :].sum(), t2[i], rtol=RTOL, atol=ATOL)
+
+            # Mesh cells in x<0 do not contribute to meshborn
+            assert_allclose(t1[:, 0].sum(), np.zeros(4), rtol=RTOL, atol=ATOL)
+            assert_allclose(t1[:, 2].sum(), np.zeros(4), rtol=RTOL, atol=ATOL)
+
+            # Consistency with total scattering
+            assert_allclose(t1.sum(), t4, rtol=RTOL, atol=ATOL)
+            assert_allclose(t2.sum(), t4, rtol=RTOL, atol=ATOL)
+            assert_allclose(t3.sum(), t4, rtol=RTOL, atol=ATOL)
+
+        super()._compare_results()
+
+
+def test_filter_meshborn(model):
+    harness = MeshBornFilterTest("statepoint.8.h5", model)
+    harness.main()

tests/unit_tests/test_filter_meshborn.py

@@ -0,0 +1,116 @@
+"""Test the meshborn filter using a fixed source calculation on a H1 sphere.
+
+"""
+
+import numpy as np
+from uncertainties import unumpy
+import openmc
+import pytest
+
+
+@pytest.fixture
+def model():
+    """Sphere of H1 with one hemisphere containing the source (x>0) and one
+    hemisphere with no source (x<0).
+
+    """
+    openmc.reset_auto_ids()
+    model = openmc.Model()
+
+    # Materials
+    h1 = openmc.Material()
+    h1.add_nuclide("H1", 1.0)
+    h1.set_density("g/cm3", 1.0)
+    model.materials = openmc.Materials([h1])
+
+    # Core geometry
+    r = 10.0
+    sphere = openmc.Sphere(r=r, boundary_type="reflective")
+    core = openmc.Cell(fill=h1, region=-sphere)
+    model.geometry = openmc.Geometry([core])
+
+    # Settings
+    model.settings.run_mode = 'fixed source'
+    model.settings.particles = 2000
+    model.settings.batches = 8
+
+    distribution = openmc.stats.Box((0., -r, -r), (r, r, r))
+    model.settings.source = openmc.IndependentSource(space=distribution)
+
+    # =============================================================================
+    # Tallies
+    # =============================================================================
+
+    mesh = openmc.RegularMesh()
+    mesh.dimension = (2, 2, 1)
+    mesh.lower_left = (-r, -r, -r)
+    mesh.upper_right = (r, r, r)
+
+    f = openmc.MeshBornFilter(mesh)
+    t_1 = openmc.Tally(name="scatter-collision")
+    t_1.filters = [f]
+    t_1.scores = ["scatter"]
+    t_1.estimator = "collision"
+
+    t_2 = openmc.Tally(name="scatter-tracklength")
+    t_2.filters = [f]
+    t_2.scores = ["scatter"]
+    t_2.estimator = "tracklength"
+
+    model.tallies = [t_1, t_2]
+
+    return model
+
+
+def test_estimator_consistency(model, run_in_tmpdir):
+    """Test that resuts obtained from a tracklength estimator are
+    consistent with results obtained from a collision estimator.
+
+    """
+    # Run OpenMC
+    sp_filename = model.run()
+
+    # Get radial flux distribution
+    with openmc.StatePoint(sp_filename) as sp:
+        scatter_collision = sp.get_tally(name="scatter-collision").mean.ravel()
+        scatter_collision_std_dev = sp.get_tally(name="scatter-collision").std_dev.ravel()
+        scatter_tracklength = sp.get_tally(name="scatter-tracklength").mean.ravel()
+        scatter_tracklength_std_dev = sp.get_tally(name="scatter-tracklength").std_dev.ravel()
+
+    collision = unumpy.uarray(scatter_collision, scatter_collision_std_dev)
+    tracklength = unumpy.uarray(scatter_tracklength, scatter_tracklength_std_dev)
+    delta = abs(collision - tracklength)
+
+    diff = unumpy.nominal_values(delta)
+    std_dev = unumpy.std_devs(delta)
+    assert np.all(diff <= 3 * std_dev)
+
+
+def test_xml_serialization():
+    """Test xml serialization of the meshborn filter."""
+    openmc.reset_auto_ids()
+
+    mesh = openmc.RegularMesh()
+    mesh.dimension = (1, 1, 1)
+    mesh.lower_left = (0.0, 0.0, 0.0)
+    mesh.upper_right = (1.0, 1.0, 1.0)
+
+    filter = openmc.MeshBornFilter(mesh)
+    filter.translation = (2.0, 2.0, 2.0)
+    assert filter.mesh.id == 1
+    assert filter.mesh.dimension == (1, 1, 1)
+    assert filter.mesh.lower_left == (0.0, 0.0, 0.0)
+    assert filter.mesh.upper_right == (1.0, 1.0, 1.0)
+
+    repr(filter)
+
+    elem = filter.to_xml_element()
+    assert elem.tag == 'filter'
+    assert elem.attrib['type'] == 'meshborn'
+    assert elem[0].text == "1"
+    assert elem.get("translation") == "2.0 2.0 2.0"
+
+    meshes = {1: mesh}
+    new_filter = openmc.Filter.from_xml_element(elem, meshes=meshes)
+    assert new_filter.bins == filter.bins
+    np.testing.assert_equal(new_filter.translation, [2.0, 2.0, 2.0])

tests/unit_tests/test_tallies.py

@@ -10,8 +10,9 @@ def test_xml_roundtrip(run_in_tmpdir):
     mesh.upper_right = (10., 10., 10.,)
     mesh.dimension = (5, 5, 5)
     mesh_filter = openmc.MeshFilter(mesh)
+    meshborn_filter = openmc.MeshBornFilter(mesh)
     tally = openmc.Tally()
-    tally.filters = [mesh_filter]
+    tally.filters = [mesh_filter, meshborn_filter]
     tally.nuclides = ['U235', 'I135', 'Li6']
     tally.scores = ['total', 'fission', 'heating']
     tally.derivative = openmc.TallyDerivative(
@@ -27,9 +28,11 @@ def test_xml_roundtrip(run_in_tmpdir):
     assert len(new_tallies) == 1
     new_tally = new_tallies[0]
     assert new_tally.id == tally.id
-    assert len(new_tally.filters) == 1
+    assert len(new_tally.filters) == 2
     assert isinstance(new_tally.filters[0], openmc.MeshFilter)
     assert np.allclose(new_tally.filters[0].mesh.lower_left, mesh.lower_left)
+    assert isinstance(new_tally.filters[1], openmc.MeshBornFilter)
+    assert np.allclose(new_tally.filters[1].mesh.lower_left, mesh.lower_left)
     assert new_tally.nuclides == tally.nuclides
     assert new_tally.scores == tally.scores
     assert new_tally.derivative.variable == tally.derivative.variable