Correctly score pulse height tally when no cell filter is present (#3821)

Co-authored-by: Paul Romano <[email protected]>

Author: GuySten

include/openmc/tallies/tally.h

@@ -212,7 +212,7 @@ extern vector<int> active_collision_tallies;
 extern vector<int> active_meshsurf_tallies;
 extern vector<int> active_surface_tallies;
 extern vector<int> active_pulse_height_tallies;
-extern vector<int> pulse_height_cells;
+extern vector<int32_t> pulse_height_cells;
 extern vector<double> time_grid;
 
 } // namespace model

src/tallies/tally.cpp

@@ -2,6 +2,7 @@
 
 #include "openmc/array.h"
 #include "openmc/capi.h"
+#include "openmc/cell.h"
 #include "openmc/constants.h"
 #include "openmc/container_util.h"
 #include "openmc/error.h"
@@ -62,7 +63,7 @@ vector<int> active_collision_tallies;
 vector<int> active_meshsurf_tallies;
 vector<int> active_surface_tallies;
 vector<int> active_pulse_height_tallies;
-vector<int> pulse_height_cells;
+vector<int32_t> pulse_height_cells;
 vector<double> time_grid;
 } // namespace model
 
@@ -632,29 +633,24 @@ void Tally::set_scores(const vector<std::string>& scores)
         estimator_ = TallyEstimator::COLLISION;
       break;
 
-    case SCORE_PULSE_HEIGHT:
+    case SCORE_PULSE_HEIGHT: {
       if (non_cell_energy_present) {
         fatal_error("Pulse-height tallies are not compatible with filters "
                     "other than CellFilter and EnergyFilter");
       }
       type_ = TallyType::PULSE_HEIGHT;
-
-      // Collecting indices of all cells covered by the filters in the pulse
-      // height tally in global variable pulse_height_cells
-      for (const auto& i_filt : filters_) {
-        auto cell_filter =
-          dynamic_cast<CellFilter*>(model::tally_filters[i_filt].get());
-        if (cell_filter) {
-          const auto& cells = cell_filter->cells();
-          for (int i = 0; i < cell_filter->n_bins(); i++) {
-            int cell_index = cells[i];
-            if (!contains(model::pulse_height_cells, cell_index)) {
-              model::pulse_height_cells.push_back(cell_index);
-            }
-          }
-        }
+      // Collect all unique cell indices covered by this tally.
+      // If no CellFilter is present, all cells in the geometry are scored.
+      const auto* cell_filter_ptr = get_filter<CellFilter>();
+      int n = cell_filter_ptr ? cell_filter_ptr->n_bins()
+                              : static_cast<int>(model::cells.size());
+      for (int i = 0; i < n; ++i) {
+        int32_t cell_index = cell_filter_ptr ? cell_filter_ptr->cells()[i] : i;
+        if (!contains(model::pulse_height_cells, cell_index))
+          model::pulse_height_cells.push_back(cell_index);
       }
       break;
+    }
 
     case SCORE_IFP_TIME_NUM:
     case SCORE_IFP_BETA_NUM:

src/tallies/tally_scoring.cpp

@@ -2672,56 +2672,49 @@ void score_pulse_height_tally(Particle& p, const vector<int>& tallies)
   for (auto i_tally : tallies) {
     auto& tally {*model::tallies[i_tally]};
 
-    // Determine all CellFilter in the tally
-    for (const auto& filter : tally.filters()) {
-      auto cell_filter =
-        dynamic_cast<CellFilter*>(model::tally_filters[filter].get());
-      if (cell_filter != nullptr) {
-
-        const auto& cells = cell_filter->cells();
-        // Loop over all cells in the CellFilter
-        for (auto cell_index = 0; cell_index < cells.size(); ++cell_index) {
-          int cell_id = cells[cell_index];
-
-          // Temporarily change cell of particle
-          p.n_coord() = 1;
-          p.coord(0).cell() = cell_id;
-
-          // Determine index of cell in model::pulse_height_cells
-          auto it = std::find(model::pulse_height_cells.begin(),
-            model::pulse_height_cells.end(), cell_id);
-          int index = std::distance(model::pulse_height_cells.begin(), it);
-
-          // Temporarily change energy of particle to pulse-height value
-          p.E_last() = p.pht_storage()[index];
-
-          // Initialize an iterator over valid filter bin combinations. If
-          // there are no valid combinations, use a continue statement to ensure
-          // we skip the assume_separate break below.
-          auto filter_iter = FilterBinIter(tally, p);
-          auto end = FilterBinIter(tally, true, &p.filter_matches());
-          if (filter_iter == end)
-            continue;
+    // Find CellFilter in the tally (if any) to determine cells to loop over
+    const auto* cell_filter = tally.get_filter<CellFilter>();
+    const auto& cells =
+      cell_filter ? cell_filter->cells() : model::pulse_height_cells;
+
+    for (auto cell_id : cells) {
+      // Temporarily change cell of particle
+      p.n_coord() = 1;
+      p.coord(0).cell() = cell_id;
+
+      // Determine index of cell in model::pulse_height_cells
+      auto it = std::find(model::pulse_height_cells.begin(),
+        model::pulse_height_cells.end(), cell_id);
+      int index = std::distance(model::pulse_height_cells.begin(), it);
+
+      // Temporarily change energy of particle to pulse-height value
+      p.E_last() = p.pht_storage()[index];
+
+      // Initialize an iterator over valid filter bin combinations. If
+      // there are no valid combinations, use a continue statement to ensure
+      // we skip the assume_separate break below.
+      auto filter_iter = FilterBinIter(tally, p);
+      auto end = FilterBinIter(tally, true, &p.filter_matches());
+      if (filter_iter == end)
+        continue;
 
-          // Loop over filter bins.
-          for (; filter_iter != end; ++filter_iter) {
-            auto filter_index = filter_iter.index_;
-            auto filter_weight = filter_iter.weight_;
+      // Loop over filter bins.
+      for (; filter_iter != end; ++filter_iter) {
+        auto filter_index = filter_iter.index_;
+        auto filter_weight = filter_iter.weight_;
 
-            // Loop over scores.
-            for (auto score_index = 0; score_index < tally.scores_.size();
-                 ++score_index) {
+        // Loop over scores.
+        for (auto score_index = 0; score_index < tally.scores_.size();
+             ++score_index) {
 #pragma omp atomic
-              tally.results_(filter_index, score_index, TallyResult::VALUE) +=
-                filter_weight;
-            }
-          }
-
-          // Reset all the filter matches for the next tally event.
-          for (auto& match : p.filter_matches())
-            match.bins_present_ = false;
+          tally.results_(filter_index, score_index, TallyResult::VALUE) +=
+            filter_weight;
         }
       }
+
+      // Reset all the filter matches for the next tally event.
+      for (auto& match : p.filter_matches())
+        match.bins_present_ = false;
     }
     // Restore cell/energy
     p.n_coord() = orig_n_coord;