Split storage optimization

include/openmc/particle.h

@@ -50,11 +50,11 @@ class Particle : public ParticleData {
   //! \return Whether a secondary particle was created
   bool create_secondary(double wgt, Direction u, double E, ParticleType type);
 
-  //! split a particle
+  //! replicate a particle
   //
-  //! creates a new particle with weight wgt
-  //! \param wgt Weight of the new particle
-  void split(double wgt);
+  //! creates replicated particles from this one
+  //! \param n_repl Number of replicas
+  void replicate(int n_repl);
 
   //! initialize from a source site
   //
@@ -64,6 +64,15 @@ class Particle : public ParticleData {
   //! \param src Source site data
   void from_source(const SourceSite* src);
 
+  //! initialize from a spendable source site
+  //
+  //! takes a particle from source site data and notifies that the site has
+  //! become spent.
+  //! \param src Source site data
+  //! \param empty_bank If false, the site does not contain remianed particle
+  //! replicas anymore
+  void from_spendable_source(SourceSite* src, bool& empty_bank);
+
   // Coarse-grained particle events
   void event_calculate_xs();
   void event_advance();

include/openmc/particle_data.h

@@ -46,6 +46,7 @@ struct SourceSite {
   double E;
   double time {0.0};
   double wgt {1.0};
+  int n_repl {1};
   int delayed_group {0};
   int surf_id {SURFACE_NONE};
   ParticleType particle;

src/initialize.cpp

@@ -157,7 +157,7 @@ void initialize_mpi(MPI_Comm intracomm)
 
   // Create bank datatype
   SourceSite b;
-  MPI_Aint disp[11];
+  MPI_Aint disp[12];
   MPI_Get_address(&b.r, &disp[0]);
   MPI_Get_address(&b.u, &disp[1]);
   MPI_Get_address(&b.E, &disp[2]);
@@ -169,14 +169,16 @@ void initialize_mpi(MPI_Comm intracomm)
   MPI_Get_address(&b.parent_nuclide, &disp[8]);
   MPI_Get_address(&b.parent_id, &disp[9]);
   MPI_Get_address(&b.progeny_id, &disp[10]);
-  for (int i = 10; i >= 0; --i) {
+  MPI_Get_address(&b.n_repl, &disp[11]);
+  for (int i = 11; i >= 0; --i) {
     disp[i] -= disp[0];
   }
 
-  int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1};
+  int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
   MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE,
-    MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_LONG, MPI_LONG};
-  MPI_Type_create_struct(11, blocks, disp, types, &mpi::source_site);
+    MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_LONG, MPI_LONG,
+    MPI_INT};
+  MPI_Type_create_struct(12, blocks, disp, types, &mpi::source_site);
   MPI_Type_commit(&mpi::source_site);
 
   CollisionTrackSite bc;

src/particle.cpp

@@ -92,16 +92,22 @@ bool Particle::create_secondary(
   return true;
 }
 
-void Particle::split(double wgt)
+void Particle::replicate(int n_repl)
 {
+  // Check for valid number of replicas
+  if (n_repl < 1)
+    return;
+
   auto& bank = secondary_bank().emplace_back();
   bank.particle = type();
-  bank.wgt = wgt;
+  bank.wgt = wgt();
   bank.r = r();
   bank.u = u();
   bank.E = settings::run_CE ? E() : g();
   bank.time = time();
 
+  bank.n_repl = n_repl;
+
   // Convert signed index to a signed surface ID
   if (surface() == SURFACE_NONE) {
     bank.surf_id = SURFACE_NONE;
@@ -157,6 +163,17 @@ void Particle::from_source(const SourceSite* src)
   }
 }
 
+void Particle::from_spendable_source(SourceSite* src, bool& empty_bank)
+{
+  // Pass data from source
+  from_source(src);
+
+  // Reduce the number of remaining particle replicas and determine if all the
+  // replicas are spent
+  src->n_repl--;
+  empty_bank = (src->n_repl == 0);
+}
+
 void Particle::event_calculate_xs()
 {
   // Set the random number stream
@@ -435,8 +452,10 @@ void Particle::event_revive_from_secondary()
     if (secondary_bank().empty())
       return;
 
-    from_source(&secondary_bank().back());
-    secondary_bank().pop_back();
+    bool spent;
+    from_spendable_source(&secondary_bank().back(), spent);
+    if (spent)
+      secondary_bank().pop_back();
     n_event() = 0;
     bank_second_E() = 0.0;
 

src/weight_windows.cpp

@@ -128,14 +128,10 @@ void apply_weight_windows(Particle& p)
     n_split = std::min(n_split, max_split);
 
     p.n_split() += n_split;
-
-    // Create secondaries and divide weight among all particles
-    int i_split = std::round(n_split);
-    for (int l = 0; l < i_split - 1; l++) {
-      p.split(weight / n_split);
-    }
     // remaining weight is applied to current particle
     p.wgt() = weight / n_split;
+    // Create secondaries from current particle
+    p.replicate(std::round(n_split) - 1);
 
   } else if (weight <= weight_window.lower_weight) {
     // if the particle weight is below the window, play Russian roulette