Make secondary_bank_ a vector and improve interfaces

include/openmc/bank.h

@@ -20,15 +20,13 @@ namespace openmc {
 
 namespace simulation {
 
-extern "C" int64_t n_bank;
-
 extern std::vector<Particle::Bank> source_bank;
 extern std::vector<Particle::Bank> fission_bank;
 #ifdef _OPENMP
 extern std::vector<Particle::Bank> master_fission_bank;
 #endif
 
-#pragma omp threadprivate(fission_bank, n_bank)
+#pragma omp threadprivate(fission_bank)
 
 } // namespace simulation
 

include/openmc/mesh.h

@@ -96,14 +96,11 @@ class RegularMesh {
 
   //! Count number of bank sites in each mesh bin / energy bin
   //
-  //! \param[in] n Number of bank sites
   //! \param[in] bank Array of bank sites
-  //! \param[in] n_energy Number of energies
-  //! \param[in] energies Array of energies
   //! \param[out] Whether any bank sites are outside the mesh
   //! \return Array indicating number of sites in each mesh/energy bin
-  xt::xarray<double> count_sites(int64_t n, const Particle::Bank* bank,
-    int n_energy, const double* energies, bool* outside) const;
+  xt::xarray<double> count_sites(const std::vector<Particle::Bank>& bank,
+    bool* outside) const;
 
   int id_ {-1};  //!< User-specified ID
   int n_dimension_; //!< Number of dimensions

include/openmc/particle.h

@@ -285,8 +285,7 @@ class Particle {
   bool write_track_ {false};
 
   // Secondary particles created
-  int64_t n_secondary_ {};
-  Bank secondary_bank_[MAX_SECONDARY];
+  std::vector<Bank> secondary_bank_;
 };
 
 } // namespace openmc

include/openmc/physics.h

@@ -7,6 +7,8 @@
 #include "openmc/position.h"
 #include "openmc/reaction.h"
 
+#include <vector>
+
 namespace openmc {
 
 //==============================================================================
@@ -47,7 +49,7 @@ int sample_nuclide(const Particle* p);
 //! Determine the average total, prompt, and delayed neutrons produced from
 //! fission and creates appropriate bank sites.
 void create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
-  Particle::Bank* bank_array, int64_t* bank_size, int64_t bank_capacity);
+  std::vector<Particle::Bank>& bank);
 
 int sample_element(Particle* p);
 

include/openmc/physics_mg.h

@@ -8,6 +8,8 @@
 #include "openmc/particle.h"
 #include "openmc/nuclide.h"
 
