Implement cell instance filter

CMakeLists.txt

@@ -253,6 +253,7 @@ list(APPEND libopenmc_SOURCES
   src/tallies/filter_cellborn.cpp
   src/tallies/filter_cellfrom.cpp
   src/tallies/filter_cell.cpp
+  src/tallies/filter_cell_instance.cpp
   src/tallies/filter_delayedgroup.cpp
   src/tallies/filter_distribcell.cpp
   src/tallies/filter_energyfunc.cpp

include/openmc/tallies/filter_cell_instance.h

@@ -0,0 +1,92 @@
+#ifndef OPENMC_TALLIES_FILTER_CELL_INSTANCE_H
+#define OPENMC_TALLIES_FILTER_CELL_INSTANCE_H
+
+#include <cstdint>
+#include <functional> // for hash
+#include <unordered_map>
+#include <vector>
+
+#include <gsl/gsl>
+
+#include "openmc/tallies/filter.h"
+
+//==============================================================================
+//! Define an instance of a particular cell
+//==============================================================================
+
+namespace openmc {
+struct CellInstance {
+  //! Check for equality
+  bool operator==(const CellInstance& other) const
+  { return index_cell == other.index_cell && instance == other.instance; }
+
+  gsl::index index_cell;
+  gsl::index instance;
+};
+} // namespace openmc
+
+namespace std {
+template<>
+struct hash<openmc::CellInstance> {
+  // Taken from https://stackoverflow.com/a/17017281
+  std::size_t operator()(const openmc::CellInstance& k) const
+  {
+    std::size_t res = 17;
+    res = 31 * res + std::hash<gsl::index>()(k.index_cell);
+    res = 31 * res + std::hash<gsl::index>()(k.instance);
+    return res;
+  }
+};
+} // namespace std
+
+
+namespace openmc {
+
+//==============================================================================
+//! Specifies cell instances that tally events reside in.
+//==============================================================================
+
+class CellInstanceFilter : public Filter {
+public:
+  //----------------------------------------------------------------------------
+  // Constructors, destructors
+
+  //CellInstanceFilter() = default;
+  //CellInstanceFilter(gsl::span<CellInstance> instances);
+  ~CellInstanceFilter() = default;
+
+  //----------------------------------------------------------------------------
+  // Methods
+
+  std::string type() const override {return "cellinstance";}
+
+  void from_xml(pugi::xml_node node) override;
+
+  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  //----------------------------------------------------------------------------
+  // Accessors
+
+  const std::vector<CellInstance>& cell_instances() const { return cell_instances_; }
+
+  void set_cell_instances(gsl::span<CellInstance> instances);
+
+private:
+  //----------------------------------------------------------------------------
+  // Data members
+
+  //! The indices of the cells binned by this filter.
+  std::vector<CellInstance> cell_instances_;
+
+  //! A map from cell/instance indices to filter bin indices.
+  std::unordered_map<CellInstance, gsl::index> map_;
+};
+
+} // namespace openmc
+
+#endif // OPENMC_TALLIES_FILTER_CELL_INSTANCE_H

src/geometry.cpp

@@ -132,8 +132,8 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
       //! Found a material cell which means this is the lowest coord level.
 
       // Find the distribcell instance number.
-      if (c.material_.size() > 1 || c.sqrtkT_.size() > 1) {
-        int offset = 0;
+      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) {
@@ -148,10 +148,8 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
             }
           }
         }
-        p->cell_instance_ = offset;
-      } else {
-        p->cell_instance_ = 0;
       }
+      p->cell_instance_ = offset;
 
       // Set the material and temperature.
       p->material_last_ = p->material_;

src/geometry_aux.cpp

@@ -18,6 +18,7 @@
 #include "openmc/settings.h"
 #include "openmc/surface.h"
 #include "openmc/tallies/filter.h"
+#include "openmc/tallies/filter_cell_instance.h"
 #include "openmc/tallies/filter_distribcell.h"
 
 
