Tally 32-bit Overflow Fix

include/openmc/random_ray/source_region.h

@@ -51,10 +51,10 @@ inline void hash_combine(size_t& seed, const size_t v)
 // every iteration.
 struct TallyTask {
   int tally_idx;
-  int filter_idx;
+  int64_t filter_idx;
   int score_idx;
   int score_type;
-  TallyTask(int tally_idx, int filter_idx, int score_idx, int score_type)
+  TallyTask(int tally_idx, int64_t filter_idx, int score_idx, int score_type)
     : tally_idx(tally_idx), filter_idx(filter_idx), score_idx(score_idx),
       score_type(score_type)
   {}
@@ -690,7 +690,7 @@ class SourceRegionContainer {
   // Private Methods
 
   // Helper function for indexing
-  inline int index(int64_t sr, int g) const { return sr * negroups_ + g; }
+  inline int64_t index(int64_t sr, int g) const { return sr * negroups_ + g; }
 };
 
 } // namespace openmc

include/openmc/tallies/tally.h

@@ -97,9 +97,9 @@ class Tally {
   //! Given already-set filters, set the stride lengths
   void set_strides();
 
-  int32_t strides(int i) const { return strides_[i]; }
+  int64_t strides(int i) const { return strides_[i]; }
 
-  int32_t n_filter_bins() const { return n_filter_bins_; }
+  int64_t n_filter_bins() const { return n_filter_bins_; }
 
   bool multiply_density() const { return multiply_density_; }
 
@@ -184,9 +184,9 @@ class Tally {
   vector<int32_t> filters_; //!< Filter indices in global filters array
 
   //! Index strides assigned to each filter to support 1D indexing.
-  vector<int32_t> strides_;
+  vector<int64_t> strides_;
 
-  int32_t n_filter_bins_ {0};
+  int64_t n_filter_bins_ {0};
 
   //! Whether to multiply by atom density for reaction rates
   bool multiply_density_ {true};

include/openmc/tallies/tally_scoring.h

@@ -41,7 +41,7 @@ class FilterBinIter {
 
   FilterBinIter& operator++();
 
-  int index_ {1};
+  int64_t index_ {1};
   double weight_ {1.};
 
   vector<FilterMatch>& filter_matches_;

src/mesh.cpp

@@ -1,8 +1,9 @@
 #include "openmc/mesh.h"
 #include <algorithm> // for copy, equal, min, min_element
 #include <cassert>
 #include <cstdint>        // for uint64_t
 #include <cstring>        // for memcpy
+#include <limits>
 #define _USE_MATH_DEFINES // to make M_PI declared in Intel and MSVC compilers
 #include <cmath>          // for ceil
 #include <cstddef>        // for size_t
@@ -1080,13 +1081,26 @@
 
 int StructuredMesh::n_bins() const
 {
-  return std::accumulate(
-    shape_.begin(), shape_.begin() + n_dimension_, 1, std::multiplies<>());
+  // Bin indices are stored as 32-bit ints in the tally system.
+  int64_t n = 1;
+  for (int i = 0; i < n_dimension_; ++i)
+    n *= shape_[i];
+  if (n > std::numeric_limits<int>::max()) {
+    fatal_error(fmt::format(
+      "Mesh {} has too many bins ({}) for 32-bit tally indexing", id_, n));
+  }
+  return static_cast<int>(n);
 }
 
 int StructuredMesh::n_surface_bins() const
 {
-  return 4 * n_dimension_ * n_bins();
+  // Surface bin indices are stored as 32-bit ints in the tally system.
+  int64_t n = static_cast<int64_t>(n_bins()) * 4 * n_dimension_;
+  if (n > std::numeric_limits<int>::max()) {
+    fatal_error(fmt::format(
+      "Mesh {} has too many surface bins ({}) for tally indexing", id_, n));
+  }
+  return static_cast<int>(n);
 }
 
 tensor::Tensor<double> StructuredMesh::count_sites(

src/random_ray/flat_source_domain.cpp

@@ -724,7 +724,7 @@ void FlatSourceDomain::random_ray_tally()
     for (int i = 0; i < model::tallies.size(); i++) {
       Tally& tally {*model::tallies[i]};
 #pragma omp parallel for
-      for (int bin = 0; bin < tally.n_filter_bins(); bin++) {
+      for (int64_t bin = 0; bin < tally.n_filter_bins(); bin++) {
         for (int score_idx = 0; score_idx < tally.n_scores(); score_idx++) {
           auto score_type = tally.scores_[score_idx];
           if (score_type == SCORE_FLUX) {

src/tallies/tally.cpp

@@ -512,7 +512,7 @@ void Tally::set_strides()
   // longest stride.
   auto n = filters_.size();
   strides_.resize(n, 0);
-  int stride = 1;
+  int64_t stride = 1;
   for (int i = n - 1; i >= 0; --i) {
     strides_[i] = stride;
     stride *= model::tally_filters[filters_[i]]->n_bins();

src/tallies/tally_scoring.cpp

@@ -158,7 +158,7 @@ void score_fission_delayed_dg(int i_tally, int d_bin, double score,
   dg_match.bins_[i_bin] = d_bin;
 
   // Determine the filter scoring index
-  auto filter_index = 0;
+  int64_t filter_index = 0;
   double filter_weight = 1.;
   for (auto i = 0; i < tally.filters().size(); ++i) {
     auto i_filt = tally.filters(i);
@@ -449,7 +449,7 @@ void score_fission_eout(Particle& p, int i_tally, int i_score, int score_bin)
         (score_bin == SCORE_PROMPT_NU_FISSION && g == 0)) {
 
       // Find the filter scoring index for this filter combination
-      int filter_index = 0;
+      int64_t filter_index = 0;
       double filter_weight = 1.0;
       for (auto j = 0; j < tally.filters().size(); ++j) {
         auto i_filt = tally.filters(j);
@@ -497,7 +497,7 @@ void score_fission_eout(Particle& p, int i_tally, int i_score, int score_bin)
       } else {
 
         // Find the filter index and weight for this filter combination
-        int filter_index = 0;
+        int64_t filter_index = 0;
         double filter_weight = 1.;
         for (auto j = 0; j < tally.filters().size(); ++j) {
           auto i_filt = tally.filters(j);
@@ -578,8 +578,8 @@ double get_nuclide_xs(const Particle& p, int i_nuclide, int score_bin)
 //! collision estimator.
 
 void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
-  int filter_index, double filter_weight, int i_nuclide, double atom_density,
-  double flux)
+  int64_t filter_index, double filter_weight, int i_nuclide,
+  double atom_density, double flux)
 {
   Tally& tally {*model::tallies[i_tally]};
 
@@ -1112,8 +1112,8 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index,
 //! is not used for analog tallies.
 
 void score_general_ce_analog(Particle& p, int i_tally, int start_index,
-  int filter_index, double filter_weight, int i_nuclide, double atom_density,
-  double flux)
+  int64_t filter_index, double filter_weight, int i_nuclide,
+  double atom_density, double flux)
 {
   Tally& tally {*model::tallies[i_tally]};
 
@@ -1615,8 +1615,8 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index,
 //! argument is really just used for filter weights.
 
 void score_general_mg(Particle& p, int i_tally, int start_index,
-  int filter_index, double filter_weight, int i_nuclide, double atom_density,
-  double flux)
+  int64_t filter_index, double filter_weight, int i_nuclide,
+  double atom_density, double flux)
 {
   auto& tally {*model::tallies[i_tally]};
 

tests/cpp_unit_tests/test_tally.cpp

@@ -1,3 +1,4 @@
+#include "openmc/tallies/filter_energy.h"
 #include "openmc/tallies/tally.h"
 #include <catch2/catch_test_macros.hpp>
 
@@ -51,5 +52,28 @@
   tally->set_filters(filters);
   REQUIRE(tally->filters().size() == 1);
   REQUIRE(model::filter_map[cell_filter->id()] == tally->filters(0));
+}
 
+// Regression test for 64-bit tally filter-bin counts (mesh x groups > 2^31).
+TEST_CASE("Tally filter-bin count does not overflow 32 bits")
+{
+  // Two energy filters whose bin counts multiply to 2.5e9, above INT32_MAX.
+  constexpr int64_t bins_per_filter = 50000;
+
+  // Only the bin count matters here, so the edge values are an arbitrary ramp.
+  std::vector<double> edges(bins_per_filter + 1);
+  for (int64_t i = 0; i < bins_per_filter + 1; ++i)
+    edges[i] = static_cast<double>(i);
+
+  Tally* tally = Tally::create();
+  for (int i = 0; i < 2; ++i) {
+    Filter* filter = Filter::create("energy");
+    dynamic_cast<EnergyFilter*>(filter)->set_bins(edges);
+    tally->add_filter(filter);
+  }
+  tally->set_strides();
+
+  // set_strides() previously accumulated this product in a 32-bit int.
+  REQUIRE(tally->n_filter_bins() == bins_per_filter * bins_per_filter);
+  REQUIRE(tally->n_filter_bins() > 2147483647);
 }