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

Weight non-fission energy deposition by keff

Browse source
This commit is contained in:
Paul Romano 2020-06-29 10:20:20 -05:00
parent f81f5b164f
commit ece6b8b15f
8 changed files with 103 additions and 46 deletions
Download patch file Download diff file Expand all files Collapse all files

2
include/openmc/nuclide.h
View file

@ -84,6 +84,8 @@ public:
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};

3
include/openmc/particle.h
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
//

2
src/bremsstrahlung.cpp
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;
}
}

5
src/nuclide.cpp
View file

@ -267,6 +267,11 @@ Nuclide::Nuclide(hid_t group, const std::vector<double>& 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");

4
src/particle.cpp
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_;

6
src/photon.cpp
View file

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

31
src/physics.cpp
View file

@ -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<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
@ -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<int>(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);
}
}

96
src/tallies/tally_scoring.cpp
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;
}
}
}