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

Use references consistently in physics modules

Browse source
This commit is contained in:
Paul Romano 2020-04-23 10:11:43 -05:00
parent 65c0dbfe46
commit 43eae72aeb
7 changed files with 276 additions and 277 deletions
Download patch file Download diff file Expand all files Collapse all files

32
include/openmc/physics.h
View file

@ -16,55 +16,55 @@ namespace openmc {
//==============================================================================
//! Sample a nuclide and reaction and then calls the appropriate routine
void collision(Particle* p);
void collision(Particle& p);
//! Samples an incident neutron reaction
void sample_neutron_reaction(Particle* p);
void sample_neutron_reaction(Particle& p);
//! Samples an element based on the macroscopic cross sections for each nuclide
//! within a material and then samples a reaction for that element and calls the
//! appropriate routine to process the physics.
void sample_photon_reaction(Particle* p);
void sample_photon_reaction(Particle& p);
//! Terminates the particle and either deposits all energy locally
//! (electron_treatment = ElectronTreatment::LED) or creates secondary bremsstrahlung
//! photons from electron deflections with charged particles (electron_treatment
//! = ElectronTreatment::TTB).
void sample_electron_reaction(Particle* p);
void sample_electron_reaction(Particle& p);
//! Terminates the particle and either deposits all energy locally
//! (electron_treatment = ElectronTreatment::LED) or creates secondary bremsstrahlung
//! photons from electron deflections with charged particles (electron_treatment
//! = ElectronTreatment::TTB). Two annihilation photons of energy MASS_ELECTRON_EV (0.511
//! MeV) are created and travel in opposite directions.
void sample_positron_reaction(Particle* p);
void sample_positron_reaction(Particle& p);
//! Sample a nuclide based on their total cross sections and densities within
//! the current material
//!
//! \param[in] p Particle
//! \return Index in the data::nuclides vector
int sample_nuclide(Particle* p);
int sample_nuclide(Particle& p);
//! Determine the average total, prompt, and delayed neutrons produced from
//! fission and creates appropriate bank sites.
void create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx);
void create_fission_sites(Particle& p, int i_nuclide, const Reaction* rx);
int sample_element(Particle* p);
int sample_element(Particle& p);
Reaction* sample_fission(int i_nuclide, Particle* p);
Reaction* sample_fission(int i_nuclide, Particle& p);
void sample_photon_product(int i_nuclide, Particle* p, int* i_rx, int* i_product);
void sample_photon_product(int i_nuclide, Particle& p, int* i_rx, int* i_product);
void absorption(Particle* p, int i_nuclide);
void absorption(Particle& p, int i_nuclide);
void scatter(Particle* p, int i_nuclide);
void scatter(Particle& p, int i_nuclide);
//! Treats the elastic scattering of a neutron with a target.
void elastic_scatter(int i_nuclide, const Reaction& rx, double kT,
Particle* p);
Particle& p);
void sab_scatter(int i_nuclide, int i_sab, Particle* p);
void sab_scatter(int i_nuclide, int i_sab, Particle& p);
//! samples the target velocity. The constant cross section free gas model is
//! the default method. Methods for correctly accounting for the energy
@ -85,9 +85,9 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in,
//! handles all reactions with a single secondary neutron (other than fission),
//! i.e. level scattering, (n,np), (n,na), etc.
void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p);
void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle& p);
void sample_secondary_photons(Particle* p, int i_nuclide);
void sample_secondary_photons(Particle& p, int i_nuclide);
} // namespace openmc

3
include/openmc/physics_common.h
View file

@ -9,8 +9,7 @@
namespace openmc {
//! \brief Performs the russian roulette operation for a particle
extern "C" void
russian_roulette(Particle* p);
void russian_roulette(Particle& p);
} // namespace openmc
#endif // OPENMC_PHYSICS_COMMON_H

10
include/openmc/physics_mg.h
View file

@ -15,7 +15,7 @@ namespace openmc {
//! \brief samples particle behavior after a collision event.
//! \param p Particle to operate on
void
collision_mg(Particle* p);
collision_mg(Particle& p);
//! \brief samples a reaction type.
//!
@ -23,23 +23,23 @@ collision_mg(Particle* p);
//! fission and disappearance are treated implicitly.
//! \param p Particle to operate on
void
sample_reaction(Particle* p);
sample_reaction(Particle& p);
//! \brief Samples the scattering event
//! \param p Particle to operate on
void
scatter(Particle* p);
scatter(Particle& p);
//! \brief Determines the average total, prompt and delayed neutrons produced
//! from fission and creates the appropriate bank sites.
//! \param p Particle to operate on
void
create_fission_sites(Particle* p);
create_fission_sites(Particle& p);
//! \brief Handles an absorption event
//! \param p Particle to operate on
void
absorption(Particle* p);
absorption(Particle& p);
} // namespace openmc
#endif // OPENMC_PHYSICS_MG_H

4
src/particle.cpp
View file

@ -289,9 +289,9 @@ Particle::event_collide()
surface_ = 0;
if (settings::run_CE) {
collision(this);
collision(*this);
} else {
collision_mg(this);
collision_mg(*this);
}
// Score collision estimator tallies -- this is done after a collision

388
src/physics.cpp
View file

