From f18caa1712a2a511ab4923be44c67bc60574d8e6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 16 Jan 2022 15:13:36 -0600 Subject: [PATCH] Get rid of wgt_absorb in ParticleData --- include/openmc/particle_data.h | 3 --- src/physics.cpp | 10 +++---- src/physics_mg.cpp | 8 +++--- src/tallies/tally_scoring.cpp | 49 +++++++++++++++++++++------------- 4 files changed, 37 insertions(+), 33 deletions(-) diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 38c0aa1df..a1804f7be 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -244,7 +244,6 @@ private: Position r_last_; //!< previous coordinates Direction u_last_; //!< previous direction coordinates double wgt_last_ {1.0}; //!< pre-collision particle weight - double wgt_absorb_ {0.0}; //!< weight absorbed for survival biasing // What event took place bool fission_ {false}; //!< did particle cause implicit fission @@ -375,8 +374,6 @@ public: const Position& u_last() const { return u_last_; } double& wgt_last() { return wgt_last_; } const double& wgt_last() const { return wgt_last_; } - double& wgt_absorb() { return wgt_absorb_; } - const double& wgt_absorb() const { return wgt_absorb_; } bool& fission() { return fission_; } TallyEvent& event() { return event_; } diff --git a/src/physics.cpp b/src/physics.cpp index f6fc31cb9..345b79947 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -134,8 +134,6 @@ void sample_neutron_reaction(Particle& p) if (p.neutron_xs(i_nuclide).absorption > 0.0) { absorption(p, i_nuclide); - } else { - p.wgt_absorb() = 0.0; } if (!p.alive()) return; @@ -626,15 +624,15 @@ void absorption(Particle& p, int i_nuclide) { if (settings::survival_biasing) { // Determine weight absorbed in survival biasing - p.wgt_absorb() = p.wgt() * p.neutron_xs(i_nuclide).absorption / - p.neutron_xs(i_nuclide).total; + double wgt_absorb = p.wgt() * p.neutron_xs(i_nuclide).absorption / + p.neutron_xs(i_nuclide).total; // Adjust weight of particle by probability of absorption - p.wgt() -= p.wgt_absorb(); + p.wgt() -= wgt_absorb; // Score implicit absorption estimate of keff if (settings::run_mode == RunMode::EIGENVALUE) { - p.keff_tally_absorption() += p.wgt_absorb() * + p.keff_tally_absorption() += wgt_absorb * p.neutron_xs(i_nuclide).nu_fission / p.neutron_xs(i_nuclide).absorption; } diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index ce7db3faf..d363c2d8e 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -55,8 +55,6 @@ void sample_reaction(Particle& p) // weight of the particle. Otherwise, it checks to see if absorption occurs. if (p.macro_xs().absorption > 0.) { absorption(p); - } else { - p.wgt_absorb() = 0.; } if (!p.alive()) return; @@ -216,14 +214,14 @@ void absorption(Particle& p) { if (settings::survival_biasing) { // Determine weight absorbed in survival biasing - p.wgt_absorb() = p.wgt() * p.macro_xs().absorption / p.macro_xs().total; + double wgt_absorb = p.wgt() * p.macro_xs().absorption / p.macro_xs().total; // Adjust weight of particle by the probability of absorption - p.wgt() -= p.wgt_absorb(); + p.wgt() -= wgt_absorb; // Score implicit absorpion estimate of keff p.keff_tally_absorption() += - p.wgt_absorb() * p.macro_xs().nu_fission / p.macro_xs().absorption; + wgt_absorb * p.macro_xs().nu_fission / p.macro_xs().absorption; } else { if (p.macro_xs().absorption > prn(p.current_seed()) * p.macro_xs().total) { p.keff_tally_absorption() += diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 8c91314aa..bd516ca18 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -546,6 +546,14 @@ void score_general_ce(Particle& p, int i_tally, int start_index, // Get the pre-collision energy of the particle. auto E = p.E_last(); + // Determine how much weight was absorbed due to survival biasing + double wgt_absorb = 0.0; + if (tally.estimator_ == TallyEstimator::ANALOG && + settings::survival_biasing) { + wgt_absorb = p.wgt_last() * p.neutron_xs(p.event_nuclide()).absorption / + p.neutron_xs(p.event_nuclide()).total; + } + using Type = ParticleType; for (auto i = 0; i < tally.scores_.size(); ++i) { @@ -614,7 +622,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index, continue; // Since only scattering events make it here, again we can use the // weight entering the collision as the estimator for the reaction rate - score = (p.wgt_last() - p.wgt_absorb()) * flux; + score = (p.wgt_last() - wgt_absorb) * flux; } else { if (i_nuclide >= 0) { if (p.type() == Type::neutron) { @@ -645,7 +653,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index, // For scattering production, we need to use the pre-collision weight // times the yield as the estimate for the number of neutrons exiting a // reaction with neutrons in the exit channel - score = (p.wgt_last() - p.wgt_absorb()) * flux; + score = (p.wgt_last() - wgt_absorb) * flux; // Don't waste time on very common reactions we know have multiplicities // of one. @@ -667,7 +675,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index, if (settings::survival_biasing) { // No absorption events actually occur if survival biasing is on -- // just use weight absorbed in survival biasing - score = p.wgt_absorb() * flux; + score = wgt_absorb * flux; } else { // Skip any event where the particle wasn't absorbed if (p.event() == TallyEvent::SCATTER) @@ -1240,7 +1248,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index, // Check if event MT matches if (p.event_mt() != ELASTIC) continue; - score = (p.wgt_last() - p.wgt_absorb()) * flux; + score = (p.wgt_last() - wgt_absorb) * flux; } else { if (i_nuclide >= 0) { if (p.neutron_xs(i_nuclide).elastic == CACHE_INVALID) @@ -1289,7 +1297,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index, // Check if the event MT matches if (p.event_mt() != score_bin) continue; - score = (p.wgt_last() - p.wgt_absorb()) * flux; + score = (p.wgt_last() - wgt_absorb) * flux; } else { int m; switch (score_bin) { @@ -1405,7 +1413,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index, // to check if it matches the MT number of the event if (p.event_mt() != score_bin) continue; - score = (p.wgt_last() - p.wgt_absorb()) * flux; + score = (p.wgt_last() - wgt_absorb) * flux; } else { // Any other cross section has to be calculated on-the-fly if (score_bin < 2) @@ -1452,10 +1460,13 @@ void score_general_mg(Particle& p, int i_tally, int start_index, // Set the direction and group to use with get_xs Direction p_u; int p_g; + double wgt_absorb = 0.0; if (tally.estimator_ == TallyEstimator::ANALOG || tally.estimator_ == TallyEstimator::COLLISION) { - if (settings::survival_biasing) { + // Determine weight that was absorbed + wgt_absorb = p.wgt_last() * p.neutron_xs(p.event_nuclide()).absorption / + p.neutron_xs(p.event_nuclide()).total; // Then we either are alive and had a scatter (and so g changed), // or are dead and g did not change @@ -1563,7 +1574,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index, continue; // Since only scattering events make it here, again we can use the // weight entering the collision as the estimator for the reaction rate - score = (p.wgt_last() - p.wgt_absorb()) * flux; + score = (p.wgt_last() - wgt_absorb) * flux; if (i_nuclide >= 0) { score *= atom_density * nuc_xs.get_xs(MgxsType::SCATTER_FMU, p.g_last(), &p.g(), @@ -1591,7 +1602,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index, // For scattering production, we need to use the pre-collision weight // times the multiplicity as the estimate for the number of neutrons // exiting a reaction with neutrons in the exit channel - score = (p.wgt_last() - p.wgt_absorb()) * flux; + score = (p.wgt_last() - wgt_absorb) * flux; // Since we transport based on material data, the angle selected // was not selected from the f(mu) for the nuclide. Therefore // adjust the score by the actual probability for that nuclide. @@ -1617,7 +1628,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index, if (settings::survival_biasing) { // No absorption events actually occur if survival biasing is on -- // just use weight absorbed in survival biasing - score = p.wgt_absorb() * flux; + score = wgt_absorb * flux; } else { // Skip any event where the particle wasn't absorbed if (p.event() == TallyEvent::SCATTER) @@ -1646,7 +1657,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // fission - score = p.wgt_absorb() * flux; + score = wgt_absorb * flux; } else { // Skip any non-absorption events if (p.event() == TallyEvent::SCATTER) @@ -1686,7 +1697,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // nu-fission - score = p.wgt_absorb() * flux; + score = wgt_absorb * flux; if (i_nuclide >= 0) { score *= atom_density * nuc_xs.get_xs(MgxsType::NU_FISSION, p_g) / macro_xs.get_xs(MgxsType::ABSORPTION, p_g); @@ -1733,7 +1744,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // prompt-nu-fission - score = p.wgt_absorb() * flux; + score = wgt_absorb * flux; if (i_nuclide >= 0) { score *= atom_density * nuc_xs.get_xs(MgxsType::PROMPT_NU_FISSION, p_g) / @@ -1794,7 +1805,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index, // Tally each delayed group bin individually for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { auto d = filt.groups()[d_bin] - 1; - score = p.wgt_absorb() * flux; + score = wgt_absorb * flux; if (i_nuclide >= 0) { score *= nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d) / @@ -1812,7 +1823,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index, // If the delayed group filter is not present, compute the score // by multiplying the absorbed weight by the fraction of the // delayed-nu-fission xs to the absorption xs - score = p.wgt_absorb() * flux; + score = wgt_absorb * flux; if (i_nuclide >= 0) { score *= nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g) / abs_xs; @@ -1909,7 +1920,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index, // Tally each delayed group bin individually for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { auto d = filt.groups()[d_bin] - 1; - score = p.wgt_absorb() * flux; + score = wgt_absorb * flux; if (i_nuclide >= 0) { score *= nuc_xs.get_xs( MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) * @@ -1935,14 +1946,14 @@ void score_general_mg(Particle& p, int i_tally, int start_index, score = 0.; for (auto d = 0; d < data::mg.num_delayed_groups_; ++d) { if (i_nuclide >= 0) { - score += p.wgt_absorb() * flux * + score += wgt_absorb * flux * nuc_xs.get_xs( MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) * nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d) / abs_xs; } else { - score += p.wgt_absorb() * flux * + score += wgt_absorb * flux * macro_xs.get_xs( MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) * macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, @@ -2050,7 +2061,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // fission scaled by the Q-value - score = p.wgt_absorb() * flux; + score = wgt_absorb * flux; } else { // Skip any non-absorption events if (p.event() == TallyEvent::SCATTER)