@@ -319,13 +320,20 @@ find_root_universe()
 void
 prepare_distribcell()
 {
-  // Find all cells listed in a DistribcellFilter.
+  // Find all cells listed in a DistribcellFilter or CellInstanceFilter
   std::unordered_set<int32_t> distribcells;
   for (auto& filt : model::tally_filters) {
     auto* distrib_filt = dynamic_cast<DistribcellFilter*>(filt.get());
     if (distrib_filt) {
       distribcells.insert(distrib_filt->cell());
     }
+
+    auto* ci_filt = dynamic_cast<CellInstanceFilter*>(filt.get());
+    if (ci_filt) {
+      for (const auto& i : ci_filt->cell_instances()) {
+        distribcells.insert(i.index_cell);
+      }
+    }
   }
 
   // Find all cells with distributed materials or temperatures.  Make sure that

src/tallies/filter.cpp

@@ -12,6 +12,7 @@
 #include "openmc/tallies/filter_cell.h"
 #include "openmc/tallies/filter_cellborn.h"
 #include "openmc/tallies/filter_cellfrom.h"
+#include "openmc/tallies/filter_cell_instance.h"
 #include "openmc/tallies/filter_delayedgroup.h"
 #include "openmc/tallies/filter_distribcell.h"
 #include "openmc/tallies/filter_energyfunc.h"
@@ -100,6 +101,8 @@ Filter* Filter::create(const std::string& type, int32_t id)
     model::tally_filters.push_back(std::make_unique<CellbornFilter>());
   } else if (type == "cellfrom") {
     model::tally_filters.push_back(std::make_unique<CellFromFilter>());
+  } else if (type == "cellinstance") {
+    model::tally_filters.push_back(std::make_unique<CellInstanceFilter>());
   } else if (type == "distribcell") {
     model::tally_filters.push_back(std::make_unique<DistribcellFilter>());
   } else if (type == "delayedgroup") {

src/tallies/filter_cell_instance.cpp

@@ -0,0 +1,94 @@
+#include "openmc/tallies/filter_cell_instance.h"
+
+#include <sstream>
+
+#include "openmc/capi.h"
+#include "openmc/cell.h"
+#include "openmc/error.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+void
+CellInstanceFilter::from_xml(pugi::xml_node node)
+{
+  // Get cell IDs/instances
+  auto cells = get_node_array<int32_t>(node, "bins");
+  Expects(cells.size() % 2 == 0);
+
+  // Convert into vector of CellInstance
+  std::vector<CellInstance> instances;
+  for (gsl::index i = 0; i < cells.size() / 2; ++i) {
+    int32_t cell_id = cells[2*i];
+    gsl::index instance = cells[2*i + 1];
+    auto search = model::cell_map.find(cell_id);
+    if (search == model::cell_map.end()) {
+      std::stringstream err_msg;
+      err_msg << "Could not find cell " << cell_id
+              << " specified on tally filter.";
+      throw std::runtime_error{err_msg.str()};
+    }
+    gsl::index index = search->second;
+    instances.push_back({index, instance});
+  }
+
+  this->set_cell_instances(instances);
+}
+
+void
+CellInstanceFilter::set_cell_instances(gsl::span<CellInstance> instances)
+{
+  // Clear existing cells
+  cell_instances_.clear();
+  cell_instances_.reserve(instances.size());
+  map_.clear();
+
+  // Update cells and mapping
+  for (auto& x : instances) {
+    Expects(x.index_cell >= 0);
+    Expects(x.index_cell < model::cells.size());
+    cell_instances_.push_back(x);
+    map_[x] = cell_instances_.size() - 1;
+  }
+
+  n_bins_ = cell_instances_.size();
+}
+
+void
+CellInstanceFilter::get_all_bins(const Particle* p, int estimator,
+                         FilterMatch& match) const
+{
+  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);
+  }
+}
+
+void
+CellInstanceFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  size_t n = cell_instances_.size();
+  xt::xtensor<size_t, 2> data({n, 2});
+  for (gsl::index i = 0; i < n; ++i) {
+    const auto& x = cell_instances_[i];
+    data(i, 0) = model::cells[x.index_cell]->id_;
+    data(i, 1) = x.instance;
+  }
+  write_dataset(filter_group, "bins", data);
+}
+
+std::string
+CellInstanceFilter::text_label(int bin) const
+{
+  const auto& x = cell_instances_[bin];
+  auto cell_id = model::cells[x.index_cell]->id_;
+  return "Cell " + std::to_string(cell_id) + ", CellInstance "
+    + std::to_string(x.instance);
+}
+
+} // namespace openmc