@ -32,13 +32,13 @@ namespace openmc {
// Non-member functions
//==============================================================================
void collision(Particle* p)
void collision(Particle& p)
{
// Add to collision counter for particle
++(p->n_collision_);
++(p.n_collision_);
// Sample reaction for the material the particle is in
switch (static_cast<Particle::Type>(p->type_)) {
switch (static_cast<Particle::Type>(p.type_)) {
case Particle::Type::neutron:
sample_neutron_reaction(p);
break;
@ -54,39 +54,39 @@ void collision(Particle* p)
}
// Kill particle if energy falls below cutoff
int type = static_cast<int>(p->type_);
if (p->E_ < settings::energy_cutoff[type]) {
p->alive_ = false;
p->wgt_ = 0.0;
int type = static_cast<int>(p.type_);
if (p.E_ < settings::energy_cutoff[type]) {
p.alive_ = false;
p.wgt_ = 0.0;
}
// Display information about collision
if (settings::verbosity >= 10 || p->trace_) {
if (settings::verbosity >= 10 || p.trace_) {
std::string msg;
if (p->event_ == TallyEvent::KILL) {
msg = fmt::format(" Killed. Energy = {} eV.", p->E_);
} else if (p->type_ == Particle::Type::neutron) {
if (p.event_ == TallyEvent::KILL) {
msg = fmt::format(" Killed. Energy = {} eV.", p.E_);
} else if (p.type_ == Particle::Type::neutron) {
msg = fmt::format(" {} with {}. Energy = {} eV.",
reaction_name(p->event_mt_), data::nuclides[p->event_nuclide_]->name_,
p->E_);
} else if (p->type_ == Particle::Type::photon) {
reaction_name(p.event_mt_), data::nuclides[p.event_nuclide_]->name_,
p.E_);
} else if (p.type_ == Particle::Type::photon) {
msg = fmt::format(" {} with {}. Energy = {} eV.",
reaction_name(p->event_mt_),
to_element(data::nuclides[p->event_nuclide_]->name_), p->E_);
reaction_name(p.event_mt_),
to_element(data::nuclides[p.event_nuclide_]->name_), p.E_);
} else {
msg = fmt::format(" Disappeared. Energy = {} eV.", p->E_);
msg = fmt::format(" Disappeared. Energy = {} eV.", p.E_);
}
write_message(msg, 1);
}
}
void sample_neutron_reaction(Particle* p)
void sample_neutron_reaction(Particle& p)
{
// Sample a nuclide within the material
int i_nuclide = sample_nuclide(p);
// Save which nuclide particle had collision with
p->event_nuclide_ = i_nuclide;
p.event_nuclide_ = i_nuclide;
// Create fission bank sites. Note that while a fission reaction is sampled,
// it never actually "happens", i.e. the weight of the particle does not
@ -105,7 +105,7 @@ void sample_neutron_reaction(Particle* p)
// Make sure particle population doesn't grow out of control for
// subcritical multiplication problems.
if (p->secondary_bank_.size() >= 10000) {
if (p.secondary_bank_.size() >= 10000) {
fatal_error("The secondary particle bank appears to be growing without "
"bound. You are likely running a subcritical multiplication problem "
"with k-effective close to or greater than one.");
@ -115,53 +115,53 @@ void sample_neutron_reaction(Particle* p)
// Create secondary photons
if (settings::photon_transport) {
p->stream_ = STREAM_PHOTON;
p.stream_ = STREAM_PHOTON;
sample_secondary_photons(p, i_nuclide);
p->stream_ = STREAM_TRACKING;
p.stream_ = STREAM_TRACKING;
}
// If survival biasing is being used, the following subroutine adjusts the
// weight of the particle. Otherwise, it checks to see if absorption occurs
if (p->neutron_xs_[i_nuclide].absorption > 0.0) {
if (p.neutron_xs_[i_nuclide].absorption > 0.0) {
absorption(p, i_nuclide);
} else {
p->wgt_absorb_ = 0.0;
p.wgt_absorb_ = 0.0;
}
if (!p->alive_) return;
if (!p.alive_) return;
// Sample a scattering reaction and determine the secondary energy of the
// exiting neutron
scatter(p, i_nuclide);
// Advance URR seed stream 'N' times after energy changes
if (p->E_ != p->E_last_) {
p->stream_ = STREAM_URR_PTABLE;
advance_prn_seed(data::nuclides.size(), p->current_seed());
p->stream_ = STREAM_TRACKING;
if (p.E_ != p.E_last_) {
p.stream_ = STREAM_URR_PTABLE;
advance_prn_seed(data::nuclides.size(), p.current_seed());
p.stream_ = STREAM_TRACKING;
}
// Play russian roulette if survival biasing is turned on
if (settings::survival_biasing) {
russian_roulette(p);
if (!p->alive_) return;
if (!p.alive_) return;
}
}
void
create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx)
create_fission_sites(Particle& p, int i_nuclide, const Reaction* rx)
{
// If uniform fission source weighting is turned on, we increase or decrease
// the expected number of fission sites produced
double weight = settings::ufs_on ? ufs_get_weight(*p) : 1.0;
double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0;
// Determine the expected number of neutrons produced
double nu_t = p->wgt_ / simulation::keff * weight * p->neutron_xs_[
i_nuclide].nu_fission / p->neutron_xs_[i_nuclide].total;
double nu_t = p.wgt_ / simulation::keff * weight * p.neutron_xs_[
i_nuclide].nu_fission / p.neutron_xs_[i_nuclide].total;
// Sample the number of neutrons produced
int nu = static_cast<int>(nu_t);
if (prn(p->current_seed()) <= (nu_t - nu)) ++nu;
if (prn(p.current_seed()) <= (nu_t - nu)) ++nu;
// Begin banking the source neutrons
// First, if our bank is full then don't continue
@ -172,9 +172,9 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx)
double nu_d[MAX_DELAYED_GROUPS] = {0.};
// Clear out particle's nu fission bank
p->nu_bank_.clear();
p.nu_bank_.clear();
p->fission_ = true;
p.fission_ = true;
int skipped = 0;
// Determine whether to place fission sites into the shared fission bank
@ -184,14 +184,14 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx)
for (int i = 0; i < nu; ++i) {
// Initialize fission site object with particle data
Particle::Bank site;
site.r = p->r();
site.r = p.r();
site.particle = Particle::Type::neutron;
site.wgt = 1. / weight;
site.parent_id = p->id_;
site.progeny_id = p->n_progeny_++;
site.parent_id = p.id_;
site.progeny_id = p.n_progeny_++;
// Sample delayed group and angle/energy for fission reaction
sample_fission_neutron(i_nuclide, rx, p->E_, &site, p->current_seed());
sample_fission_neutron(i_nuclide, rx, p.E_, &site, p.current_seed());
// Store fission site in bank
if (use_fission_bank) {
@ -203,21 +203,21 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx)
break;
}
} else {
p->secondary_bank_.push_back(site);
p.secondary_bank_.push_back(site);
}
// Set the delayed group on the particle as well
p->delayed_group_ = site.delayed_group;
p.delayed_group_ = site.delayed_group;
// Increment the number of neutrons born delayed
if (p->delayed_group_ > 0) {
nu_d[p->delayed_group_-1]++;
if (p.delayed_group_ > 0) {
nu_d[p.delayed_group_-1]++;
}
// Write fission particles to nuBank
if (use_fission_bank) {
p->nu_bank_.emplace_back();
Particle::NuBank* nu_bank_entry = &p->nu_bank_.back();
p.nu_bank_.emplace_back();
Particle::NuBank* nu_bank_entry = &p.nu_bank_.back();
nu_bank_entry->wgt = site.wgt;
nu_bank_entry->E = site.E;
nu_bank_entry->delayed_group = site.delayed_group;
@ -227,7 +227,7 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx)
// If shared fission bank was full, and no fissions could be added,
// set the particle fission flag to false.
if (nu == skipped) {
p->fission_ = false;
p.fission_ = false;
return;
}
@ -236,45 +236,45 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx)
nu -= skipped;
// Store the total weight banked for analog fission tallies
p->n_bank_ = nu;
p->wgt_bank_ = nu / weight;
p.n_bank_ = nu;
p.wgt_bank_ = nu / weight;
for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) {
p->n_delayed_bank_[d] = nu_d[d];
p.n_delayed_bank_[d] = nu_d[d];
}
}
void sample_photon_reaction(Particle* p)
void sample_photon_reaction(Particle& p)
{
// Kill photon if below energy cutoff -- an extra check is made here because
// photons with energy below the cutoff may have been produced by neutrons
// reactions or atomic relaxation
int photon = static_cast<int>(Particle::Type::photon);
if (p->E_ < settings::energy_cutoff[photon]) {
p->E_ = 0.0;
p->alive_ = false;
if (p.E_ < settings::energy_cutoff[photon]) {
p.E_ = 0.0;
p.alive_ = false;
return;
}
// Sample element within material
int i_element = sample_element(p);
const auto& micro {p->photon_xs_[i_element]};
const auto& micro {p.photon_xs_[i_element]};
const auto& element {data::elements[i_element]};
// Calculate photon energy over electron rest mass equivalent
double alpha = p->E_/MASS_ELECTRON_EV;
double alpha = p.E_/MASS_ELECTRON_EV;
// For tallying purposes, this routine might be called directly. In that
// case, we need to sample a reaction via the cutoff variable
double prob = 0.0;
double cutoff = prn(p->current_seed()) * micro.total;
double cutoff = prn(p.current_seed()) * micro.total;
// Coherent (Rayleigh) scattering
prob += micro.coherent;
if (prob > cutoff) {
double mu = element.rayleigh_scatter(alpha, p->current_seed());
p->u() = rotate_angle(p->u(), mu, nullptr, p->current_seed());
p->event_ = TallyEvent::SCATTER;
p->event_mt_ = COHERENT;
double mu = element.rayleigh_scatter(alpha, p.current_seed());
p.u() = rotate_angle(p.u(), mu, nullptr, p.current_seed());
p.event_ = TallyEvent::SCATTER;
p.event_mt_ = COHERENT;
return;
}
@ -283,7 +283,7 @@ void sample_photon_reaction(Particle* p)
if (prob > cutoff) {
double alpha_out, mu;
int i_shell;
element.compton_scatter(alpha, true, &alpha_out, &mu, &i_shell, p->current_seed());
element.compton_scatter(alpha, true, &alpha_out, &mu, &i_shell, p.current_seed());
// Determine binding energy of shell. The binding energy is 0.0 if
// doppler broadening is not used.
@ -295,14 +295,14 @@ void sample_photon_reaction(Particle* p)
}
// Create Compton electron
double phi = 2.0*PI*prn(p->current_seed());
double phi = 2.0*PI*prn(p.current_seed());
double E_electron = (alpha - alpha_out)*MASS_ELECTRON_EV - e_b;
int electron = static_cast<int>(Particle::Type::electron);
if (E_electron >= settings::energy_cutoff[electron]) {
double mu_electron = (alpha - alpha_out*mu)
/ std::sqrt(alpha*alpha + alpha_out*alpha_out - 2.0*alpha*alpha_out*mu);
Direction u = rotate_angle(p->u(), mu_electron, &phi, p->current_seed());
p->create_secondary(u, E_electron, Particle::Type::electron);
Direction u = rotate_angle(p.u(), mu_electron, &phi, p.current_seed());
p.create_secondary(u, E_electron, Particle::Type::electron);
}
// TODO: Compton subshell data does not match atomic relaxation data
@ -310,14 +310,14 @@ void sample_photon_reaction(Particle* p)
// and fluorescent photons
if (i_shell >= 0) {
const auto& shell = element.shells_[i_shell];
element.atomic_relaxation(shell, *p);
element.atomic_relaxation(shell, p);
}
phi += PI;
p->E_ = alpha_out*MASS_ELECTRON_EV;
p->u() = rotate_angle(p->u(), mu, &phi, p->current_seed());
p->event_ = TallyEvent::SCATTER;
p->event_mt_ = INCOHERENT;
p.E_ = alpha_out*MASS_ELECTRON_EV;
p.u() = rotate_angle(p.u(), mu, &phi, p.current_seed());
p.event_ = TallyEvent::SCATTER;
p.event_mt_ = INCOHERENT;
return;
}
@ -340,15 +340,15 @@ void sample_photon_reaction(Particle* p)
prob += xs;
if (prob > cutoff) {
double E_electron = p->E_ - shell.binding_energy;
double E_electron = p.E_ - shell.binding_energy;
// Sample mu using non-relativistic Sauter distribution.
// See Eqns 3.19 and 3.20 in "Implementing a photon physics
// model in Serpent 2" by Toni Kaltiaisenaho
double mu;
while (true) {
double r = prn(p->current_seed());
if (4.0*(1.0 - r)*r >= prn(p->current_seed())) {
double r = prn(p.current_seed());
if (4.0*(1.0 - r)*r >= prn(p.current_seed())) {
double rel_vel = std::sqrt(E_electron * (E_electron +
2.0*MASS_ELECTRON_EV)) / (E_electron + MASS_ELECTRON_EV);
mu = (2.0*r + rel_vel - 1.0) / (2.0*rel_vel*r - rel_vel + 1.0);
@ -356,22 +356,22 @@ void sample_photon_reaction(Particle* p)
}
}
double phi = 2.0*PI*prn(p->current_seed());
double phi = 2.0*PI*prn(p.current_seed());
Direction u;
u.x = mu;
u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
// Create secondary electron
p->create_secondary(u, E_electron, Particle::Type::electron);
p.create_secondary(u, E_electron, Particle::Type::electron);
// Allow electrons to fill orbital and produce auger electrons
// and fluorescent photons
element.atomic_relaxation(shell, *p);
p->event_ = TallyEvent::ABSORB;
p->event_mt_ = 533 + shell.index_subshell;
p->alive_ = false;
p->E_ = 0.0;
element.atomic_relaxation(shell, p);
p.event_ = TallyEvent::ABSORB;
p.event_mt_ = 533 + shell.index_subshell;
p.alive_ = false;
p.E_ = 0.0;
return;
}
}
@ -384,70 +384,70 @@ void sample_photon_reaction(Particle* p)
double E_electron, E_positron;
double mu_electron, mu_positron;
element.pair_production(alpha, &E_electron, &E_positron,
&mu_electron, &mu_positron, p->current_seed());
&mu_electron, &mu_positron, p.current_seed());
// Create secondary electron
Direction u = rotate_angle(p->u(), mu_electron, nullptr, p->current_seed());
p->create_secondary(u, E_electron, Particle::Type::electron);
Direction u = rotate_angle(p.u(), mu_electron, nullptr, p.current_seed());
p.create_secondary(u, E_electron, Particle::Type::electron);
// Create secondary positron
u = rotate_angle(p->u(), mu_positron, nullptr, p->current_seed());
p->create_secondary(u, E_positron, Particle::Type::positron);
u = rotate_angle(p.u(), mu_positron, nullptr, p.current_seed());
p.create_secondary(u, E_positron, Particle::Type::positron);
p->event_ = TallyEvent::ABSORB;
p->event_mt_ = PAIR_PROD;
p->alive_ = false;
p->E_ = 0.0;
p.event_ = TallyEvent::ABSORB;
p.event_mt_ = PAIR_PROD;
p.alive_ = false;
p.E_ = 0.0;
}
}
void sample_electron_reaction(Particle* p)
void sample_electron_reaction(Particle& p)
{
// TODO: create reaction types
if (settings::electron_treatment == ElectronTreatment::TTB) {
double E_lost;
thick_target_bremsstrahlung(*p, &E_lost);
thick_target_bremsstrahlung(p, &E_lost);
}
p->E_ = 0.0;
p->alive_ = false;
p->event_ = TallyEvent::ABSORB;
p.E_ = 0.0;
p.alive_ = false;
p.event_ = TallyEvent::ABSORB;
}
void sample_positron_reaction(Particle* p)
void sample_positron_reaction(Particle& p)
{
// TODO: create reaction types
if (settings::electron_treatment == ElectronTreatment::TTB) {
double E_lost;
thick_target_bremsstrahlung(*p, &E_lost);
thick_target_bremsstrahlung(p, &E_lost);
}
// Sample angle isotropically
double mu = 2.0*prn(p->current_seed()) - 1.0;
double phi = 2.0*PI*prn(p->current_seed());
double mu = 2.0*prn(p.current_seed()) - 1.0;
double phi = 2.0*PI*prn(p.current_seed());
Direction u;
u.x = mu;
u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
// Create annihilation photon pair traveling in opposite directions
p->create_secondary(u, MASS_ELECTRON_EV, Particle::Type::photon);
p->create_secondary(-u, MASS_ELECTRON_EV, Particle::Type::photon);
p.create_secondary(u, MASS_ELECTRON_EV, Particle::Type::photon);
p.create_secondary(-u, MASS_ELECTRON_EV, Particle::Type::photon);
p->E_ = 0.0;
p->alive_ = false;
p->event_ = TallyEvent::ABSORB;
p.E_ = 0.0;
p.alive_ = false;
p.event_ = TallyEvent::ABSORB;
}
int sample_nuclide(Particle* p)
int sample_nuclide(Particle& p)
{
// Sample cumulative distribution function
double cutoff = prn(p->current_seed()) * p->macro_xs_.total;
double cutoff = prn(p.current_seed()) * p.macro_xs_.total;
// Get pointers to nuclide/density arrays
const auto& mat {model::materials[p->material_]};
const auto& mat {model::materials[p.material_]};
int n = mat->nuclide_.size();
double prob = 0.0;
@ -457,22 +457,22 @@ int sample_nuclide(Particle* p)
double atom_density = mat->atom_density_[i];
// Increment probability to compare to cutoff
prob += atom_density * p->neutron_xs_[i_nuclide].total;
prob += atom_density * p.neutron_xs_[i_nuclide].total;
if (prob >= cutoff) return i_nuclide;
}
// If we reach here, no nuclide was sampled
p->write_restart();
p.write_restart();
throw std::runtime_error{"Did not sample any nuclide during collision."};
}
int sample_element(Particle* p)
int sample_element(Particle& p)
{
// Sample cumulative distribution function
double cutoff = prn(p->current_seed()) * p->macro_xs_.total;
double cutoff = prn(p.current_seed()) * p.macro_xs_.total;
// Get pointers to elements, densities
const auto& mat {model::materials[p->material_]};
const auto& mat {model::materials[p.material_]};
double prob = 0.0;
for (int i = 0; i < mat->element_.size(); ++i) {
@ -481,24 +481,24 @@ int sample_element(Particle* p)
double atom_density = mat->atom_density_[i];
// Determine microscopic cross section
double sigma = atom_density * p->photon_xs_[i_element].total;
double sigma = atom_density * p.photon_xs_[i_element].total;
// Increment probability to compare to cutoff
prob += sigma;
if (prob > cutoff) {
// Save which nuclide particle had collision with for tally purpose
p->event_nuclide_ = mat->nuclide_[i];
p.event_nuclide_ = mat->nuclide_[i];
return i_element;
}
}
// If we made it here, no element was sampled
p->write_restart();
p.write_restart();
fatal_error("Did not sample any element during collision.");
}
Reaction* sample_fission(int i_nuclide, Particle* p)
Reaction* sample_fission(int i_nuclide, Particle& p)
{
// Get pointer to nuclide
const auto& nuc {data::nuclides[i_nuclide]};
@ -506,23 +506,23 @@ Reaction* sample_fission(int i_nuclide, Particle* p)
// If we're in the URR, by default use the first fission reaction. We also
// default to the first reaction if we know that there are no partial fission
// reactions
if (p->neutron_xs_[i_nuclide].use_ptable || !nuc->has_partial_fission_) {
if (p.neutron_xs_[i_nuclide].use_ptable || !nuc->has_partial_fission_) {
return nuc->fission_rx_[0];
}
// Check to see if we are in a windowed multipole range. WMP only supports
// the first fission reaction.
if (nuc->multipole_) {
if (p->E_ >= nuc->multipole_->E_min_ && p->E_ <= nuc->multipole_->E_max_) {
if (p.E_ >= nuc->multipole_->E_min_ && p.E_ <= nuc->multipole_->E_max_) {
return nuc->fission_rx_[0];
}
}
// Get grid index and interpolatoin factor and sample fission cdf
int i_temp = p->neutron_xs_[i_nuclide].index_temp;
int i_grid = p->neutron_xs_[i_nuclide].index_grid;
double f = p->neutron_xs_[i_nuclide].interp_factor;
double cutoff = prn(p->current_seed()) * p->neutron_xs_[i_nuclide].fission;
int i_temp = p.neutron_xs_[i_nuclide].index_temp;
int i_grid = p.neutron_xs_[i_nuclide].index_grid;
double f = p.neutron_xs_[i_nuclide].interp_factor;
double cutoff = prn(p.current_seed()) * p.neutron_xs_[i_nuclide].fission;
double prob = 0.0;
// Loop through each partial fission reaction type
@ -543,13 +543,13 @@ Reaction* sample_fission(int i_nuclide, Particle* p)
throw std::runtime_error{"No fission reaction was sampled for " + nuc->name_};
}
void sample_photon_product(int i_nuclide, Particle* p, int* i_rx, int* i_product)
void sample_photon_product(int i_nuclide, Particle& p, int* i_rx, int* i_product)
{
// Get grid index and interpolation factor and sample photon production cdf
int i_temp = p->neutron_xs_[i_nuclide].index_temp;
int i_grid = p->neutron_xs_[i_nuclide].index_grid;
double f = p->neutron_xs_[i_nuclide].interp_factor;
double cutoff = prn(p->current_seed()) * p->neutron_xs_[i_nuclide].photon_prod;
int i_temp = p.neutron_xs_[i_nuclide].index_temp;
int i_grid = p.neutron_xs_[i_nuclide].index_grid;
double f = p.neutron_xs_[i_nuclide].interp_factor;
double cutoff = prn(p.current_seed()) * p.neutron_xs_[i_nuclide].photon_prod;
double prob = 0.0;
// Loop through each reaction type
@ -568,7 +568,7 @@ void sample_photon_product(int i_nuclide, Particle* p, int* i_rx, int* i_product
for (int j = 0; j < rx->products_.size(); ++j) {
if (rx->products_[j].particle_ == Particle::Type::photon) {
// add to cumulative probability
prob += (*rx->products_[j].yield_)(p->E_) * xs;
prob += (*rx->products_[j].yield_)(p.E_) * xs;
*i_rx = i;
*i_product = j;
@ -578,59 +578,59 @@ void sample_photon_product(int i_nuclide, Particle* p, int* i_rx, int* i_product
}
}
void absorption(Particle* p, int i_nuclide)
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;
p.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_ -= p.wgt_absorb_;
p.wgt_last_ = p.wgt_;
// Score implicit absorption estimate of keff
if (settings::run_mode == RunMode::EIGENVALUE) {
p->keff_tally_absorption_ += p->wgt_absorb_ * p->neutron_xs_[
i_nuclide].nu_fission / p->neutron_xs_[i_nuclide].absorption;
p.keff_tally_absorption_ += p.wgt_absorb_ * p.neutron_xs_[
i_nuclide].nu_fission / p.neutron_xs_[i_nuclide].absorption;
}
} else {
// See if disappearance reaction happens
if (p->neutron_xs_[i_nuclide].absorption >
prn(p->current_seed()) * p->neutron_xs_[i_nuclide].total) {
if (p.neutron_xs_[i_nuclide].absorption >
prn(p.current_seed()) * p.neutron_xs_[i_nuclide].total) {
// Score absorption estimate of keff
if (settings::run_mode == RunMode::EIGENVALUE) {
p->keff_tally_absorption_ += p->wgt_ * p->neutron_xs_[
i_nuclide].nu_fission / p->neutron_xs_[i_nuclide].absorption;
p.keff_tally_absorption_ += p.wgt_ * p.neutron_xs_[
i_nuclide].nu_fission / p.neutron_xs_[i_nuclide].absorption;
}
p->alive_ = false;
p->event_ = TallyEvent::ABSORB;
p->event_mt_ = N_DISAPPEAR;
p.alive_ = false;
p.event_ = TallyEvent::ABSORB;
p.event_mt_ = N_DISAPPEAR;
}
}
}
void scatter(Particle* p, int i_nuclide)
void scatter(Particle& p, int i_nuclide)
{
// copy incoming direction
Direction u_old {p->u()};
Direction u_old {p.u()};
// Get pointer to nuclide and grid index/interpolation factor
const auto& nuc {data::nuclides[i_nuclide]};
const auto& micro {p->neutron_xs_[i_nuclide]};
const auto& micro {p.neutron_xs_[i_nuclide]};
int i_temp = micro.index_temp;
int i_grid = micro.index_grid;
double f = micro.interp_factor;
// For tallying purposes, this routine might be called directly. In that
// case, we need to sample a reaction via the cutoff variable
double cutoff = prn(p->current_seed()) * (micro.total - micro.absorption);
double cutoff = prn(p.current_seed()) * (micro.total - micro.absorption);
bool sampled = false;
// Calculate elastic cross section if it wasn't precalculated
if (micro.elastic == CACHE_INVALID) {
nuc->calculate_elastic_xs(*p);
nuc->calculate_elastic_xs(p);
}
double prob = micro.elastic - micro.thermal;
@ -639,12 +639,12 @@ void scatter(Particle* p, int i_nuclide)
// NON-S(A,B) ELASTIC SCATTERING
// Determine temperature
double kT = nuc->multipole_ ? p->sqrtkT_*p->sqrtkT_ : nuc->kTs_[i_temp];
double kT = nuc->multipole_ ? p.sqrtkT_*p.sqrtkT_ : nuc->kTs_[i_temp];
// Perform collision physics for elastic scattering
elastic_scatter(i_nuclide, *nuc->reactions_[0], kT, p);
p->event_mt_ = ELASTIC;
p.event_mt_ = ELASTIC;
sampled = true;
}
@ -655,7 +655,7 @@ void scatter(Particle* p, int i_nuclide)
sab_scatter(i_nuclide, micro.index_sab, p);
p->event_mt_ = ELASTIC;
p.event_mt_ = ELASTIC;
sampled = true;
}
@ -671,7 +671,7 @@ void scatter(Particle* p, int i_nuclide)
// Check to make sure inelastic scattering reaction sampled
if (i >= nuc->reactions_.size()) {
p->write_restart();
p.write_restart();
fatal_error("Did not sample any reaction for nuclide " + nuc->name_);
}
@ -687,47 +687,47 @@ void scatter(Particle* p, int i_nuclide)
// Perform collision physics for inelastic scattering
const auto& rx {nuc->reactions_[i]};
inelastic_scatter(nuc.get(), rx.get(), p);
p->event_mt_ = rx->mt_;
p.event_mt_ = rx->mt_;
}
// Set event component
p->event_ = TallyEvent::SCATTER;
p.event_ = TallyEvent::SCATTER;
// Sample new outgoing angle for isotropic-in-lab scattering
const auto& mat {model::materials[p->material_]};
const auto& mat {model::materials[p.material_]};
if (!mat->p0_.empty()) {
int i_nuc_mat = mat->mat_nuclide_index_[i_nuclide];
if (mat->p0_[i_nuc_mat]) {
// Sample isotropic-in-lab outgoing direction
double mu = 2.0*prn(p->current_seed()) - 1.0;
double phi = 2.0*PI*prn(p->current_seed());
double mu = 2.0*prn(p.current_seed()) - 1.0;
double phi = 2.0*PI*prn(p.current_seed());
// Change direction of particle
p->u().x = mu;
p->u().y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
p->u().z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
p->mu_ = u_old.dot(p->u());
p.u().x = mu;
p.u().y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
p.u().z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
p.mu_ = u_old.dot(p.u());
}
}
}
void elastic_scatter(int i_nuclide, const Reaction& rx, double kT,
Particle* p)
Particle& p)
{
// get pointer to nuclide
const auto& nuc {data::nuclides[i_nuclide]};
double vel = std::sqrt(p->E_);
double vel = std::sqrt(p.E_);
double awr = nuc->awr_;
// Neutron velocity in LAB
Direction v_n = vel*p->u();
Direction v_n = vel*p.u();
// Sample velocity of target nucleus
Direction v_t {};
if (!p->neutron_xs_[i_nuclide].use_ptable) {
v_t = sample_target_velocity(nuc.get(), p->E_, p->u(), v_n,
p->neutron_xs_[i_nuclide].elastic, kT, p->current_seed());
if (!p.neutron_xs_[i_nuclide].use_ptable) {
v_t = sample_target_velocity(nuc.get(), p.E_, p.u(), v_n,
p.neutron_xs_[i_nuclide].elastic, kT, p.current_seed());
}
// Velocity of center-of-mass
@ -745,9 +745,9 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT,
auto& d = rx.products_[0].distribution_[0];
auto d_ = dynamic_cast<UncorrelatedAngleEnergy*>(d.get());
if (!d_->angle().empty()) {
mu_cm = d_->angle().sample(p->E_, p->current_seed());
mu_cm = d_->angle().sample(p.E_, p.current_seed());
} else {
mu_cm = 2.0*prn(p->current_seed()) - 1.0;
mu_cm = 2.0*prn(p.current_seed()) - 1.0;
}
// Determine direction cosines in CM
@ -756,40 +756,40 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT,
// Rotate neutron velocity vector to new angle -- note that the speed of the
// neutron in CM does not change in elastic scattering. However, the speed
// will change when we convert back to LAB
v_n = vel * rotate_angle(u_cm, mu_cm, nullptr, p->current_seed());
v_n = vel * rotate_angle(u_cm, mu_cm, nullptr, p.current_seed());
// Transform back to LAB frame
v_n += v_cm;
p->E_ = v_n.dot(v_n);
vel = std::sqrt(p->E_);
p.E_ = v_n.dot(v_n);
vel = std::sqrt(p.E_);
// compute cosine of scattering angle in LAB frame by taking dot product of
// neutron's pre- and post-collision angle
p->mu_ = p->u().dot(v_n) / vel;
p.mu_ = p.u().dot(v_n) / vel;
// Set energy and direction of particle in LAB frame
p->u() = v_n / vel;
p.u() = v_n / vel;
// Because of floating-point roundoff, it may be possible for mu_lab to be
// outside of the range [-1,1). In these cases, we just set mu_lab to exactly
// -1 or 1
if (std::abs(p->mu_) > 1.0) p->mu_ = std::copysign(1.0, p->mu_);
if (std::abs(p.mu_) > 1.0) p.mu_ = std::copysign(1.0, p.mu_);
}
void sab_scatter(int i_nuclide, int i_sab, Particle* p)
void sab_scatter(int i_nuclide, int i_sab, Particle& p)
{
// Determine temperature index
const auto& micro {p->neutron_xs_[i_nuclide]};
const auto& micro {p.neutron_xs_[i_nuclide]};
int i_temp = micro.index_temp_sab;
// Sample energy and angle
double E_out;
data::thermal_scatt[i_sab]->data_[i_temp].sample(micro, p->E_, &E_out, &p->mu_, p->current_seed());
data::thermal_scatt[i_sab]->data_[i_temp].sample(micro, p.E_, &E_out, &p.mu_, p.current_seed());
// Set energy to outgoing, change direction of particle
p->E_ = E_out;
p->u() = rotate_angle(p->u(), p->mu_, nullptr, p->current_seed());
p.E_ = E_out;
p.u() = rotate_angle(p.u(), p.mu_, nullptr, p.current_seed());
}
Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
@ -1085,15 +1085,15 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Part
}
}
void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p)
void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle& p)
{
// copy energy of neutron
double E_in = p->E_;
double E_in = p.E_;
// sample outgoing energy and scattering cosine
double E;
double mu;
rx->products_[0].sample(E_in, E, mu, p->current_seed());
rx->products_[0].sample(E_in, E, mu, p.current_seed());
// if scattering system is in center-of-mass, transfer cosine of scattering
// angle and outgoing energy from CM to LAB
@ -1115,32 +1115,32 @@ void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p)
if (std::abs(mu) > 1.0) mu = std::copysign(1.0, mu);
// Set outgoing energy and scattering angle
p->E_ = E;
p->mu_ = mu;
p.E_ = E;
p.mu_ = mu;
// change direction of particle
p->u() = rotate_angle(p->u(), mu, nullptr, p->current_seed());
p.u() = rotate_angle(p.u(), mu, nullptr, p.current_seed());
// evaluate yield
double yield = (*rx->products_[0].yield_)(E_in);
if (std::floor(yield) == yield) {
// If yield is integral, create exactly that many secondary particles
for (int i = 0; i < static_cast<int>(std::round(yield)) - 1; ++i) {
p->create_secondary(p->u(), p->E_, Particle::Type::neutron);
p.create_secondary(p.u(), p.E_, Particle::Type::neutron);
}
} else {
// Otherwise, change weight of particle based on yield
p->wgt_ *= yield;
p.wgt_ *= yield;
}
}
void sample_secondary_photons(Particle* p, int i_nuclide)
void sample_secondary_photons(Particle& p, int i_nuclide)
{
// Sample the number of photons produced
double y_t = p->wgt_ * p->neutron_xs_[i_nuclide].photon_prod /
p->neutron_xs_[i_nuclide].total;
double y_t = p.wgt_ * p.neutron_xs_[i_nuclide].photon_prod /
p.neutron_xs_[i_nuclide].total;
int y = static_cast<int>(y_t);
if (prn(p->current_seed()) <= y_t - y) ++y;
if (prn(p.current_seed()) <= y_t - y) ++y;
// Sample each secondary photon
for (int i = 0; i < y; ++i) {
@ -1153,13 +1153,13 @@ void sample_secondary_photons(Particle* p, int i_nuclide)
auto& rx = data::nuclides[i_nuclide]->reactions_[i_rx];
double E;
double mu;
rx->products_[i_product].sample(p->E_, E, mu, p->current_seed());
rx->products_[i_product].sample(p.E_, E, mu, p.current_seed());
// Sample the new direction
Direction u = rotate_angle(p->u(), mu, nullptr, p->current_seed());
Direction u = rotate_angle(p.u(), mu, nullptr, p.current_seed());
// Create the secondary photon
p->create_secondary(u, E, Particle::Type::photon);
p.create_secondary(u, E, Particle::Type::photon);
}
}

