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

Fix off-by-one for index_grid

Browse source
This commit is contained in:
Paul Romano 2019-02-25 09:55:28 -06:00
parent 4a89075890
commit 6ce85a9f0c
4 changed files with 18 additions and 22 deletions
Download patch file Download diff file Expand all files Collapse all files

13
src/nuclide.cpp
View file

@ -282,8 +282,7 @@ void Nuclide::create_derived()
reaction_index_[rx->mt_] = i;
for (int t = 0; t < kTs_.size(); ++t) {
// TODO: off-by-one
int j = rx->xs_[t].threshold - 1;
int j = rx->xs_[t].threshold;
int n = rx->xs_[t].value.size();
auto xs = xt::adapt(rx->xs_[t].value);
@ -467,7 +466,7 @@ void Nuclide::calculate_elastic_xs() const
// Get temperature index, grid index, and interpolation factor
auto& micro = simulation::micro_xs[i_nuclide_];
int i_temp = micro.index_temp;
int i_grid = micro.index_grid - 1;
int i_grid = micro.index_grid;
double f = micro.interp_factor;
if (i_temp >= 0) {
@ -549,7 +548,7 @@ void Nuclide::calculate_xs(int i_sab, double E, int i_log_union,
// set to -1 to force a segfault in case a developer messes up and tries
// to use it with multipole.
micro_xs.index_temp = -1;
micro_xs.index_grid = 0;
micro_xs.index_grid = -1;
micro_xs.interp_factor = 0.0;
} else {
@ -613,8 +612,7 @@ void Nuclide::calculate_xs(int i_sab, double E, int i_log_union,
(grid.energy[i_grid + 1]- grid.energy[i_grid]);
micro_xs.index_temp = i_temp;
// TODO: off-by-one
micro_xs.index_grid = i_grid + 1;
micro_xs.index_grid = i_grid;
micro_xs.interp_factor = f;
// Calculate microscopic nuclide total cross section
@ -665,8 +663,7 @@ void Nuclide::calculate_xs(int i_sab, double E, int i_log_union,
continue;
}
// TODO: off-by-one
int threshold = rx->xs_[i_temp].threshold - 1;
int threshold = rx->xs_[i_temp].threshold;
if (i_grid >= threshold) {
micro_xs.reaction[j] = (1.0 - f)*rx_xs[i_grid - threshold] +
f*rx_xs[i_grid - threshold + 1];

9
src/physics.cpp
View file

@ -597,7 +597,7 @@ void scatter(Particle* p, int i_nuclide)
const auto& nuc {data::nuclides[i_nuclide]};
const auto& micro {simulation::micro_xs[i_nuclide]};
int i_temp = micro.index_temp;
int i_grid = micro.index_grid - 1;
int i_grid = micro.index_grid;
double f = micro.interp_factor;
// For tallying purposes, this routine might be called directly. In that
@ -655,12 +655,11 @@ void scatter(Particle* p, int i_nuclide)
// if energy is below threshold for this reaction, skip it
const auto& xs {nuc->reactions_[i]->xs_[i_temp]};
int threshold = xs.threshold - 1;
if (i_grid < threshold) continue;
if (i_grid < xs.threshold) continue;
// add to cumulative probability
prob += (1.0 - f)*xs.value[i_grid - threshold] +
f*xs.value[i_grid - threshold + 1];
prob += (1.0 - f)*xs.value[i_grid - xs.threshold] +
f*xs.value[i_grid - xs.threshold + 1];
}
// Perform collision physics for inelastic scattering

2
src/reaction.cpp
View file

@ -42,6 +42,8 @@ Reaction::Reaction(hid_t group, const std::vector<int>& temperatures)
// Get threshold index
TemperatureXS xs;
read_attribute(dset, "threshold_idx", xs.threshold);
// TODO: change HDF5 format so that threshold_idx is 0-based
--xs.threshold;
// Read cross section values
read_dataset(dset, xs.value);

16
src/tallies/tally_scoring.cpp
View file

@ -1175,13 +1175,12 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
const auto& rxn {*nuc.reactions_[m]};
auto i_temp = simulation::micro_xs[i_nuclide].index_temp;
if (i_temp >= 0) { // Can be false due to multipole
auto i_grid = simulation::micro_xs[i_nuclide].index_grid - 1;
auto i_grid = simulation::micro_xs[i_nuclide].index_grid;
auto f = simulation::micro_xs[i_nuclide].interp_factor;
const auto& xs {rxn.xs_[i_temp]};
auto threshold = xs.threshold - 1;
if (i_grid >= xs.threshold) {
score = ((1.0 - f) * xs.value[i_grid-threshold]
+ f * xs.value[i_grid-threshold+1]) * atom_density * flux;
score = ((1.0 - f) * xs.value[i_grid-xs.threshold]
+ f * xs.value[i_grid-xs.threshold+1]) * atom_density * flux;
}
}
} else {
@ -1196,13 +1195,12 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
const auto& rxn {*nuc.reactions_[m]};
auto i_temp = simulation::micro_xs[j_nuclide].index_temp;
if (i_temp >= 0) { // Can be false due to multipole
auto i_grid = simulation::micro_xs[j_nuclide].index_grid - 1;
auto i_grid = simulation::micro_xs[j_nuclide].index_grid;
auto f = simulation::micro_xs[j_nuclide].interp_factor;
const auto& xs {rxn.xs_[i_temp]};
auto threshold = xs.threshold - 1;
if (i_grid >= threshold) {
score += ((1.0 - f) * xs.value[i_grid-threshold]
+ f * xs.value[i_grid-threshold+1]) * atom_density
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;
}
}