Merge pull request #1863 from pshriwise/cell_inst_filter

Enable cell instance tallies for universe/lattice-filled cells

include/openmc/geometry.h

@@ -40,6 +40,17 @@ inline bool coincident(double d1, double d2) {
 
 bool check_cell_overlap(Particle& p, bool error=true);
 
+//==============================================================================
+//! Get the cell instance for a particle at the specified universe level
+//!
+//! \param p A particle for which to compute the instance using
+//!   its coordinates
+//! \param level The level (zero indexed) of the geometry where the instance should be computed.
+//! \return The instance of the cell at the specified level.
+//==============================================================================
+
+int cell_instance_at_level(const Particle& p, int level);
+
 //==============================================================================
 //! Locate a particle in the geometry tree and set its geometry data fields.
 //!

include/openmc/tallies/filter_cell_instance.h

@@ -53,8 +53,14 @@ private:
   //! The indices of the cells binned by this filter.
   vector<CellInstance> cell_instances_;
 
+  //! The set of cells used in this filter
+  std::unordered_set<int32_t> cells_;
+
   //! A map from cell/instance indices to filter bin indices.
   std::unordered_map<CellInstance, gsl::index, CellInstanceHash> map_;
+
+  //! Indicates if filter uses only material-filled cells
+  bool material_cells_only_;
 };
 
 } // namespace openmc

src/geometry.cpp

@@ -65,6 +65,37 @@ bool check_cell_overlap(Particle& p, bool error)
 
 //==============================================================================
 