16
src/physics_common.cpp
View file

@ -9,16 +9,16 @@ namespace openmc {
// RUSSIAN_ROULETTE
//==============================================================================
void russian_roulette(Particle* p)
void russian_roulette(Particle& p)
{
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_;
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;
p.wgt_ = 0.;
p.wgt_last_ = 0.;
p.alive_ = false;
}
}
}

100
src/physics_mg.cpp
View file

@ -22,29 +22,29 @@
namespace openmc {
void
collision_mg(Particle* p)
collision_mg(Particle& p)
{
// Add to the collision counter for the particle
p->n_collision_++;
p.n_collision_++;
// Sample the reaction type
sample_reaction(p);
// Display information about collision
if ((settings::verbosity >= 10) || p->trace_) {
write_message(fmt::format(" Energy Group = {}", p->g_), 1);
if ((settings::verbosity >= 10) || p.trace_) {
write_message(fmt::format(" Energy Group = {}", p.g_), 1);
}
}
void
sample_reaction(Particle* p)
sample_reaction(Particle& p)
{
// Create fission bank sites. Note that while a fission reaction is sampled,
// it never actually "happens", i.e. the weight of the particle does not
// change when sampling fission sites. The following block handles all
// absorption (including fission)
if (model::materials[p->material_]->fissionable_) {
if (model::materials[p.material_]->fissionable_) {
if (settings::run_mode == RunMode::EIGENVALUE ||
(settings::run_mode == RunMode::FIXED_SOURCE &&
settings::create_fission_neutrons)) {
@ -54,12 +54,12 @@ sample_reaction(Particle* p)
// If survival biasing is being used, the following subroutine adjusts the
// weight of the particle. Otherwise, it checks to see if absorption occurs.
if (p->macro_xs_.absorption > 0.) {
if (p.macro_xs_.absorption > 0.) {
absorption(p);
} else {
p->wgt_absorb_ = 0.;
p.wgt_absorb_ = 0.;
}
if (!p->alive_) return;
if (!p.alive_) return;
// Sample a scattering event to determine the energy of the exiting neutron
scatter(p);
@ -67,40 +67,40 @@ 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.alive_) return;
}
}
void
scatter(Particle* p)
scatter(Particle& p)
{
data::mg.macro_xs_[p->material_].sample_scatter(p->g_last_, p->g_, p->mu_,
p->wgt_, p->current_seed());
data::mg.macro_xs_[p.material_].sample_scatter(p.g_last_, p.g_, p.mu_,
p.wgt_, p.current_seed());
// Rotate the angle
p->u() = rotate_angle(p->u(), p->mu_, nullptr, p->current_seed());
p.u() = rotate_angle(p.u(), p.mu_, nullptr, p.current_seed());
// Update energy value for downstream compatability (in tallying)
p->E_ = data::mg.energy_bin_avg_[p->g_];
p.E_ = data::mg.energy_bin_avg_[p.g_];
// Set event component
p->event_ = TallyEvent::SCATTER;
p.event_ = TallyEvent::SCATTER;
}
void
create_fission_sites(Particle* p)
create_fission_sites(Particle& p)
{
// If uniform fission source weighting is turned on, we increase or decrease
// the expected number of fission sites produced
double weight = settings::ufs_on ? ufs_get_weight(*p) : 1.0;
double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0;
// Determine the expected number of neutrons produced
double nu_t = p->wgt_ / simulation::keff * weight *
p->macro_xs_.nu_fission / p->macro_xs_.total;
double nu_t = p.wgt_ / simulation::keff * weight *
p.macro_xs_.nu_fission / p.macro_xs_.total;
// Sample the number of neutrons produced
int nu = static_cast<int>(nu_t);
if (prn(p->current_seed()) <= (nu_t - int(nu_t))) {
if (prn(p.current_seed()) <= (nu_t - int(nu_t))) {
nu++;
}
@ -113,9 +113,9 @@ create_fission_sites(Particle* p)
double nu_d[MAX_DELAYED_GROUPS] = {0.};
// Clear out particle's nu fission bank
p->nu_bank_.clear();
p.nu_bank_.clear();
p->fission_ = true;
p.fission_ = true;
int skipped = 0;
// Determine whether to place fission sites into the shared fission bank
@ -125,18 +125,18 @@ create_fission_sites(Particle* p)
for (int i = 0; i < nu; ++i) {
// Initialize fission site object with particle data
Particle::Bank site;
site.r = p->r();
site.r = p.r();
site.particle = Particle::Type::neutron;
site.wgt = 1. / weight;
site.parent_id = p->id_;
site.progeny_id = p->n_progeny_++;
site.parent_id = p.id_;
site.progeny_id = p.n_progeny_++;
// Sample the cosine of the angle, assuming fission neutrons are emitted
// isotropically
double mu = 2.*prn(p->current_seed()) - 1.;
double mu = 2.*prn(p.current_seed()) - 1.;
// Sample the azimuthal angle uniformly in [0, 2.pi)
double phi = 2. * PI * prn(p->current_seed() );
double phi = 2. * PI * prn(p.current_seed() );
site.u.x = mu;
site.u.y = std::sqrt(1. - mu * mu) * std::cos(phi);
site.u.z = std::sqrt(1. - mu * mu) * std::sin(phi);
@ -144,8 +144,8 @@ create_fission_sites(Particle* p)
// Sample secondary energy distribution for the fission reaction
int dg;
int gout;
data::mg.macro_xs_[p->material_].sample_fission_energy(p->g_, dg, gout,
p->current_seed());
data::mg.macro_xs_[p.material_].sample_fission_energy(p.g_, dg, gout,
p.current_seed());
// Store the energy and delayed groups on the fission bank
site.E = gout;
@ -164,21 +164,21 @@ create_fission_sites(Particle* p)
break;
}
} else {
p->secondary_bank_.push_back(site);
p.secondary_bank_.push_back(site);
}
// Set the delayed group on the particle as well
p->delayed_group_ = dg + 1;
p.delayed_group_ = dg + 1;
// Increment the number of neutrons born delayed
if (p->delayed_group_ > 0) {
if (p.delayed_group_ > 0) {
nu_d[dg]++;
}
// Write fission particles to nuBank
if (use_fission_bank) {
p->nu_bank_.emplace_back();
Particle::NuBank* nu_bank_entry = &p->nu_bank_.back();
p.nu_bank_.emplace_back();
Particle::NuBank* nu_bank_entry = &p.nu_bank_.back();
nu_bank_entry->wgt = site.wgt;
nu_bank_entry->E = site.E;
nu_bank_entry->delayed_group = site.delayed_group;
@ -188,7 +188,7 @@ create_fission_sites(Particle* p)
// If shared fission bank was full, and no fissions could be added,
// set the particle fission flag to false.
if (nu == skipped) {
p->fission_ = false;
p.fission_ = false;
return;
}
@ -197,33 +197,33 @@ create_fission_sites(Particle* p)
nu -= skipped;
// Store the total weight banked for analog fission tallies
p->n_bank_ = nu;
p->wgt_bank_ = nu / weight;
p.n_bank_ = nu;
p.wgt_bank_ = nu / weight;
for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) {
p->n_delayed_bank_[d] = nu_d[d];
p.n_delayed_bank_[d] = nu_d[d];
}
}
void
absorption(Particle* p)
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;
p.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_ -= p.wgt_absorb_;
p.wgt_last_ = p.wgt_;
// Score implicit absorpion estimate of keff
p->keff_tally_absorption_ += p->wgt_absorb_ * p->macro_xs_.nu_fission /
p->macro_xs_.absorption;
p.keff_tally_absorption_ += p.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_ += p->wgt_ * p->macro_xs_.nu_fission /
p->macro_xs_.absorption;
p->alive_ = false;
p->event_ = TallyEvent::ABSORB;
if (p.macro_xs_.absorption > prn(p.current_seed()) * p.macro_xs_.total) {
p.keff_tally_absorption_ += p.wgt_ * p.macro_xs_.nu_fission /
p.macro_xs_.absorption;
p.alive_ = false;
p.event_ = TallyEvent::ABSORB;
}
}