crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 13:45:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge pull request #1974 from paulromano/survival-biasing-fix

Browse source 
Remove wgt_absorb data member from Particle
This commit is contained in:
Gavin Ridley 2022-02-18 10:31:01 -05:00 • committed by GitHub
commit 831c8d1c50
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
10 changed files with 142 additions and 196 deletions
Download patch file Download diff file Expand all files Collapse all files

3
include/openmc/particle_data.h
View file

@ -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_; }

4
include/openmc/physics_common.h
View file

@ -9,7 +9,9 @@
namespace openmc {
//! \brief Performs the russian roulette operation for a particle
void russian_roulette(Particle& p);
//! \param[in,out] p Particle object
//! \param[in] weight_survive Weight assigned to particles that survive
void russian_roulette(Particle& p, double weight_survive);
} // namespace openmc
#endif // OPENMC_PHYSICS_COMMON_H

2
openmc/surface.py
View file

@ -2633,7 +2633,7 @@ class Halfspace(Region):
memo = {}
# If rotated surface not in memo, add it
key = (self.surface, tuple(rotation), tuple(pivot), order, inplace)
key = (self.surface, tuple(np.ravel(rotation)), tuple(pivot), order, inplace)
if key not in memo:
memo[key] = self.surface.rotate(rotation, pivot=pivot, order=order,
inplace=inplace)

27
src/physics.cpp
View file

