Make MaterialFilter support void material

openmc/filter.py

@@ -468,7 +468,7 @@ class UniverseFilter(WithIDFilter):
     expected_type = UniverseBase
 
 
-class MaterialFilter(WithIDFilter):
+class MaterialFilter(Filter):
     """Bins tally event locations based on the Material they occurred in.
 
     Parameters
@@ -489,7 +489,34 @@ class MaterialFilter(WithIDFilter):
         The number of filter bins
 
     """
-    expected_type = Material
+    expected_type = Material | str | None
+
+    def __init__(self, bins, filter_id=None):
+        bins = np.atleast_1d(bins)
+
+        # Make sure bins are either integers or appropriate objects
+        cv.check_iterable_type('filter bins', bins,
+                               (Integral, self.expected_type))
+        for i,b in enumerate(bins):
+            if isinstance(b, str):
+              cv.check_value(f'filter bins[{i}]', b, "void")
+
+        def get_id(b):
+            if isinstance(b, Integral):
+                return b
+            elif b == 'void' or b is None:
+                return -1
+            else:
+                return b.id
+
+        # Extract ID values
+        bins = np.array([get_id(b) for b in bins])
+        super().__init__(bins, filter_id)
+
+    def check_bins(self, bins):
+        # Check the bin values.
+        for edge in bins:
+            cv.check_greater_than('filter bin', edge, -1, equality=True)        
 
 
 class MaterialFromFilter(WithIDFilter):

src/tallies/filter_material.cpp

@@ -15,6 +15,10 @@ void MaterialFilter::from_xml(pugi::xml_node node)
   // Get material IDs and convert to indices in the global materials vector
   auto mats = get_node_array<int32_t>(node, "bins");
   for (auto& m : mats) {
+    // If void material skip check
+    if (m == MATERIAL_VOID)
+      continue;
+
     auto search = model::material_map.find(m);
     if (search == model::material_map.end()) {
       throw std::runtime_error {fmt::format(
@@ -33,10 +37,12 @@ void MaterialFilter::set_materials(span<const int32_t> materials)
   materials_.reserve(materials.size());
   map_.clear();
 
+  int32_t size = model::materials.size();
+
   // Update materials and mapping
   for (auto& index : materials) {
-    assert(index >= 0);
-    assert(index < model::materials.size());
+    assert(index >= MATERIAL_VOID);
+    assert(index < size);
     materials_.push_back(index);
     map_[index] = materials_.size() - 1;
   }
@@ -58,14 +64,21 @@ void MaterialFilter::to_statepoint(hid_t filter_group) const
 {
   Filter::to_statepoint(filter_group);
   vector<int32_t> material_ids;
-  for (auto c : materials_)
-    material_ids.push_back(model::materials[c]->id_);
+  for (auto c : materials_) {
+    auto id = MATERIAL_VOID;
+    if (c != MATERIAL_VOID)
+      id = model::materials[c]->id_;
+    material_ids.push_back(id);
+  }
   write_dataset(filter_group, "bins", material_ids);
 }
 
 std::string MaterialFilter::text_label(int bin) const
 {
-  return fmt::format("Material {}", model::materials[materials_[bin]]->id_);
+  std::string label = "Material Void";
+  if (materials_[bin] != MATERIAL_VOID)
+    label = fmt::format("Material {}", model::materials[materials_[bin]]->id_);
+  return label;
 }
 
 //==============================================================================

tests/unit_tests/test_filters.py

@@ -334,6 +334,16 @@ def test_weight():
     assert new_f.id == f.id
     assert np.allclose(new_f.bins, f.bins)
 
+def test_material():
+    filt1 = openmc.MaterialFilter([None])
+    filt2 = openmc.MaterialFilter(['void'])
+
+    assert filt1.bins == [-1]
+    assert filt2.bins == [-1]
+
+    with raises(ValueError):
+        openmc.MaterialFilter(['void2'])
+
 
 def test_mesh_material():
     mat1 = openmc.Material()

tests/unit_tests/test_void.py

@@ -0,0 +1,50 @@
+import numpy as np
+import openmc
+import pytest
+
+
+@pytest.fixture
+def model():
+    openmc.reset_auto_ids()
+    model = openmc.model.Model()
+
+    zn = openmc.Material()
+    zn.set_density('g/cm3', 7.14)
+    zn.add_nuclide('Zn64', 1.0)
+    model.materials.append(zn)
+
+    radii = np.linspace(1.0, 100.0)
+    surfs = [openmc.Sphere(r=r) for r in radii]
+    surfs[-1].boundary_type = 'vacuum'
+    cells = [openmc.Cell(fill=(zn if i % 2 == 0 else None), region=region)
+             for i, region in enumerate(openmc.model.subdivide(surfs))]
+    model.geometry = openmc.Geometry(cells)
+
+    model.settings.run_mode = 'fixed source'
+    model.settings.batches = 3
+    model.settings.particles = 1000
+    model.settings.source = openmc.IndependentSource(space=openmc.stats.Point())
+
+    cell_filter = openmc.CellFilter(cells[1::2])
+    material_filter = openmc.MaterialFilter([None])
+
+    tally1 = openmc.Tally()
+    tally1.filters = [cell_filter]
+    tally1.scores = ['flux']
+
+    tally2 = openmc.Tally()
+    tally2.filters = [material_filter]
+    tally2.scores = ['flux']
+
+    model.tallies.append(tally1)
+    model.tallies.append(tally2)
+
+    return model
+
+def test_equivalent_void_specification(model, run_in_tmpdir):
+    sp_file = model.run()
+    with openmc.StatePoint(sp_file) as sp:
+        tally1 = sp.tallies[1]
+        tally2 = sp.tallies[2]
+
+    assert np.isclose(tally1.mean.sum(),tally2.mean.sum(), rtol=1e-10, atol=0)