+int cell_instance_at_level(const Particle& p, int level) {
+  // throw error if the requested level is too deep for the geometry
+  if (level > model::n_coord_levels) {
+    fatal_error(fmt::format("Cell instance at level {} requested, but only {} levels exist in the geometry.", level, p.n_coord()));
+  }
+
+  // determine the cell instance
+  Cell& c {*model::cells[p.coord(level).cell]};
+
+  // quick exit if this cell doesn't have distribcell instances
+  if (c.distribcell_index_ == C_NONE) return C_NONE;
+
+  // compute the cell's instance
+  int instance = 0;
+  for (int i = 0; i < level; i++) {
+    const auto& c_i {*model::cells[p.coord(i).cell]};
+    if (c_i.type_ == Fill::UNIVERSE) {
+      instance += c_i.offset_[c.distribcell_index_];
+    } else if (c_i.type_ == Fill::LATTICE) {
+      auto& lat {*model::lattices[p.coord(i + 1).lattice]};
+      const auto& i_xyz {p.coord(i + 1).lattice_i};
+      if (lat.are_valid_indices(i_xyz)) {
+        instance += lat.offset(c.distribcell_index_, i_xyz);
+      }
+    }
+  }
+  return instance;
+}
+
+//==============================================================================
+
 bool
 find_cell_inner(Particle& p, const NeighborList* neighbor_list)
 {
@@ -129,23 +160,11 @@ find_cell_inner(Particle& p, const NeighborList* neighbor_list)
     if (c.type_ == Fill::MATERIAL) {
       // Found a material cell which means this is the lowest coord level.
 
+      p.cell_instance() = 0;
       // Find the distribcell instance number.
-      int offset = 0;
       if (c.distribcell_index_ >= 0) {
-        for (int i = 0; i < p.n_coord(); i++) {
-          const auto& c_i {*model::cells[p.coord(i).cell]};
-          if (c_i.type_ == Fill::UNIVERSE) {
-            offset += c_i.offset_[c.distribcell_index_];
-          } else if (c_i.type_ == Fill::LATTICE) {
-            auto& lat {*model::lattices[p.coord(i + 1).lattice]};
-            const auto& i_xyz {p.coord(i + 1).lattice_i};
-            if (lat.are_valid_indices(i_xyz)) {
-              offset += lat.offset(c.distribcell_index_, i_xyz);
-            }
-          }
-        }
+        p.cell_instance() = cell_instance_at_level(p, p.n_coord() - 1);
       }
-      p.cell_instance() = offset;
 
       // Set the material and temperature.
       p.material_last() = p.material();

src/geometry_aux.cpp

@@ -319,9 +319,13 @@ prepare_distribcell()
   std::unordered_set<int32_t> distribcells;
   for (auto& filt : model::tally_filters) {
     auto* distrib_filt = dynamic_cast<DistribcellFilter*>(filt.get());
+    auto* cell_inst_filt = dynamic_cast<CellInstanceFilter*>(filt.get());
     if (distrib_filt) {
       distribcells.insert(distrib_filt->cell());
     }
+    if (cell_inst_filt) {
+      for (const auto& c_inst : cell_inst_filt->cell_instances()) distribcells.insert(c_inst.index_cell);
+    }
   }
 
   // By default, add material cells to the list of distributed cells

src/plot.cpp

@@ -47,7 +47,7 @@ IdData::set_value(size_t y, size_t x, const Particle& p, int level) {
     data_(y, x, 1) = NOT_FOUND;
   } else {
     data_(y, x, 0) = model::cells.at(p.coord(level).cell)->id_;
-    data_(y, x, 1) = p.cell_instance();
+    data_(y, x, 1) = level == p.n_coord() - 1 ?  p.cell_instance() : cell_instance_at_level(p, level);
   }
 
   // set material data

src/tallies/filter_cell_instance.cpp

@@ -7,6 +7,7 @@
 #include "openmc/capi.h"
 #include "openmc/cell.h"
 #include "openmc/error.h"
+#include "openmc/geometry.h"
 #include "openmc/xml_interface.h"
 
 namespace openmc {
@@ -46,23 +47,27 @@ CellInstanceFilter::set_cell_instances(gsl::span<CellInstance> instances)
   // Clear existing cells
   cell_instances_.clear();
   cell_instances_.reserve(instances.size());
+  cells_.clear();
   map_.clear();
 
   // Update cells and mapping
   for (auto& x : instances) {
     Expects(x.index_cell >= 0);
     Expects(x.index_cell < model::cells.size());
-    const auto& c {model::cells[x.index_cell]};
-    if (c->type_ != Fill::MATERIAL) {
-      throw std::invalid_argument{fmt::format(
-        "Cell {} is not filled with a material. Only material cells can be "
-        "used in a cell instance filter.", c->id_)};
-    }
     cell_instances_.push_back(x);
+    cells_.insert(x.index_cell);
     map_[x] = cell_instances_.size() - 1;
   }
 
   n_bins_ = cell_instances_.size();
+
+  material_cells_only_ = true;
+  for (const auto& cell_inst : cell_instances_) {
+    const auto& c = *model::cells[cell_inst.index_cell];
+    if (c.type_ == Fill::MATERIAL) continue;
+    material_cells_only_ = false;
+    break;
+  }
 }
 
 void
@@ -71,11 +76,30 @@ CellInstanceFilter::get_all_bins(const Particle& p, TallyEstimator estimator,
 {
   gsl::index index_cell = p.coord(p.n_coord() - 1).cell;
   gsl::index instance = p.cell_instance();
-  auto search = map_.find({index_cell, instance});
-  if (search != map_.end()) {
-    int index_bin = search->second;
-    match.bins_.push_back(index_bin);
-    match.weights_.push_back(1.0);
+
+  if (cells_.count(index_cell) > 0) {
+    auto search = map_.find({index_cell, instance});
+    if (search != map_.end()) {
+      int index_bin = search->second;
+      match.bins_.push_back(index_bin);
+      match.weights_.push_back(1.0);
+    }
+  }
+
+  if (material_cells_only_) return;
+
+  for (int i = 0; i < p.n_coord() - 1; i++) {
+    gsl::index index_cell = p.coord(i).cell;
+    // if this cell isn't used on the filter, move on
+    if (cells_.count(index_cell) == 0) continue;
+
+    // if this cell is used in the filter, check the instance as well
+    gsl::index instance = cell_instance_at_level(p, i);
+    auto search = map_.find({index_cell, instance});
+    if (search != map_.end()) {
+      match.bins_.push_back(search->second);
+      match.weights_.push_back(1.0);
+    }
   }
 }
 

tests/regression_tests/filter_cellinstance/inputs_true.dat

@@ -1,19 +1,20 @@
 <?xml version='1.0' encoding='utf-8'?>
 <geometry>
-  <cell id="1" material="1" region="-1" universe="1" />
+  <cell id="1" material="1" region="-1" universe="2" />
   <cell id="2" material="2" region="1" universe="1" />
-  <cell id="3" material="1" region="-2" universe="2" />
-  <cell id="4" material="2" region="2" universe="2" />
-  <cell fill="3" id="5" region="3 -4 5 -6" universe="4" />
-  <lattice id="3">
+  <cell fill="1" id="3" universe="2" />
+  <cell id="4" material="1" region="-2" universe="3" />
+  <cell id="5" material="2" region="2" universe="3" />
+  <cell fill="4" id="6" region="3 -4 5 -6" universe="5" />
+  <lattice id="4">
     <pitch>2 2</pitch>
     <dimension>4 4</dimension>
     <lower_left>-4 -4</lower_left>
     <universes>
-1 2 2 2 
-2 1 2 2 
-2 2 1 2 
-2 2 2 1 </universes>
+2 3 3 3 
+3 2 3 3 
+3 3 2 3 
+3 3 3 2 </universes>
   </lattice>
   <surface coeffs="0.0 0.0 0.7" id="1" type="z-cylinder" />
   <surface coeffs="0.0 0.0 0.5" id="2" type="z-cylinder" />
@@ -48,10 +49,10 @@
 <?xml version='1.0' encoding='utf-8'?>
 <tallies>
   <filter id="1" type="cellinstance">
-    <bins>3 0 3 1 3 2 3 3 3 4 3 5 3 6 3 7 3 8 3 9 3 10 3 11 2 0 2 1 2 2 2 3</bins>
+    <bins>4 0 4 1 4 2 4 3 4 4 4 5 4 6 4 7 4 8 4 9 4 10 4 11 2 0 2 1 2 2 2 3 3 0 3 1 3 2 3 3</bins>
   </filter>
   <filter id="2" type="cellinstance">
-    <bins>2 3 2 2 2 1 2 0 3 11 3 10 3 9 3 8 3 7 3 6 3 5 3 4 3 3 3 2 3 1 3 0</bins>
+    <bins>3 3 3 2 3 1 3 0 2 3 2 2 2 1 2 0 4 11 4 10 4 9 4 8 4 7 4 6 4 5 4 4 4 3 4 2 4 1 4 0</bins>
   </filter>
   <tally id="1">
     <filters>1</filters>

tests/regression_tests/filter_cellinstance/results_true.dat

@@ -33,6 +33,14 @@ tally 1:
 2.015188E+02
 9.951936E+00
 2.345815E+01
+1.061333E+01
+2.508486E+01
+2.767786E+01
+1.781264E+02
+2.898957E+01
+2.015188E+02
+9.951936E+00
+2.345815E+01
 tally 2:
 9.951936E+00
 2.345815E+01
@@ -42,6 +50,14 @@ tally 2:
 1.781264E+02
 1.061333E+01
 2.508486E+01
+9.951936E+00
+2.345815E+01
+2.898957E+01
+2.015188E+02
+2.767786E+01
+1.781264E+02
+1.061333E+01
+2.508486E+01
 6.794395E-02
 1.135006E-03
 1.299309E-01

tests/regression_tests/filter_cellinstance/test.py

@@ -1,3 +1,4 @@
+from numpy.testing import assert_array_almost_equal
 import openmc
 import openmc.model
 import pytest
@@ -5,6 +6,24 @@ import pytest
 from tests.testing_harness import PyAPITestHarness
 
 
+class CellInstanceFilterTest(PyAPITestHarness):
+
+    def _compare_results(self):
+        with openmc.StatePoint(self.statepoint_name) as sp:
+            # we expect the tally results for the instances of
+            # cells 2 and 3 to be the same as 2 is nested
+            # in a universe directly under 3
+            t1 = sp.tallies[1]
+            f1 = sp.filters[1]
+            c2_bins = [tuple(tuple(i) for i in f1.bins if i[0] == 2)]
+            c2_mean = t1.get_values(filters=[openmc.CellInstanceFilter], filter_bins=c2_bins)
+            c3_bins = [tuple(tuple(i) for i in f1.bins if i[0] == 3)]
+            c3_mean = t1.get_values(filters=[openmc.CellInstanceFilter], filter_bins=c3_bins)
+            assert_array_almost_equal(c2_mean, c3_mean)
+
+        return super()._compare_results()
+
+
 @pytest.fixture
 def model():
     model = openmc.model.Model()
@@ -22,21 +41,24 @@ def model():
     cyl1 = openmc.ZCylinder(r=0.7)
     c1 = openmc.Cell(fill=m1, region=-cyl1)
     c2 = openmc.Cell(fill=m2, region=+cyl1)
-    u1 = openmc.Universe(cells=[c1, c2])
+    # intermediate universe containing only cell 2
+    u1 = openmc.Universe(cells=[c2])
+    c3 = openmc.Cell(fill=u1)
+    u2 = openmc.Universe(cells=[c1, c3])
 
     cyl2 = openmc.ZCylinder(r=0.5)
-    c3 = openmc.Cell(fill=m1, region=-cyl2)
-    c4 = openmc.Cell(fill=m2, region=+cyl2)
-    u2 = openmc.Universe(cells=[c3, c4])
+    c4 = openmc.Cell(fill=m1, region=-cyl2)
+    c5 = openmc.Cell(fill=m2, region=+cyl2)
+    u3 = openmc.Universe(cells=[c4, c5])
 
     lat = openmc.RectLattice()
     lat.lower_left = (-4, -4)
     lat.pitch = (2, 2)
     lat.universes = [
-        [u1, u2, u2, u2],
-        [u2, u1, u2, u2],
-        [u2, u2, u1, u2],
-        [u2, u2, u2, u1]
+        [u2, u3, u3, u3],
+        [u3, u2, u3, u3],
+        [u3, u3, u2, u3],
+        [u3, u3, u3, u2]
     ]
     box = openmc.model.rectangular_prism(8.0, 8.0, boundary_type='reflective')
     main_cell = openmc.Cell(fill=lat, region=box)
@@ -49,8 +71,9 @@ def model():
     model.settings.particles = 1000
     model.settings.source = openmc.Source(space=openmc.stats.Point())
 
-    instances = ([(c3, i) for i in range(c3.num_instances)] +
-                 [(c2, i) for i in range(c2.num_instances)])
+    instances = ([(c4, i) for i in range(c4.num_instances)] +
+                 [(c2, i) for i in range(c2.num_instances)] +
+                 [(c3, i) for i in range(c3.num_instances)])
     f1 = openmc.CellInstanceFilter(instances)
     f2 = openmc.CellInstanceFilter(instances[::-1])
     t1 = openmc.Tally()
@@ -65,5 +88,5 @@ def model():
 
 
 def test_cell_instance(model):
-    harness = PyAPITestHarness('statepoint.5.h5', model)
+    harness = CellInstanceFilterTest('statepoint.5.h5', model)
     harness.main()

tests/testing_harness.py

@@ -115,6 +115,10 @@ class TestHarness:
 
         return outstr
 
+    @property
+    def statepoint_name(self):
+        return self._sp_name
+
     def _write_results(self, results_string):
         """Write the results to an ASCII file."""
         with open('results_test.dat', 'w') as fh: