From 6ce85a9f0ce67a858c8a5234ad3e88ccab905a2c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 25 Feb 2019 09:55:28 -0600 Subject: [PATCH] Fix off-by-one for index_grid --- src/nuclide.cpp | 13 +++++-------- src/physics.cpp | 9 ++++----- src/reaction.cpp | 2 ++ src/tallies/tally_scoring.cpp | 16 +++++++--------- 4 files changed, 18 insertions(+), 22 deletions(-) diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 592bcceb5..a279abf2b 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -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]; diff --git a/src/physics.cpp b/src/physics.cpp index 8b65fe852..b93c2c232 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -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 diff --git a/src/reaction.cpp b/src/reaction.cpp index 66b841ef4..0fcac2e1d 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -42,6 +42,8 @@ Reaction::Reaction(hid_t group, const std::vector& 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); diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index cdd80adba..a0c2883b0 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -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; } }