Tally 32-bit Overflow Fix (#3960)

Co-authored-by: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

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 @@ private:
   // 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 @@ public:
   //! 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 @@ private:
   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 @@ public:
 
   FilterBinIter& operator++();
 
-  int index_ {1};
+  int64_t index_ {1};
   double weight_ {1.};
 
   vector<FilterMatch>& filter_matches_;

src/mesh.cpp

@@ -6,7 +6,8 @@
 #define _USE_MATH_DEFINES // to make M_PI declared in Intel and MSVC compilers
 #include <cmath>          // for ceil
 #include <cstddef>        // for size_t
-#include <numeric>        // for accumulate
+#include <limits>
+#include <numeric> // for accumulate
 #include <string>
 
 #ifdef _MSC_VER
@@ -1080,13 +1081,26 @@ int StructuredMesh::get_bin(Position r) const
 
 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

@@ -725,7 +725,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();
@@ -887,7 +887,7 @@ void Tally::accumulate()
     if (higher_moments_) {
 #pragma omp parallel for
       // filter bins (specific cell, energy bins)
-      for (int i = 0; i < results_.shape(0); ++i) {
+      for (int64_t i = 0; i < results_.shape(0); ++i) {
         // score bins (flux, total reaction rate, fission reaction rate, etc.)
         for (int j = 0; j < results_.shape(1); ++j) {
           double val = results_(i, j, TallyResult::VALUE) * norm;
@@ -902,7 +902,7 @@ void Tally::accumulate()
     } else {
 #pragma omp parallel for
       // filter bins (specific cell, energy bins)
-      for (int i = 0; i < results_.shape(0); ++i) {
+      for (int64_t i = 0; i < results_.shape(0); ++i) {
         // score bins (flux, total reaction rate, fission reaction rate, etc.)
         for (int j = 0; j < results_.shape(1); ++j) {
           double val = results_(i, j, TallyResult::VALUE) * norm;

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]};
 

src/tallies/trigger.cpp

@@ -28,7 +28,7 @@ KTrigger keff_trigger;
 //==============================================================================
 
 std::pair<double, double> get_tally_uncertainty(
-  int i_tally, int score_index, int filter_index)
+  int i_tally, int score_index, int64_t filter_index)
 {
   const auto& tally {model::tallies[i_tally]};
 
@@ -71,7 +71,7 @@ void check_tally_triggers(double& ratio, int& tally_id, int& score)
         continue;
 
       const auto& results = t.results_;
-      for (auto filter_index = 0; filter_index < results.shape(0);
+      for (int64_t filter_index = 0; filter_index < results.shape(0);
            ++filter_index) {
         // Compute the tally uncertainty metrics.
         auto uncert_pair =

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>
 
@@ -46,10 +47,34 @@ TEST_CASE("Test add/set_filter")
   REQUIRE(model::filter_map[cell_filter->id()] == tally->filters(0));
   REQUIRE(model::filter_map[particle_filter->id()] == tally->filters(1));
 
-  // set filters with a duplicate filter, should only add the filter to the tally once
+  // set filters with a duplicate filter, should only add the filter to the
+  // tally once
   filters = {cell_filter, cell_filter};
   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);
 }