@ -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;
@ -153,9 +151,9 @@ void sample_neutron_reaction(Particle& p)
// Play russian roulette if survival biasing is turned on
if (settings::survival_biasing) {
russian_roulette(p);
if (!p.alive())
return;
if (p.wgt() < settings::weight_cutoff) {
russian_roulette(p, settings::weight_survive);
}
}
}
@ -626,16 +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;
const 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_last() = p.wgt();
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;
}
@ -952,9 +949,8 @@ Direction sample_target_velocity(const Nuclide& nuc, double E, Direction u,
// cdf value at upper bound attainable energy
double m = (nuc.xs_cdf_[i_E_up + 1] - nuc.xs_cdf_[i_E_up]) /
(nuc.energy_0K_[i_E_up + 1] - nuc.energy_0K_[i_E_up]);
double cdf_up =
nuc.xs_cdf_[i_E_up] + m * (E_up - nuc.energy_0K_[i_E_up]);
(nuc.energy_0K_[i_E_up + 1] - nuc.energy_0K_[i_E_up]);
double cdf_up = nuc.xs_cdf_[i_E_up] + m * (E_up - nuc.energy_0K_[i_E_up]);
while (true) {
// directly sample Maxwellian
@ -962,9 +958,8 @@ Direction sample_target_velocity(const Nuclide& nuc, double E, Direction u,
// sample a relative energy using the xs cdf
double cdf_rel = cdf_low + prn(seed) * (cdf_up - cdf_low);
int i_E_rel = lower_bound_index(nuc.xs_cdf_.begin() + i_E_low,
nuc.xs_cdf_.begin() + i_E_up+2,
cdf_rel);
int i_E_rel = lower_bound_index(nuc.xs_cdf_.begin() + i_E_low,
nuc.xs_cdf_.begin() + i_E_up + 2, cdf_rel);
double E_rel = nuc.energy_0K_[i_E_low + i_E_rel];
double m = (nuc.xs_cdf_[i_E_low + i_E_rel + 1] -
nuc.xs_cdf_[i_E_low + i_E_rel]) /

16
src/physics_common.cpp
View file

@ -9,17 +9,13 @@ namespace openmc {
// RUSSIAN_ROULETTE
//==============================================================================
void russian_roulette(Particle& p)
void russian_roulette(Particle& p, double weight_survive)
{
if (p.wgt() < settings::weight_cutoff) {
if (prn(p.current_seed()) < p.wgt() / settings::weight_survive) {
p.wgt() = settings::weight_survive;
p.wgt_last() = p.wgt();
} else {
p.wgt() = 0.;
p.wgt_last() = 0.;
p.alive() = false;
}
if (weight_survive * prn(p.current_seed()) < p.wgt()) {
p.wgt() = weight_survive;
} else {
p.wgt() = 0.;
p.alive() = false;
}
}

15
src/physics_mg.cpp
View file

@ -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;
@ -66,9 +64,9 @@ void sample_reaction(Particle& p)
// Play Russian roulette if survival biasing is turned on
if (settings::survival_biasing) {
russian_roulette(p);
if (!p.alive())
return;
if (p.wgt() < settings::weight_cutoff) {
russian_roulette(p, settings::weight_survive);
}
}
}
@ -216,15 +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_last() = p.wgt();
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() +=

246
src/tallies/tally_scoring.cpp
View file

@ -203,10 +203,10 @@ double score_fission_q(const Particle& p, int score_bin, const Tally& tally,
// 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
if (p.neutron_xs(p.event_nuclide()).absorption > 0) {
return p.wgt_absorb() * get_nuc_fission_q(nuc, p, score_bin) *
if (p.neutron_xs(p.event_nuclide()).total > 0) {
return p.wgt_last() * get_nuc_fission_q(nuc, p, score_bin) *
p.neutron_xs(p.event_nuclide()).fission * flux /
p.neutron_xs(p.event_nuclide()).absorption;
p.neutron_xs(p.event_nuclide()).total;
}
} else {
// Skip any non-absorption events
@ -291,7 +291,6 @@ double get_nuclide_neutron_heating(
double score_neutron_heating(const Particle& p, const Tally& tally, double flux,
int rxn_bin, int i_nuclide, double atom_density)
{
double score;
// Get heating macroscopic "cross section"
double heating_xs;
if (i_nuclide >= 0) {
@ -318,17 +317,12 @@ double score_neutron_heating(const Particle& p, const Tally& tally, double flux,
}
}
}
score = heating_xs * flux;
double score = heating_xs * flux;
if (tally.estimator_ == TallyEstimator::ANALOG) {
// All events score to a heating tally bin. We actually use a
// collision estimator in place of an analog one since there is no
// reaction-wise heating cross section
if (settings::survival_biasing) {
// Account for the fact that some weight has been absorbed
score *= p.wgt_last() + p.wgt_absorb();
} else {
score *= p.wgt_last();
}
score *= p.wgt_last();
}
return score;
}
@ -552,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) {
@ -565,16 +567,8 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
// All events score to a flux bin. We actually use a collision estimator
// in place of an analog one since there is no way to count 'events'
// exactly for the flux
if (settings::survival_biasing) {
// We need to account for the fact that some weight was already
// absorbed
score = p.wgt_last() + p.wgt_absorb();
} else {
score = p.wgt_last();
}
if (p.type() == Type::neutron || p.type() == Type::photon) {
score *= flux / p.macro_xs().total;
score = flux * p.wgt_last() / p.macro_xs().total;
} else {
score = 0.;
}
@ -587,13 +581,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
if (tally.estimator_ == TallyEstimator::ANALOG) {
// All events will score to the total reaction rate. We can just use
// use the weight of the particle entering the collision as the score
if (settings::survival_biasing) {
// We need to account for the fact that some weight was already
// absorbed
score = (p.wgt_last() + p.wgt_absorb()) * flux;
} else {
score = p.wgt_last() * flux;
}
score = p.wgt_last() * flux;
} else {
if (i_nuclide >= 0) {
@ -616,14 +604,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
// All events score to an inverse velocity bin. We actually use a
// collision estimator in place of an analog one since there is no way
// to count 'events' exactly for the inverse velocity
if (settings::survival_biasing) {
// We need to account for the fact that some weight was already
// absorbed
score = p.wgt_last() + p.wgt_absorb();
} else {
score = p.wgt_last();
}
score *= flux / p.macro_xs().total;
score = flux * p.wgt_last() / p.macro_xs().total;
} else {
score = flux;
}
@ -641,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() * flux;
score = (p.wgt_last() - wgt_absorb) * flux;
} else {
if (i_nuclide >= 0) {
if (p.type() == Type::neutron) {
@ -672,17 +653,17 @@ 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
if (p.event_mt() == ELASTIC || p.event_mt() == N_LEVEL ||
(p.event_mt() >= N_N1 && p.event_mt() <= N_NC)) {
// Don't waste time on very common reactions we know have
// multiplicities of one.
score = p.wgt_last() * flux;
} else {
score = (p.wgt_last() - wgt_absorb) * flux;
// Don't waste time on very common reactions we know have multiplicities
// of one.
if (p.event_mt() != ELASTIC && p.event_mt() != N_LEVEL &&
!(p.event_mt() >= N_N1 && p.event_mt() <= N_NC)) {
// Get yield and apply to score
auto m =
data::nuclides[p.event_nuclide()]->reaction_index_[p.event_mt()];
const auto& rxn {*data::nuclides[p.event_nuclide()]->reactions_[m]};
score = p.wgt_last() * flux * (*rxn.products_[0].yield_)(E);
score *= (*rxn.products_[0].yield_)(E);
}
break;
@ -694,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)
@ -725,16 +706,15 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
break;
case SCORE_FISSION:
if (p.macro_xs().absorption == 0)
if (p.macro_xs().fission == 0)
continue;
if (tally.estimator_ == TallyEstimator::ANALOG) {
if (settings::survival_biasing) {
// No fission events occur if survival biasing is on -- need to
// calculate fraction of absorptions that would have resulted in
// fission
if (p.neutron_xs(p.event_nuclide()).absorption > 0) {
score = p.wgt_absorb() * p.neutron_xs(p.event_nuclide()).fission /
p.neutron_xs(p.event_nuclide()).absorption * flux;
// No fission events occur if survival biasing is on -- use collision
// estimator instead
if (p.neutron_xs(p.event_nuclide()).total > 0) {
score = p.wgt_last() * p.neutron_xs(p.event_nuclide()).fission /
p.neutron_xs(p.event_nuclide()).total * flux;
} else {
score = 0.;
}
@ -758,7 +738,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
break;
case SCORE_NU_FISSION:
if (p.macro_xs().absorption == 0)
if (p.macro_xs().fission == 0)
continue;
if (tally.estimator_ == TallyEstimator::ANALOG) {
if (settings::survival_biasing || p.fission()) {
@ -770,13 +750,11 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
}
}
if (settings::survival_biasing) {
// No fission events occur if survival biasing is on -- need to
// calculate fraction of absorptions that would have resulted in
// nu-fission
if (p.neutron_xs(p.event_nuclide()).absorption > 0) {
score = p.wgt_absorb() *
p.neutron_xs(p.event_nuclide()).nu_fission /
p.neutron_xs(p.event_nuclide()).absorption * flux;
// No fission events occur if survival biasing is on -- use collision
// estimator instead
if (p.neutron_xs(p.event_nuclide()).total > 0) {
score = p.wgt_last() * p.neutron_xs(p.event_nuclide()).nu_fission /
p.neutron_xs(p.event_nuclide()).total * flux;
} else {
score = 0.;
}
@ -801,7 +779,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
break;
case SCORE_PROMPT_NU_FISSION:
if (p.macro_xs().absorption == 0)
if (p.macro_xs().fission == 0)
continue;
if (tally.estimator_ == TallyEstimator::ANALOG) {
if (settings::survival_biasing || p.fission()) {
@ -816,11 +794,11 @@ void score_general_ce(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
if (p.neutron_xs(p.event_nuclide()).absorption > 0) {
score = p.wgt_absorb() * p.neutron_xs(p.event_nuclide()).fission *
if (p.neutron_xs(p.event_nuclide()).total > 0) {
score = p.wgt_last() * p.neutron_xs(p.event_nuclide()).fission *
data::nuclides[p.event_nuclide()]->nu(
E, ReactionProduct::EmissionMode::prompt) /
p.neutron_xs(p.event_nuclide()).absorption * flux;
p.neutron_xs(p.event_nuclide()).total * flux;
} else {
score = 0.;
}
@ -863,7 +841,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
break;
case SCORE_DELAYED_NU_FISSION:
if (p.macro_xs().absorption == 0)
if (p.macro_xs().fission == 0)
continue;
if (tally.estimator_ == TallyEstimator::ANALOG) {
if (settings::survival_biasing || p.fission()) {
@ -878,7 +856,7 @@ void score_general_ce(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
// delayed-nu-fission
if (p.neutron_xs(p.event_nuclide()).absorption > 0 &&
if (p.neutron_xs(p.event_nuclide()).total > 0 &&
data::nuclides[p.event_nuclide()]->fissionable_) {
if (tally.delayedgroup_filter_ != C_NONE) {
auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_];
@ -890,9 +868,9 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
auto dg = filt.groups()[d_bin];
auto yield = data::nuclides[p.event_nuclide()]->nu(
E, ReactionProduct::EmissionMode::delayed, dg);
score = p.wgt_absorb() * yield *
score = p.wgt_last() * yield *
p.neutron_xs(p.event_nuclide()).fission /
p.neutron_xs(p.event_nuclide()).absorption * flux;
p.neutron_xs(p.event_nuclide()).total * flux;
score_fission_delayed_dg(
i_tally, d_bin, score, score_index, p.filter_matches());
}
@ -901,10 +879,10 @@ void score_general_ce(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() * p.neutron_xs(p.event_nuclide()).fission *
score = p.wgt_last() * p.neutron_xs(p.event_nuclide()).fission *
data::nuclides[p.event_nuclide()]->nu(
E, ReactionProduct::EmissionMode::delayed) /
p.neutron_xs(p.event_nuclide()).absorption * flux;
p.neutron_xs(p.event_nuclide()).total * flux;
}
}
} else {
@ -1007,7 +985,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
break;
case SCORE_DECAY_RATE:
if (p.macro_xs().absorption == 0)
if (p.macro_xs().fission == 0)
continue;
if (tally.estimator_ == TallyEstimator::ANALOG) {
if (settings::survival_biasing) {
@ -1015,8 +993,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
// calculate fraction of absorptions that would have resulted in
// delayed-nu-fission
const auto& nuc {*data::nuclides[p.event_nuclide()]};
if (p.neutron_xs(p.event_nuclide()).absorption > 0 &&
nuc.fissionable_) {
if (p.neutron_xs(p.event_nuclide()).total > 0 && nuc.fissionable_) {
const auto& rxn {*nuc.fission_rx_[0]};
if (tally.delayedgroup_filter_ != C_NONE) {
auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_];
@ -1029,10 +1006,9 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
auto yield =
nuc.nu(E, ReactionProduct::EmissionMode::delayed, d);
auto rate = rxn.products_[d].decay_rate_;
score = p.wgt_absorb() * yield *
score = p.wgt_last() * yield *
p.neutron_xs(p.event_nuclide()).fission /
p.neutron_xs(p.event_nuclide()).absorption * rate *
flux;
p.neutron_xs(p.event_nuclide()).total * rate * flux;
score_fission_delayed_dg(
i_tally, d_bin, score, score_index, p.filter_matches());
}
@ -1048,13 +1024,17 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
// rxn.products_ array to be exceeded. Hence, we use the size
// of this array and not the MAX_DELAYED_GROUPS constant for
// this loop.
for (auto d = 0; d < rxn.products_.size() - 2; ++d) {
for (auto d = 1; d < rxn.products_.size(); ++d) {
const auto& product = rxn.products_[d];
if (product.particle_ != Type::neutron)
continue;
auto yield =
nuc.nu(E, ReactionProduct::EmissionMode::delayed, d + 1);
auto rate = rxn.products_[d + 1].decay_rate_;
score += rate * p.wgt_absorb() *
nuc.nu(E, ReactionProduct::EmissionMode::delayed, d);
auto rate = product.decay_rate_;
score += rate * p.wgt_last() *
p.neutron_xs(p.event_nuclide()).fission * yield /
p.neutron_xs(p.event_nuclide()).absorption * flux;
p.neutron_xs(p.event_nuclide()).total * flux;
}
}
}
@ -1123,10 +1103,13 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
// rxn.products_ array to be exceeded. Hence, we use the size
// of this array and not the MAX_DELAYED_GROUPS constant for
// this loop.
for (auto d = 0; d < rxn.products_.size() - 2; ++d) {
auto yield =
nuc.nu(E, ReactionProduct::EmissionMode::delayed, d + 1);
auto rate = rxn.products_[d + 1].decay_rate_;
for (auto d = 1; d < rxn.products_.size(); ++d) {
const auto& product = rxn.products_[d];
if (product.particle_ != Type::neutron)
continue;
auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d);
auto rate = product.decay_rate_;
score += p.neutron_xs(i_nuclide).fission * flux * yield *
atom_density * rate;
}
@ -1174,10 +1157,14 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
// rxn.products_ array to be exceeded. Hence, we use the size
// of this array and not the MAX_DELAYED_GROUPS constant for
// this loop.
for (auto d = 0; d < rxn.products_.size() - 2; ++d) {
for (auto d = 1; d < rxn.products_.size(); ++d) {
const auto& product = rxn.products_[d];
if (product.particle_ != Type::neutron)
continue;
auto yield =
nuc.nu(E, ReactionProduct::EmissionMode::delayed, d + 1);
auto rate = rxn.products_[d + 1].decay_rate_;
nuc.nu(E, ReactionProduct::EmissionMode::delayed, d);
auto rate = product.decay_rate_;
score += p.neutron_xs(j_nuclide).fission * yield *
atom_density * flux * rate;
}
@ -1190,7 +1177,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
break;
case SCORE_KAPPA_FISSION:
if (p.macro_xs().absorption == 0.)
if (p.macro_xs().fission == 0.)
continue;
score = 0.;
// Kappa-fission values are determined from the Q-value listed for the
@ -1201,12 +1188,11 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
// calculate fraction of absorptions that would have resulted in
// fission scaled by the Q-value
const auto& nuc {*data::nuclides[p.event_nuclide()]};
if (p.neutron_xs(p.event_nuclide()).absorption > 0 &&
nuc.fissionable_) {
if (p.neutron_xs(p.event_nuclide()).total > 0 && nuc.fissionable_) {
const auto& rxn {*nuc.fission_rx_[0]};
score = p.wgt_absorb() * rxn.q_value_ *
score = p.wgt_last() * rxn.q_value_ *
p.neutron_xs(p.event_nuclide()).fission /
p.neutron_xs(p.event_nuclide()).absorption * flux;
p.neutron_xs(p.event_nuclide()).total * flux;
}
} else {
// Skip any non-absorption events
@ -1262,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() * flux;
score = (p.wgt_last() - wgt_absorb) * flux;
} else {
if (i_nuclide >= 0) {
if (p.neutron_xs(i_nuclide).elastic == CACHE_INVALID)
@ -1286,7 +1272,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
case SCORE_FISS_Q_PROMPT:
case SCORE_FISS_Q_RECOV:
if (p.macro_xs().absorption == 0.)
if (p.macro_xs().fission == 0.)
continue;
score =
score_fission_q(p, score_bin, tally, flux, i_nuclide, atom_density);
@ -1311,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() * flux;
score = (p.wgt_last() - wgt_absorb) * flux;
} else {
int m;
switch (score_bin) {
@ -1423,12 +1409,11 @@ void score_general_ce(Particle& p, int i_tally, int start_index,
continue;
if (tally.estimator_ == TallyEstimator::ANALOG) {
// Any other score is assumed to be a MT number. Thus, we just need
// to check if it matches the MT number of the event
if (p.event_mt() != score_bin)
continue;
score = p.wgt_last() * flux;
score = (p.wgt_last() - wgt_absorb) * flux;
} else {
// Any other cross section has to be calculated on-the-fly
if (score_bin < 2)
@ -1475,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
@ -1532,14 +1520,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
// All events score to a flux bin. We actually use a collision estimator
// in place of an analog one since there is no way to count 'events'
// exactly for the flux
if (settings::survival_biasing) {
// We need to account for the fact that some weight was already
// absorbed
score = p.wgt_last() + p.wgt_absorb();
} else {
score = p.wgt_last();
}
score *= flux / p.macro_xs().total;
score = flux * p.wgt_last() / p.macro_xs().total;
} else {
score = flux;
}
@ -1549,16 +1530,9 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
if (tally.estimator_ == TallyEstimator::ANALOG) {
// All events will score to the total reaction rate. We can just use
// use the weight of the particle entering the collision as the score
if (settings::survival_biasing) {
// We need to account for the fact that some weight was already
// absorbed
score = p.wgt_last() + p.wgt_absorb();
} else {
score = p.wgt_last();
}
// TODO: should flux be multiplied in above instead of below?
score = flux * p.wgt_last();
if (i_nuclide >= 0) {
score *= flux * atom_density * nuc_xs.get_xs(MgxsType::TOTAL, p_g) /
score *= atom_density * nuc_xs.get_xs(MgxsType::TOTAL, p_g) /
macro_xs.get_xs(MgxsType::TOTAL, p_g);
}
} else {
@ -1576,18 +1550,12 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
// All events score to an inverse velocity bin. We actually use a
// collision estimator in place of an analog one since there is no way
// to count 'events' exactly for the inverse velocity
if (settings::survival_biasing) {
// We need to account for the fact that some weight was already
// absorbed
score = p.wgt_last() + p.wgt_absorb();
} else {
score = p.wgt_last();
}
score = flux * p.wgt_last();
if (i_nuclide >= 0) {
score *= flux * nuc_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g) /
score *= nuc_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g) /
macro_xs.get_xs(MgxsType::TOTAL, p_g);
} else {
score *= flux * macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g) /
score *= macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g) /
macro_xs.get_xs(MgxsType::TOTAL, p_g);
}
} else {
@ -1606,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() * 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(),
@ -1634,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() * 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.
@ -1660,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)
@ -1689,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)
@ -1729,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);
@ -1776,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) /
@ -1837,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) /
@ -1855,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;
@ -1952,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) *
@ -1978,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,
@ -2093,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)
@ -2362,11 +2330,7 @@ void score_collision_tally(Particle& p)
// Determine the collision estimate of the flux
double flux = 0.0;
if (p.type() == ParticleType::neutron || p.type() == ParticleType::photon) {
if (!settings::survival_biasing) {
flux = p.wgt_last() / p.macro_xs().total;
} else {
flux = (p.wgt_last() + p.wgt_absorb()) / p.macro_xs().total;
}
flux = p.wgt_last() / p.macro_xs().total;
}
for (auto i_tally : model::active_collision_tallies) {

11
src/weight_windows.cpp
View file

@ -5,6 +5,7 @@
#include "openmc/hdf5_interface.h"
#include "openmc/particle.h"
#include "openmc/particle_data.h"
#include "openmc/physics_common.h"
#include "openmc/search.h"
#include "openmc/xml_interface.h"
@ -94,14 +95,8 @@ void apply_weight_windows(Particle& p)
// if the particle weight is below the window, play Russian roulette
double weight_survive =
std::min(weight * weight_window.max_split, weight_window.survival_weight);
if (weight_survive * prn(p.current_seed()) <= weight) {
p.wgt() = weight_survive;
} else {
p.alive() = false;
p.wgt() = 0.0;
}
// else particle is in the window, continue as normal
}
russian_roulette(p, weight_survive);
} // else particle is in the window, continue as normal
}
void free_memory_weight_windows()

12
tests/regression_tests/survival_biasing/results_true.dat
View file

@ -1,10 +1,10 @@
k-combined:
9.826269E-01 1.626276E-02
tally 1:
4.212303E+01
3.550855E+02
1.760388E+01
6.205773E+01
4.209907E+01
3.546281E+02
1.759415E+01
6.197412E+01
2.165888E+00
9.392975E-01
1.873459E+00
@ -16,5 +16,5 @@ tally 1:
3.628554E+08
2.635633E+16
tally 2:
1.760388E+01
6.205773E+01
1.759415E+01
6.197412E+01

2
tests/regression_tests/tallies/results_true.dat
View file

@ -1 +1 @@
de09325b517aad9f58940e5fcd53b005c57ea008a15699769ab91aba723de141014101e55be2621e34c369beefb6db8c6cf847e241ee26e73e394c0f155138b0
683a82a10c4254c7e503f9c86e192070224c0d19c1817e51da7acd74b9529a475020a6cbd658644f7cd18b244cc477d1810521aa720c4fa9354cb6955bce13f6