+#include <vector>
+
 namespace openmc {
 
 //! \brief samples particle behavior after a collision event.
@@ -31,12 +33,9 @@ scatter(Particle* p);
 //! \brief Determines the average total, prompt and delayed neutrons produced
 //! from fission and creates the appropriate bank sites.
 //! \param p Particle to operate on
-//! \param bank_array The particle bank to populate
-//! \param size_bank Number of particles currently in the bank
-//! \param bank_array_size Allocated size of the bank
+//! \param bank The particle bank to populate
 void
-create_fission_sites(Particle* p, Particle::Bank* bank_array, int64_t* size_bank,
-     int64_t bank_array_size);
+create_fission_sites(Particle* p, std::vector<Particle::Bank>& bank);
 
 //! \brief Handles an absorption event
 //! \param p Particle to operate on

src/bank.cpp

@@ -16,8 +16,6 @@ namespace openmc {
 
 namespace simulation {
 
-int64_t n_bank;
-
 std::vector<Particle::Bank> source_bank;
 std::vector<Particle::Bank> fission_bank;
 #ifdef _OPENMP

src/eigenvalue.cpp

@@ -23,6 +23,7 @@
 #include <algorithm> // for min
 #include <array>
 #include <cmath> // for sqrt, abs, pow
+#include <iterator> // for back_inserter
 #include <string>
 
 namespace openmc {
@@ -81,29 +82,30 @@ void synchronize_bank()
 
 #ifdef OPENMC_MPI
   int64_t start = 0;
-  MPI_Exscan(&simulation::n_bank, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm);
+  int64_t n_bank = simulation::fission_bank.size();
+  MPI_Exscan(&n_bank, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm);
 
   // While we would expect the value of start on rank 0 to be 0, the MPI
   // standard says that the receive buffer on rank 0 is undefined and not
   // significant
   if (mpi::rank == 0) start = 0;
 
-  int64_t finish = start + simulation::n_bank;
+  int64_t finish = start + simulation::fission_bank.size();
   int64_t total = finish;
   MPI_Bcast(&total, 1, MPI_INT64_T, mpi::n_procs - 1, mpi::intracomm);
 
 #else
   int64_t start  = 0;
-  int64_t finish = simulation::n_bank;
-  int64_t total  = simulation::n_bank;
+  int64_t finish = simulation::fission_bank.size();
+  int64_t total  = simulation::fission_bank.size();
 #endif
 
   // If there are not that many particles per generation, it's possible that no
   // fission sites were created at all on a single processor. Rather than add
   // extra logic to treat this circumstance, we really want to ensure the user
   // runs enough particles to avoid this in the first place.
 
-  if (simulation::n_bank == 0) {
+  if (simulation::fission_bank.empty()) {
     fatal_error("No fission sites banked on MPI rank " + std::to_string(mpi::rank));
   }
 
@@ -134,21 +136,21 @@ void synchronize_bank()
   int64_t index_temp = 0;
   std::vector<Particle::Bank> temp_sites(3*simulation::work);
 
-  for (int64_t i = 0; i < simulation::n_bank; ++i) {
+  for (const auto& site : simulation::fission_bank) {
     // If there are less than n_particles particles banked, automatically add
     // int(n_particles/total) sites to temp_sites. For example, if you need
     // 1000 and 300 were banked, this would add 3 source sites per banked site
     // and the remaining 100 would be randomly sampled.
     if (total < settings::n_particles) {
       for (int64_t j = 1; j <= settings::n_particles / total; ++j) {
-        temp_sites[index_temp] = simulation::fission_bank[i];
+        temp_sites[index_temp] = site;
         ++index_temp;
       }
     }
 
     // Randomly sample sites needed
     if (prn() < p_sample) {
-      temp_sites[index_temp] = simulation::fission_bank[i];
+      temp_sites[index_temp] = site;
       ++index_temp;
     }
   }
@@ -189,7 +191,7 @@ void synchronize_bank()
       // fission bank
       sites_needed = settings::n_particles - finish;
       for (int i = 0; i < sites_needed; ++i) {
-        int i_bank = simulation::n_bank - sites_needed + i;
+        int i_bank = simulation::fission_bank.size() - sites_needed + i;
         temp_sites[index_temp] = simulation::fission_bank[i_bank];
         ++index_temp;
       }
@@ -346,38 +348,30 @@ void calculate_average_keff()
 #ifdef _OPENMP
 void join_bank_from_threads()
 {
-  // Initialize the total number of fission bank sites
-  int64_t total = 0;
-
-#pragma omp parallel
+  #pragma omp parallel
   {
     // Copy thread fission bank sites to one shared copy
-#pragma omp for ordered schedule(static)
+    #pragma omp for ordered schedule(static)
     for (int i = 0; i < simulation::n_threads; ++i) {
-#pragma omp ordered
+      #pragma omp ordered
       {
         std::copy(
-          &simulation::fission_bank[0],
-          &simulation::fission_bank[0] + simulation::n_bank,
-          &simulation::master_fission_bank[total]
+          simulation::fission_bank.cbegin(),
+          simulation::fission_bank.cend(),
+          std::back_inserter(simulation::master_fission_bank)
         );
-        total += simulation::n_bank;
       }
     }
 
     // Make sure all threads have made it to this point
-#pragma omp barrier
+    #pragma omp barrier
 
     // Now copy the shared fission bank sites back to the master thread's copy.
     if (simulation::thread_id == 0) {
-      simulation::n_bank = total;
-      std::copy(
-        &simulation::master_fission_bank[0],
-        &simulation::master_fission_bank[0] + simulation::n_bank,
-        &simulation::fission_bank[0]
-      );
+      simulation::fission_bank = simulation::master_fission_bank;
+      simulation::master_fission_bank.clear();
     } else {
-      simulation::n_bank = 0;
+      simulation::fission_bank.clear();
     }
   }
 }
@@ -537,8 +531,8 @@ void shannon_entropy()
 
   // Get source weight in each mesh bin
   bool sites_outside;
-  xt::xtensor<double, 1> p = m->count_sites(simulation::n_bank,
-    simulation::fission_bank.data(), 0, nullptr, &sites_outside);
+  xt::xtensor<double, 1> p = m->count_sites(simulation::fission_bank,
+    &sites_outside);
 
   // display warning message if there were sites outside entropy box
   if (sites_outside) {
@@ -577,8 +571,8 @@ void ufs_count_sites()
   } else {
     // count number of source sites in each ufs mesh cell
     bool sites_outside;
-    simulation::source_frac = m->count_sites(simulation::work,
-      simulation::source_bank.data(), 0, nullptr, &sites_outside);
+    simulation::source_frac = m->count_sites(simulation::source_bank,
+      &sites_outside);
 
     // Check for sites outside of the mesh
     if (mpi::master && sites_outside) {

src/mesh.cpp

@@ -732,45 +732,30 @@ void RegularMesh::to_hdf5(hid_t group) const
   close_group(mesh_group);
 }
 
-xt::xarray<double> RegularMesh::count_sites(int64_t n, const Particle::Bank* bank,
-  int n_energy, const double* energies, bool* outside) const
+xt::xarray<double>
+RegularMesh::count_sites(const std::vector<Particle::Bank>& bank,
+  bool* outside) const
 {
   // Determine shape of array for counts
   std::size_t m = xt::prod(shape_)();
-  std::vector<std::size_t> shape;
-  if (n_energy > 0) {
-    shape = {m, static_cast<std::size_t>(n_energy - 1)};
-  } else {
-    shape = {m};
-  }
+  std::vector<std::size_t> shape = {m};
 
   // Create array of zeros
   xt::xarray<double> cnt {shape, 0.0};
   bool outside_ = false;
 
-  for (int64_t i = 0; i < n; ++i) {
+  for (const auto& site : bank) {
     // determine scoring bin for entropy mesh
-    int mesh_bin = get_bin(bank[i].r);
+    int mesh_bin = get_bin(site.r);
 
     // if outside mesh, skip particle
     if (mesh_bin < 0) {
       outside_ = true;
       continue;
     }
 
-    if (n_energy > 0) {
-      double E = bank[i].E;
-      if (E >= energies[0] && E <= energies[n_energy - 1]) {
-        // determine energy bin
-        int e_bin = lower_bound_index(energies, energies + n_energy, E);
-
-        // Add to appropriate bin
-        cnt(mesh_bin, e_bin) += bank[i].wgt;
-      }
-    } else {
-      // Add to appropriate bin
-      cnt(mesh_bin) += bank[i].wgt;
-    }
+    // Add to appropriate bin
+    cnt(mesh_bin) += site.wgt;
   }
 
   // Create copy of count data

src/particle.cpp

@@ -75,17 +75,14 @@ Particle::clear()
 void
 Particle::create_secondary(Direction u, double E, Type type)
 {
-  if (n_secondary_ == MAX_SECONDARY) {
-    fatal_error("Too many secondary particles created.");
-  }
-
-  int64_t n = n_secondary_;
-  secondary_bank_[n].particle = type;
-  secondary_bank_[n].wgt = wgt_;
-  secondary_bank_[n].r = this->r();
-  secondary_bank_[n].u = u;
-  secondary_bank_[n].E = settings::run_CE ? E : g_;
-  ++n_secondary_;
+  secondary_bank_.emplace_back();
+
+  auto& bank {secondary_bank_.back()};
+  bank.particle = type;
+  bank.wgt = wgt_;
+  bank.r = this->r();
+  bank.u = u;
+  bank.E = settings::run_CE ? E : g_;
 }
 
 void
@@ -360,10 +357,10 @@ Particle::transport()
     // Check for secondary particles if this particle is dead
     if (!alive_) {
       // If no secondary particles, break out of event loop
-      if (n_secondary_ == 0) break;
+      if (secondary_bank_.empty()) break;
 
-      this->from_source(&secondary_bank_[n_secondary_ - 1]);
-      --n_secondary_;
+      this->from_source(&secondary_bank_.back());
+      secondary_bank_.pop_back();
       n_event = 0;
 
       // Enter new particle in particle track file