diff --git a/docs/source/methods/neutron_physics.rst b/docs/source/methods/neutron_physics.rst index 793e852b395..70a57598ad9 100644 --- a/docs/source/methods/neutron_physics.rst +++ b/docs/source/methods/neutron_physics.rst @@ -1688,12 +1688,12 @@ This parameter can be used to make the adjustable parameters of survival biasing It is currently only implemented for phase-space sources (MCPL file sources and HDF5 surface sources). For these, the normalization of the parameters -is done per source particle. That is for each history the parameters, +is done per source particle. That is, for each history, the parameters :math:`w_c` and :math:`w_s`, are multiplied by the start weight of the current history. This normalization is recommended for problems where numerous source particles are -initialized below the survival biasing weight cutoff, :math:`w_c`; to prevent them +initialized below the survival biasing weight cutoff :math:`w_c`; to prevent them from being immediately rouletted. This can notably be the case for biased sources. Weight Windows diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 430786d3cc5..d2af4b47d9d 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -264,8 +264,7 @@ class ParticleData { // Other physical data double wgt_ {1.0}; //!< particle weight - double wgt_cutoff_; //!< particle weight cutoff - double wgt_survive_; //!< particle weight survive + double wgt0_ {1.0}; //!< particle start weight double mu_; //!< angle of scatter double time_ {0.0}; //!< time in [s] double time_last_ {0.0}; //!< previous time in [s] @@ -399,12 +398,9 @@ class ParticleData { double& wgt() { return wgt_; } double wgt() const { return wgt_; } - double& wgt_cutoff() { return wgt_cutoff_; } - double wgt_cutoff() const { return wgt_cutoff_; } - void wgt_cutoff(double cutoff) { wgt_cutoff_ = cutoff; } - double& wgt_survive() { return wgt_survive_; } - double wgt_survive() const { return wgt_survive_; } - void wgt_survive(double survive) { wgt_survive_ = survive; } + double& wgt0() { return wgt0_; } + double wgt0() const { return wgt0_; } + double& mu() { return mu_; } const double& mu() const { return mu_; } double& time() { return time_; } diff --git a/src/physics.cpp b/src/physics.cpp index 7e142d2fb5d..db5ec2f56b5 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -157,13 +157,13 @@ void sample_neutron_reaction(Particle& p) // Play russian roulette if survival biasing is turned on if (settings::survival_biasing) { - double weight_cutoff = settings::weight_cutoff; - // if survival normalization is applicable, use normalized weight cutoff - if((settings::source_file || settings::surf_source_read)&&(settings::survival_normalization)){ - weight_cutoff = p.wgt_cutoff(); - } - if (p.wgt() < weight_cutoff) { - russian_roulette(p, settings::weight_survive); + // if survival normalization is applicable, use normalized weight cutoff and normalized weight survive + if ((settings::source_file || settings::surf_source_read)&&(settings::survival_normalization)) { + if (p.wgt() < settings::weight_cutoff*p.wgt0()) { + russian_roulette(p, settings::weight_survive*p.wgt0()); + } + } else if (p.wgt() < settings::weight_cutoff) { + russian_roulette(p, settings::weight_survive); } } } diff --git a/src/physics_common.cpp b/src/physics_common.cpp index a1e768f2522..e25ae6b97ab 100644 --- a/src/physics_common.cpp +++ b/src/physics_common.cpp @@ -11,15 +11,11 @@ namespace openmc { void russian_roulette(Particle& p, double weight_survive) { - // if survival normalization is applicable, use normalized weight survive - if((settings::source_file || settings::surf_source_read)&&(settings::survival_normalization && settings::survival_biasing)){ - weight_survive = p.wgt_survive(); - } if (weight_survive * prn(p.current_seed()) < p.wgt()) { p.wgt() = weight_survive; } else { p.wgt() = 0.; - } + } } } // namespace openmc diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index ffa8c9e1390..43e3cfa7455 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -64,12 +64,7 @@ void sample_reaction(Particle& p) // Play Russian roulette if survival biasing is turned on if (settings::survival_biasing) { - double weight_cutoff = settings::weight_cutoff; - // if survival normalization is applicable, use normalized weight cutoff - if((settings::source_file || settings::surf_source_read)&&(settings::survival_normalization)){ - weight_cutoff = p.wgt_cutoff(); - } - if (p.wgt() < weight_cutoff) { + if (p.wgt() < settings::weight_cutoff) { russian_roulette(p, settings::weight_survive); } } diff --git a/src/simulation.cpp b/src/simulation.cpp index 1f619141506..bcf3d2795fc 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -533,19 +533,15 @@ void initialize_history(Particle& p, int64_t index_source) auto site = sample_external_source(&seed); p.from_source(&site); } - // normalize biasing weight parameters; multiply them by start weight of initialized history - // applicable only for MCPL and HDF5 phase-space sources. - if(settings::source_file || settings::surf_source_read){ - if(settings::survival_normalization && settings::survival_biasing){ - p.wgt_cutoff(settings::weight_cutoff * p.wgt()); - p.wgt_survive(settings::weight_survive * p.wgt()); - } - } + p.current_work() = index_source; // set identifier for particle p.id() = simulation::work_index[mpi::rank] + index_source; + // set particle history start weight + p.wgt0(p.wgt()); + // set progeny count to zero p.n_progeny() = 0; @@ -582,7 +578,7 @@ void initialize_history(Particle& p, int64_t index_source) write_message("Simulating Particle {}", p.id()); } -// Add paricle's starting weight to count for normalizing tallies later +// Add particle's starting weight to count for normalizing tallies later #pragma omp atomic simulation::total_weight += p.wgt();