diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 0f37719b94f..de82a7d319b 100644 --- a/include/openmc/particle.h +++ b/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(); diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index fdacfa765b3..c93d798a0a9 100644 --- a/include/openmc/particle_data.h +++ b/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; diff --git a/src/initialize.cpp b/src/initialize.cpp index e2a5b974338..4bda392f779 100644 --- a/src/initialize.cpp +++ b/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; diff --git a/src/particle.cpp b/src/particle.cpp index 2d70d715e22..0b997138689 100644 --- a/src/particle.cpp +++ b/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; diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index 9333800bc36..04c0f2fc637 100644 --- a/src/weight_windows.cpp +++ b/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