From ba1a39e28175e74af849298d77fa77427cfc9a23 Mon Sep 17 00:00:00 2001 From: amandalund Date: Sun, 3 Oct 2021 18:35:34 -0500 Subject: [PATCH] Make atomic_relaxation iterative instead of recursive --- include/openmc/photon.h | 5 ++- src/photon.cpp | 91 +++++++++++++++++++++++------------------ src/physics.cpp | 9 ++-- 3 files changed, 60 insertions(+), 45 deletions(-) diff --git a/include/openmc/photon.h b/include/openmc/photon.h index 4dc4068b52..f5a3d8e81b 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -57,7 +57,7 @@ public: void pair_production(double alpha, double* E_electron, double* E_positron, double* mu_electron, double* mu_positron, uint64_t* seed) const; - void atomic_relaxation(const ElectronSubshell& shell, Particle& p) const; + void atomic_relaxation(int i_shell, Particle& p) const; // Data members std::string name_; //!< Name of element, e.g. "Zr" @@ -98,6 +98,9 @@ public: // Bremsstrahlung scaled DCS xt::xtensor dcs_; + // Constant data + static constexpr int MAX_STACK_SIZE = + 7; //!< maximum possible size of atomic relaxation stack private: void compton_doppler( double alpha, double mu, double* E_out, int* i_shell, uint64_t* seed) const; diff --git a/src/photon.cpp b/src/photon.cpp index 18f79dd0c1..b1720019f3 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -684,49 +684,60 @@ void PhotonInteraction::pair_production(double alpha, double* E_electron, *mu_positron = (rn + beta) / (rn * beta + 1.0); } -void PhotonInteraction::atomic_relaxation( - const ElectronSubshell& shell, Particle& p) const +void PhotonInteraction::atomic_relaxation(int i_shell, Particle& p) const { - // If no transitions, assume fluorescent photon from captured free electron - if (shell.transitions.empty()) { + // Stack for unprocessed holes left by transitioning electrons + int n_holes = 0; + array holes; + + // Push the initial hole onto the stack + holes[n_holes++] = i_shell; + + while (n_holes > 0) { + // Pop the next hole off the stack + int i_hole = holes[--n_holes]; + const auto& shell {shells_[i_hole]}; + + // If no transitions, assume fluorescent photon from captured free electron + if (shell.transitions.empty()) { + Direction u = isotropic_direction(p.current_seed()); + double E = shell.binding_energy; + p.create_secondary(p.wgt(), u, E, ParticleType::photon); + continue; + } + + // Sample transition + double c = -prn(p.current_seed()); + int i_trans; + for (i_trans = 0; i_trans < shell.transitions.size(); ++i_trans) { + c += shell.transitions[i_trans].probability; + if (c > 0) + break; + } + const auto& transition = shell.transitions[i_trans]; + + // Sample angle isotropically Direction u = isotropic_direction(p.current_seed()); - double E = shell.binding_energy; - p.create_secondary(p.wgt(), u, E, ParticleType::photon); - return; + + // Push the hole created by the electron transitioning to the photoelectron + // hole onto the stack + holes[n_holes++] = transition.primary_subshell; + + if (transition.secondary_subshell != -1) { + // Non-radiative transition -- Auger/Coster-Kronig effect + + // Push the hole left by emitted auger electron onto the stack + holes[n_holes++] = transition.secondary_subshell; + + // Create auger electron + p.create_secondary(p.wgt(), u, transition.energy, ParticleType::electron); + } else { + // Radiative transition -- get X-ray energy + + // Create fluorescent photon + p.create_secondary(p.wgt(), u, transition.energy, ParticleType::photon); + } } - - // Sample transition - double c = -prn(p.current_seed()); - int i_trans; - for (i_trans = 0; i_trans < shell.transitions.size(); ++i_trans) { - c += shell.transitions[i_trans].probability; - if (c > 0) - break; - } - const auto& transition = shell.transitions[i_trans]; - - // Sample angle isotropically - Direction u = isotropic_direction(p.current_seed()); - - if (transition.secondary_subshell != -1) { - // Non-radiative transition -- Auger/Coster-Kronig effect - - // Create auger electron - p.create_secondary(p.wgt(), u, transition.energy, ParticleType::electron); - - // Fill hole left by emitted auger electron - const auto& hole = shells_[transition.secondary_subshell]; - this->atomic_relaxation(hole, p); - } else { - // Radiative transition -- get X-ray energy - - // Create fluorescent photon - p.create_secondary(p.wgt(), u, transition.energy, ParticleType::photon); - } - - // Fill hole created by electron transitioning to the photoelectron hole - const auto& hole = shells_[transition.primary_subshell]; - this->atomic_relaxation(hole, p); } //============================================================================== diff --git a/src/physics.cpp b/src/physics.cpp index a466c0f9d2..3a52396e1e 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -326,8 +326,7 @@ void sample_photon_reaction(Particle& p) // Allow electrons to fill orbital and produce auger electrons // and fluorescent photons if (i_shell >= 0) { - const auto& shell = element.shells_[i_shell]; - element.atomic_relaxation(shell, p); + element.atomic_relaxation(i_shell, p); } phi += PI; @@ -341,7 +340,9 @@ void sample_photon_reaction(Particle& p) // Photoelectric effect double prob_after = prob + micro.photoelectric; if (prob_after > cutoff) { - for (const auto& shell : element.shells_) { + for (int i_shell = 0; i_shell < element.shells_.size(); ++i_shell) { + const auto& shell {element.shells_[i_shell]}; + // Get grid index and interpolation factor int i_grid = micro.index_grid; double f = micro.interp_factor; @@ -387,7 +388,7 @@ void sample_photon_reaction(Particle& p) // Allow electrons to fill orbital and produce auger electrons // and fluorescent photons - element.atomic_relaxation(shell, p); + element.atomic_relaxation(i_shell, p); p.event() = TallyEvent::ABSORB; p.event_mt() = 533 + shell.index_subshell; p.alive() = false;