diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index f3c811b51..3ea3a05a7 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -84,6 +84,8 @@ public: std::unique_ptr fission_q_recov_; //!< Recoverable fission energy release std::unique_ptr prompt_photons_; //!< Prompt photon energy release std::unique_ptr delayed_photons_; //!< Delayed photon energy release + std::unique_ptr fragments_; //!< Fission fragment energy release + std::unique_ptr betas_; //!< Delayed beta energy release // Resonance scattering information bool resonant_ {false}; diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 85641f7c4..e551261a7 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -212,10 +212,11 @@ public: // //! stores the current phase space attributes of the particle in the //! secondary bank and increments the number of sites in the secondary bank. + //! \param wgt Weight of the secondary particle //! \param u Direction of the secondary particle //! \param E Energy of the secondary particle in [eV] //! \param type Particle type - void create_secondary(Direction u, double E, Type type); + void create_secondary(double wgt, Direction u, double E, Type type); //! initialize from a source site // diff --git a/src/bremsstrahlung.cpp b/src/bremsstrahlung.cpp index e7476124d..a7fae50c1 100644 --- a/src/bremsstrahlung.cpp +++ b/src/bremsstrahlung.cpp @@ -108,7 +108,7 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost) if (w > settings::energy_cutoff[photon]) { // Create secondary photon - p.create_secondary(p.u(), w, Particle::Type::photon); + p.create_secondary(p.wgt_, p.u(), w, Particle::Type::photon); *E_lost += w; } } diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 8e199cb12..9b4500e63 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -267,6 +267,11 @@ Nuclide::Nuclide(hid_t group, const std::vector& temperature) fission_q_prompt_ = read_function(fer_group, "q_prompt"); fission_q_recov_ = read_function(fer_group, "q_recoverable"); + // Read fission fragment and delayed beta energy release. This is needed for + // energy normalization in k-eigenvalue calculations + fragments_ = read_function(fer_group, "fragments"); + betas_ = read_function(fer_group, "betas"); + // We need prompt/delayed photon energy release for scaling fission photon // production prompt_photons_ = read_function(fer_group, "prompt_photons"); diff --git a/src/particle.cpp b/src/particle.cpp index c4458a1e4..ae75d91b4 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -84,13 +84,13 @@ Particle::clear() } void -Particle::create_secondary(Direction u, double E, Type type) +Particle::create_secondary(double wgt, Direction u, double E, Type type) { secondary_bank_.emplace_back(); auto& bank {secondary_bank_.back()}; bank.particle = type; - bank.wgt = wgt_; + bank.wgt = wgt; bank.r = this->r(); bank.u = u; bank.E = settings::run_CE ? E : g_; diff --git a/src/photon.cpp b/src/photon.cpp index c3cd46010..6333f9b29 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -669,7 +669,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi); u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi); double E = shell.binding_energy; - p.create_secondary(u, E, Particle::Type::photon); + p.create_secondary(p.wgt_, u, E, Particle::Type::photon); return; } @@ -701,7 +701,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl // Non-radiative transition -- Auger/Coster-Kronig effect // Create auger electron - p.create_secondary(u, E, Particle::Type::electron); + p.create_secondary(p.wgt_, u, E, Particle::Type::electron); // Fill hole left by emitted auger electron int i_hole = shell_map_.at(secondary); @@ -711,7 +711,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl // Radiative transition -- get X-ray energy // Create fluorescent photon - p.create_secondary(u, E, Particle::Type::photon); + p.create_secondary(p.wgt_, u, E, Particle::Type::photon); } // Fill hole created by electron transitioning to the photoelectron hole diff --git a/src/physics.cpp b/src/physics.cpp index e3ede4d61..2636207fb 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -303,7 +303,7 @@ void sample_photon_reaction(Particle& p) 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); + p.create_secondary(p.wgt_, u, E_electron, Particle::Type::electron); } // TODO: Compton subshell data does not match atomic relaxation data @@ -364,7 +364,7 @@ void sample_photon_reaction(Particle& p) 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(p.wgt_, u, E_electron, Particle::Type::electron); // Allow electrons to fill orbital and produce auger electrons // and fluorescent photons @@ -389,11 +389,11 @@ void sample_photon_reaction(Particle& p) // Create secondary electron Direction u = rotate_angle(p.u(), mu_electron, nullptr, p.current_seed()); - p.create_secondary(u, E_electron, Particle::Type::electron); + p.create_secondary(p.wgt_, 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); + p.create_secondary(p.wgt_, u, E_positron, Particle::Type::positron); p.event_ = TallyEvent::ABSORB; p.event_mt_ = PAIR_PROD; @@ -434,8 +434,8 @@ void sample_positron_reaction(Particle& p) 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(p.wgt_, u, MASS_ELECTRON_EV, Particle::Type::photon); + p.create_secondary(p.wgt_, -u, MASS_ELECTRON_EV, Particle::Type::photon); p.E_ = 0.0; p.alive_ = false; @@ -1140,7 +1140,7 @@ void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p) 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.wgt_, p.u(), p.E_, Particle::Type::neutron); } } else { // Otherwise, change weight of particle based on yield @@ -1151,7 +1151,7 @@ void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p) void sample_secondary_photons(Particle& p, int i_nuclide) { // Sample the number of photons produced - double y_t = p.neutron_xs_[i_nuclide].photon_prod / + double y_t = 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; @@ -1172,8 +1172,21 @@ void sample_secondary_photons(Particle& p, int i_nuclide) // Sample the new direction Direction u = rotate_angle(p.u(), mu, nullptr, p.current_seed()); + // In a k-eigenvalue simulation, it's necessary to provide higher weight to + // secondary photons from non-fission reactions to properly balance energy + // release and deposition. See D. P. Griesheimer, S. J. Douglass, and M. H. + // Stedry, "Self-consistent energy normalization for quasistatic reactor + // calculations", Proc. PHYSOR, Cambridge, UK, Mar 29-Apr 2, 2020. + double wgt; + if (settings::run_mode == RunMode::EIGENVALUE && !is_fission(rx->mt_)) { + wgt = simulation::keff * p.wgt_; + } else { + wgt = p.wgt_; + } + // Create the secondary photon - p.create_secondary(u, E, Particle::Type::photon); + p.create_secondary(wgt, u, E, Particle::Type::photon); + } } diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 3e551ff9e..049c7be6b 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -259,9 +259,27 @@ double get_nuclide_neutron_heating(const Particle& p, const Nuclide& nuc, auto i_grid = p.neutron_xs_[i_nuclide].index_grid; if (i_grid < xs.threshold) return 0.0; + // Determine total kerma auto f = p.neutron_xs_[i_nuclide].interp_factor; - return (1.0 - f) * xs.value[i_grid-xs.threshold] + double kerma = (1.0 - f) * xs.value[i_grid-xs.threshold] + f * xs.value[i_grid-xs.threshold+1]; + + if (settings::run_mode == RunMode::EIGENVALUE) { + // Determine kerma for fission as (EFR + EB)*sigma_f + double kerma_fission = nuc.fragments_ ? + ((*nuc.fragments_)(p.E_last_) + (*nuc.betas_)(p.E_last_)) + *p.neutron_xs_[i_nuclide].fission : 0.0; + + // Determine non-fission kerma as difference + double kerma_non_fission = kerma - kerma_fission; + + // Re-weight non-fission kerma by keff to properly balance energy release + // and deposition. See D. P. Griesheimer, S. J. Douglass, and M. H. Stedry, + // "Self-consistent energy normalization for quasistatic reactor + // calculations", Proc. PHYSOR, Cambridge, UK, Mar 29-Apr 2, 2020. + kerma = simulation::keff * kerma_non_fission + kerma_fission; + } + return kerma; } //! Helper function to obtain neutron heating [eV] @@ -452,6 +470,51 @@ score_fission_eout(Particle& p, int i_tally, int i_score, int score_bin) p.filter_matches_[i_eout_filt].bins_[i_bin] = bin_energyout; } +double get_nuclide_xs(const Particle& p, int i_nuclide, int score_bin) { + const auto& nuc {*data::nuclides[i_nuclide]}; + + // Get reaction object, or return 0 if reaction is not present + auto m = nuc.reaction_index_[score_bin]; + if (m == C_NONE) return 0.0; + const auto& rxn {*nuc.reactions_[m]}; + + auto i_temp = p.neutron_xs_[i_nuclide].index_temp; + if (i_temp >= 0) { // Can be false due to multipole + // Get index on energy grid and interpolation factor + auto i_grid = p.neutron_xs_[i_nuclide].index_grid; + auto f = p.neutron_xs_[i_nuclide].interp_factor; + + // Calculate interpolated cross section + const auto& xs {rxn.xs_[i_temp]}; + double value; + if (i_grid >= xs.threshold) { + value = ((1.0 - f) * xs.value[i_grid-xs.threshold] + + f * xs.value[i_grid-xs.threshold+1]); + } else { + value = 0.0; + } + + if (settings::run_mode == RunMode::EIGENVALUE && score_bin == HEATING_LOCAL) { + // Determine kerma for fission as (EFR + EGP + EGD + EB)*sigma_f + double kerma_fission = nuc.fragments_ ? + ((*nuc.fragments_)(p.E_last_) + (*nuc.betas_)(p.E_last_) + + (*nuc.prompt_photons_)(p.E_last_) + (*nuc.delayed_photons_)(p.E_last_)) + * p.neutron_xs_[i_nuclide].fission : 0.0; + + // Determine non-fission kerma as difference + double kerma_non_fission = value - kerma_fission; + + // Re-weight non-fission kerma by keff to properly balance energy release + // and deposition. See D. P. Griesheimer, S. J. Douglass, and M. H. Stedry, + // "Self-consistent energy normalization for quasistatic reactor + // calculations", Proc. PHYSOR, Cambridge, UK, Mar 29-Apr 2, 2020. + value = simulation::keff * kerma_non_fission + kerma_fission; + } + return value; + } + return 0.0; +} + //! Update tally results for continuous-energy tallies with any estimator. // //! For analog tallies, the flux estimate depends on the score type so the flux @@ -1252,40 +1315,13 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, + std::to_string(tally.id_)); score = 0.; if (i_nuclide >= 0) { - const auto& nuc {*data::nuclides[i_nuclide]}; - auto m = nuc.reaction_index_[score_bin]; - if (m == C_NONE) continue; - const auto& rxn {*nuc.reactions_[m]}; - auto i_temp = p.neutron_xs_[i_nuclide].index_temp; - if (i_temp >= 0) { // Can be false due to multipole - auto i_grid = p.neutron_xs_[i_nuclide].index_grid; - auto f = p.neutron_xs_[i_nuclide].interp_factor; - const auto& xs {rxn.xs_[i_temp]}; - if (i_grid >= xs.threshold) { - score = ((1.0 - f) * xs.value[i_grid-xs.threshold] - + f * xs.value[i_grid-xs.threshold+1]) * atom_density * flux; - } - } + score = get_nuclide_xs(p, i_nuclide, score_bin) * atom_density * flux; } else if (p.material_ != MATERIAL_VOID) { const Material& material {*model::materials[p.material_]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); - const auto& nuc {*data::nuclides[j_nuclide]}; - auto m = nuc.reaction_index_[score_bin]; - if (m == C_NONE) continue; - const auto& rxn {*nuc.reactions_[m]}; - auto i_temp = p.neutron_xs_[j_nuclide].index_temp; - if (i_temp >= 0) { // Can be false due to multipole - auto i_grid = p.neutron_xs_[j_nuclide].index_grid; - auto f = p.neutron_xs_[j_nuclide].interp_factor; - const auto& xs {rxn.xs_[i_temp]}; - if (i_grid >= xs.threshold) { - score += ((1.0 - f) * xs.value[i_grid-xs.threshold] - + f * xs.value[i_grid-xs.threshold+1]) * atom_density - * flux; - } - } + score += get_nuclide_xs(p, j_nuclide, score_bin) * atom_density * flux; } } }