From 9c7dcbba633d6bf74ea8b1135ff7525775a60f01 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Thu, 6 Jan 2022 14:51:38 +0100 Subject: [PATCH] Off by 1 error on Nuclide::cs_cdf Error occured in regression_tests/resonance_scattering/test.py --- src/nuclide.cpp | 2 +- src/physics.cpp | 27 ++++++++++++++------------- 2 files changed, 15 insertions(+), 14 deletions(-) diff --git a/src/nuclide.cpp b/src/nuclide.cpp index ae3085aff..564e27d44 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -455,7 +455,7 @@ void Nuclide::create_derived( xs_cdf_sum += (std::sqrt(E[i]) * xs[i] + std::sqrt(E[i + 1]) * xs[i + 1]) / 2.0 * (E[i + 1] - E[i]); - xs_cdf_[i] = xs_cdf_sum; + xs_cdf_[i+1] = xs_cdf_sum; } } } diff --git a/src/physics.cpp b/src/physics.cpp index 36161f3ea..bb9873bd1 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -936,18 +936,18 @@ Direction sample_target_velocity(const Nuclide& nuc, double E, Direction u, } else if (sampling_method == ResScatMethod::rvs) { // interpolate xs CDF since we're not exactly at the energy indices // cdf value at lower bound attainable energy - double m = (nuc.xs_cdf_[i_E_low] - nuc.xs_cdf_[i_E_low - 1]) / - (nuc.energy_0K_[i_E_low + 1] - nuc.energy_0K_[i_E_low]); - double cdf_low = - nuc.xs_cdf_[i_E_low - 1] + m * (E_low - nuc.energy_0K_[i_E_low]); - if (E_low <= nuc.energy_0K_.front()) - cdf_low = 0.0; + double cdf_low = 0.0; + if (E_low > nuc.energy_0K_.front()) { + double m = (nuc.xs_cdf_[i_E_low + 1] - nuc.xs_cdf_[i_E_low]) / + (nuc.energy_0K_[i_E_low + 1] - nuc.energy_0K_[i_E_low]); + cdf_low = nuc.xs_cdf_[i_E_low] + m * (E_low - nuc.energy_0K_[i_E_low]); + } // cdf value at upper bound attainable energy - m = (nuc.xs_cdf_[i_E_up] - nuc.xs_cdf_[i_E_up - 1]) / + double m = (nuc.xs_cdf_[i_E_up + 1] - nuc.xs_cdf_[i_E_up]) / (nuc.energy_0K_[i_E_up + 1] - nuc.energy_0K_[i_E_up]); double cdf_up = - nuc.xs_cdf_[i_E_up - 1] + m * (E_up - nuc.energy_0K_[i_E_up]); + nuc.xs_cdf_[i_E_up] + m * (E_up - nuc.energy_0K_[i_E_up]); while (true) { // directly sample Maxwellian @@ -955,14 +955,15 @@ Direction sample_target_velocity(const Nuclide& nuc, double E, Direction u, // sample a relative energy using the xs cdf double cdf_rel = cdf_low + prn(seed) * (cdf_up - cdf_low); - int i_E_rel = lower_bound_index( - &nuc.xs_cdf_[i_E_low - 1], &nuc.xs_cdf_[i_E_up + 1], cdf_rel); + int i_E_rel = lower_bound_index(nuc.xs_cdf_.begin() + i_E_low, + nuc.xs_cdf_.begin() + i_E_up+2, + cdf_rel); double E_rel = nuc.energy_0K_[i_E_low + i_E_rel]; - double m = (nuc.xs_cdf_[i_E_low + i_E_rel] - - nuc.xs_cdf_[i_E_low + i_E_rel - 1]) / + double m = (nuc.xs_cdf_[i_E_low + i_E_rel + 1] - + nuc.xs_cdf_[i_E_low + i_E_rel]) / (nuc.energy_0K_[i_E_low + i_E_rel + 1] - nuc.energy_0K_[i_E_low + i_E_rel]); - E_rel += (cdf_rel - nuc.xs_cdf_[i_E_low + i_E_rel - 1]) / m; + E_rel += (cdf_rel - nuc.xs_cdf_[i_E_low + i_E_rel]) / m; // perform rejection sampling on cosine between // neutron and target velocities