Merge pull request #1599 from paulromano/energy-dep-fixes

Two fixes for energy deposition calculations
This commit is contained in:
Amanda Lund 2020-07-04 18:18:50 -05:00 committed by GitHub
commit 557385c2bf
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
14 changed files with 173 additions and 93 deletions

View file

@ -12,19 +12,18 @@ is deposited for a specific reaction is referred to as
"heating numbers" and can be computed using a program like NJOY with the
``heatr`` module.
These heating rate is the product of reaction-specific coefficients and
a reaction cross section
The heating rate is the product of reaction-specific coefficients and a reaction
cross section
.. math::
H(E) = \phi(E)\sum_i\rho_i\sum_rk_{i, r}(E),
and has units energy per time, typically eV / s.
Here, :math:`k_{i, r}` are the KERMA (Kinetic Energy Release in Materials)
[Mack97]_ coefficients for reaction :math:`r` of isotope :math:`i`.
The KERMA coefficients have units energy :math:`\times` cross-section, e.g.
eV-barn, and can be used much like a reaction cross section for the purpose
of tallying energy deposition.
and has units energy per time, typically eV/s. Here, :math:`k_{i, r}` are the
KERMA (Kinetic Energy Release in Materials) [Mack97]_ coefficients for reaction
:math:`r` of isotope :math:`i`. The KERMA coefficients have units of energy
:math:`\times` cross-section (e.g., eV-barn) and can be used much like a reaction
cross section for the purpose of tallying energy deposition.
KERMA coefficients can be computed using the energy-balance method with
a nuclear data processing code like NJOY, which performs the following
@ -56,11 +55,11 @@ broken up into the following categories:
- :math:`E_{\beta}` - energy of released :math:`\beta` particles
- :math:`E_{\nu}` - energy of neutrinos
These components are defined in MF=1,MT=458 data in a standard ENDF/B-6 formatted
file. All these quantities may depend upon incident neutron energy,
but this dependence is not shown to make the following demonstrations cleaner.
As neutrinos scarcely interact with matter, the recoverable energy from
fission is defined as
These components are defined in MF=1, MT=458 data in a standard ENDF-6 formatted
file. All these quantities may depend upon incident neutron energy, but this
dependence is not shown to make the following demonstrations cleaner. As
neutrinos scarcely interact with matter, the recoverable energy from fission is
defined as
.. math::

View file

@ -149,6 +149,16 @@ Geometry and Visualization
Miscellaneous
-------------
- Shikhar Kumar, Benoit Forget, and Kord Smith, "`Stationarity Diagnostic using
Functional Expansion Tallies
<https://doi.org/10.1016/j.anucene.2020.107388>`_", *Ann. Nucl. Energy*,
**143**, 107388 (2020).
- T. Eade, B. Colling, J. Naish, L. W. Packer, and A. Valentine, "`Shutdown dose
rate benchmarking using modern particle transport codes
<https://doi.org/10.1088/1741-4326/ab8181>`_, *Nucl. Fusion*, **60**, 056024
(2020).
- Jiankai Yu, Qiudong Wang, Ding She, and Benoit Forget, "Modelling of the
HTR-PM Pebble-bed Reactor using OpenMC", *Trans. Am. Nucl. Soc.*, **122**,
643-646 (2020).

View file

@ -82,6 +82,10 @@ public:
std::unique_ptr<Function1D> total_nu_; //!< Total neutron yield
std::unique_ptr<Function1D> fission_q_prompt_; //!< Prompt fission energy release
std::unique_ptr<Function1D> fission_q_recov_; //!< Recoverable fission energy release
std::unique_ptr<Function1D> prompt_photons_; //!< Prompt photon energy release
std::unique_ptr<Function1D> delayed_photons_; //!< Delayed photon energy release
std::unique_ptr<Function1D> fragments_; //!< Fission fragment energy release
std::unique_ptr<Function1D> betas_; //!< Delayed beta energy release
// Resonance scattering information
bool resonant_ {false};

View file

@ -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
//

View file

@ -120,8 +120,8 @@ namespace data {
extern xt::xtensor<double, 1> compton_profile_pz; //! Compton profile momentum grid
//! Photon interaction data for each element
extern std::vector<std::unique_ptr<PhotonInteraction>> elements;
extern std::unordered_map<std::string, int> element_map;
extern std::vector<std::unique_ptr<PhotonInteraction>> elements;
} // namespace data

View file

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

View file

