From 43eae72aeb33ff8cc9c6148e33ebec3efa76de90 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 23 Apr 2020 10:11:43 -0500 Subject: [PATCH] Use references consistently in physics modules --- include/openmc/physics.h | 32 +-- include/openmc/physics_common.h | 3 +- include/openmc/physics_mg.h | 10 +- src/particle.cpp | 4 +- src/physics.cpp | 388 ++++++++++++++++---------------- src/physics_common.cpp | 16 +- src/physics_mg.cpp | 100 ++++---- 7 files changed, 276 insertions(+), 277 deletions(-) diff --git a/include/openmc/physics.h b/include/openmc/physics.h index ac83dbc70..670a90bfb 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -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 diff --git a/include/openmc/physics_common.h b/include/openmc/physics_common.h index 33fb8ed61..1398557dc 100644 --- a/include/openmc/physics_common.h +++ b/include/openmc/physics_common.h @@ -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 diff --git a/include/openmc/physics_mg.h b/include/openmc/physics_mg.h index 5cb7226d9..c73bd2eb5 100644 --- a/include/openmc/physics_mg.h +++ b/include/openmc/physics_mg.h @@ -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 diff --git a/src/particle.cpp b/src/particle.cpp index f7c7f3c0e..8008310fd 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -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 diff --git a/src/physics.cpp b/src/physics.cpp index 068fff7ce..ef28f8cdb 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -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(p->type_)) { + switch (static_cast(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(p->type_); - if (p->E_ < settings::energy_cutoff[type]) { - p->alive_ = false; - p->wgt_ = 0.0; + int type = static_cast(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(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(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(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(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(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(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); } } diff --git a/src/physics_common.cpp b/src/physics_common.cpp index 1ea40a896..32e00056d 100644 --- a/src/physics_common.cpp +++ b/src/physics_common.cpp @@ -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; } } } diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 254379d59..b8e8fea45 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -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(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; } }