@ -224,13 +224,11 @@ double ContinuousTabular::sample(double E, uint64_t* seed) const
double E_l_k = distribution_[l].e_out[k];
double p_l_k = distribution_[l].p[k];
double E_out;
double E_out = E_l_k;
if (distribution_[l].interpolation == Interpolation::histogram) {
// Histogram interpolation
if (p_l_k > 0.0 && k >= n_discrete) {
E_out = E_l_k + (r1 - c_k)/p_l_k;
} else {
E_out = E_l_k;
}
} else if (distribution_[l].interpolation == Interpolation::lin_lin) {
@ -238,12 +236,14 @@ double ContinuousTabular::sample(double E, uint64_t* seed) const
double E_l_k1 = distribution_[l].e_out[k+1];
double p_l_k1 = distribution_[l].p[k+1];
double frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k);
if (frac == 0.0) {
E_out = E_l_k + (r1 - c_k)/p_l_k;
} else {
E_out = E_l_k + (std::sqrt(std::max(0.0, p_l_k*p_l_k +
2.0*frac*(r1 - c_k))) - p_l_k)/frac;
if (E_l_k != E_l_k1) {
double frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k);
if (frac == 0.0) {
E_out = E_l_k + (r1 - c_k)/p_l_k;
} else {
E_out = E_l_k + (std::sqrt(std::max(0.0, p_l_k*p_l_k +
2.0*frac*(r1 - c_k))) - p_l_k)/frac;
}
}
} else {
throw std::runtime_error{"Unexpected interpolation for continuous energy "

View file

@ -262,21 +262,24 @@ Nuclide::Nuclide(hid_t group, const std::vector<double>& temperature)
}
// Read fission energy release data if present
std::unique_ptr<Function1D> prompt_photons;
std::unique_ptr<Function1D> delayed_photons;
if (object_exists(group, "fission_energy_release")) {
hid_t fer_group = open_group(group, "fission_energy_release");
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");
delayed_photons = read_function(fer_group, "delayed_photons");
prompt_photons_ = read_function(fer_group, "prompt_photons");
delayed_photons_ = read_function(fer_group, "delayed_photons");
close_group(fer_group);
}
this->create_derived(prompt_photons.get(), delayed_photons.get());
this->create_derived(prompt_photons_.get(), delayed_photons_.get());
}
Nuclide::~Nuclide()

View file

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

View file

@ -27,8 +27,8 @@ namespace data {
xt::xtensor<double, 1> compton_profile_pz;
std::vector<std::unique_ptr<PhotonInteraction>> elements;
std::unordered_map<std::string, int> element_map;
std::vector<std::unique_ptr<PhotonInteraction>> elements;
} // namespace data
@ -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

View file

@ -4,6 +4,7 @@
#include "openmc/bremsstrahlung.h"
#include "openmc/constants.h"
#include "openmc/eigenvalue.h"
#include "openmc/endf.h"
#include "openmc/error.h"
#include "openmc/material.h"
#include "openmc/math_functions.h"
@ -302,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
@ -363,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
@ -388,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;
@ -433,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;
@ -567,8 +568,21 @@ 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) {
// For fission, artificially increase the photon yield to account
// for delayed photons
double f = 1.0;
if (settings::delayed_photon_scaling) {
if (is_fission(rx->mt_)) {
if (nuc->prompt_photons_ && nuc->delayed_photons_) {
double energy_prompt = (*nuc->prompt_photons_)(p.E_);
double energy_delayed = (*nuc->delayed_photons_)(p.E_);
f = (energy_prompt + energy_delayed)/(energy_prompt);
}
}
}
// add to cumulative probability
prob += (*rx->products_[j].yield_)(p.E_) * xs;
prob += f * (*rx->products_[j].yield_)(p.E_) * xs;
*i_rx = i;
*i_product = j;
@ -1126,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<int>(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
@ -1137,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<int>(y_t);
if (prn(p.current_seed()) <= y_t - y) ++y;
@ -1158,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);
}
}

View file

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

View file

@ -20,50 +20,50 @@ tally 1:
tally 2:
2.602512E+00
2.258176E+00
4.394176E+08
6.437709E+16
4.185969E+08
5.842526E+16
0.000000E+00
0.000000E+00
2.602512E+00
2.258176E+00
4.394176E+08
6.437709E+16
4.185969E+08
5.842526E+16
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.921131E+06
2.853657E+12
2.366791E+06
1.872101E+12
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.921131E+06
2.853657E+12
2.366791E+06
1.872101E+12
0.000000E+00
0.000000E+00
tally 3:
2.663535E+00
2.364845E+00
4.394176E+08
6.437709E+16
4.185969E+08
5.842526E+16
0.000000E+00
0.000000E+00
2.663535E+00
2.364845E+00
4.394176E+08
6.437709E+16
4.185969E+08
5.842526E+16
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.780823E+06
2.586737E+12
2.185301E+06
1.596625E+12
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.921131E+06
2.853657E+12
2.366791E+06
1.872101E+12
0.000000E+00
0.000000E+00

View file

@ -1 +1 @@
4d9e4a6d4891d02fac2222e0c8665033f0c955ac0082e0c6851c7674a6e967721bad8b77bf0d6c243736128a77a274b4664af54ab10462416ef6c4d18e7218ca
ba4369130ac65812939cb417c7088ae6c8cdec118fd7dce5645b40a20a3ae47072cd72c24aa091ccc09f833a5bfc0bd72407a0dcde7d3120e39379942614a75b