From 3de5b4a9f84ff167f456871bbd5bc04c2a734d65 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 2 Sep 2026 00:53:56 +0300 Subject: [PATCH 1/4] fix --- src/distribution_energy.cpp | 9 ++++-- src/secondary_correlated.cpp | 54 +++++++++++++++++++++--------------- src/secondary_kalbach.cpp | 54 +++++++++++++++++++++--------------- 3 files changed, 71 insertions(+), 46 deletions(-) diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp index 7712c0d763f..68784d8555c 100644 --- a/src/distribution_energy.cpp +++ b/src/distribution_energy.cpp @@ -213,8 +213,13 @@ double ContinuousTabular::sample(double E, uint64_t* seed) const double E_l_k = distribution_[l].e_out[k]; if (k < n_discrete) { - // Discrete case - return E_l_k; + // Discrete case. Discrete lines correspond by index between adjacent + // incident energies, so interpolate the line energy rather than returning + // the value tabulated at grid l alone. Histogram interpolation on the + // incident energy grid pins l to i, so no interpolation is applied there. + const auto& e_lo = distribution_[i].e_out; + const auto& e_hi = distribution_[i + 1].e_out; + return e_lo[k] + (histogram_interp ? 0.0 : r) * (e_hi[k] - e_lo[k]); } else { // Continuous case double p_l_k = distribution_[l].p[k]; diff --git a/src/secondary_correlated.cpp b/src/secondary_correlated.cpp index 642ecd82b78..16424783cea 100644 --- a/src/secondary_correlated.cpp +++ b/src/secondary_correlated.cpp @@ -166,23 +166,9 @@ Distribution& CorrelatedAngleEnergy::sample_dist( // Sample between the ith and [i+1]th bin int l = r > prn(seed) ? i + 1 : i; - // Interpolation for energy E1 and EK - int n_energy_out = distribution_[i].e_out.size(); - int n_discrete = distribution_[i].n_discrete; - double E_i_1 = distribution_[i].e_out[n_discrete]; - double E_i_K = distribution_[i].e_out[n_energy_out - 1]; - - n_energy_out = distribution_[i + 1].e_out.size(); - n_discrete = distribution_[i + 1].n_discrete; - double E_i1_1 = distribution_[i + 1].e_out[n_discrete]; - double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1]; - - double E_1 = E_i_1 + r * (E_i1_1 - E_i_1); - double E_K = E_i_K + r * (E_i1_K - E_i_K); - // Determine outgoing energy bin - n_energy_out = distribution_[l].e_out.size(); - n_discrete = distribution_[l].n_discrete; + int n_energy_out = distribution_[l].e_out.size(); + int n_discrete = distribution_[l].n_discrete; double r1 = prn(seed); double c_k = distribution_[l].c[0]; int k = 0; @@ -211,8 +197,12 @@ Distribution& CorrelatedAngleEnergy::sample_dist( double E_l_k = distribution_[l].e_out[k]; if (k < n_discrete) { - // Discrete case - E_out = E_l_k; + // Discrete case. Discrete lines correspond by index between adjacent + // incident energies, so interpolate the line energy rather than using the + // value tabulated at grid l alone. + const auto& e_lo = distribution_[i].e_out; + const auto& e_hi = distribution_[i + 1].e_out; + E_out = e_lo[k] + r * (e_hi[k] - e_lo[k]); } else { // Continuous case double p_l_k = distribution_[l].p[k]; @@ -240,11 +230,31 @@ Distribution& CorrelatedAngleEnergy::sample_dist( } } + // Unit-base interpolation between the continua at i and i + 1. The ACE + // format requires the number of discrete lines to match at every incident + // energy, but does not require a continuum to be present at each one, so + // check for each side and fall back to whichever one has a continuum. + const int n_i = distribution_[i].e_out.size(); + const int nd_i = distribution_[i].n_discrete; + const int n_i1 = distribution_[i + 1].e_out.size(); + const int nd_i1 = distribution_[i + 1].n_discrete; + const bool cont_i = nd_i < n_i; + const bool cont_i1 = nd_i1 < n_i1; + + const double r_c = cont_i ? (cont_i1 ? r : 0.0) : 1.0; + const double E_i_1 = cont_i ? distribution_[i].e_out[nd_i] : 0.0; + const double E_i_K = cont_i ? distribution_[i].e_out[n_i - 1] : 0.0; + const double E_i1_1 = cont_i1 ? distribution_[i + 1].e_out[nd_i1] : 0.0; + const double E_i1_K = cont_i1 ? distribution_[i + 1].e_out[n_i1 - 1] : 0.0; + + const double E_1 = E_i_1 + r_c * (E_i1_1 - E_i_1); + const double E_K = E_i_K + r_c * (E_i1_K - E_i_K); + // Now interpolate between incident energy bins i and i + 1 - if (l == i) { - E_out = E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1); - } else { - E_out = E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1); + const double E_lo = (l == i) ? E_i_1 : E_i1_1; + const double E_hi = (l == i) ? E_i_K : E_i1_K; + if (E_hi != E_lo) { + E_out = E_1 + (E_out - E_lo) * (E_K - E_1) / (E_hi - E_lo); } } diff --git a/src/secondary_kalbach.cpp b/src/secondary_kalbach.cpp index a5d9b9ad0b4..680ca0bc72f 100644 --- a/src/secondary_kalbach.cpp +++ b/src/secondary_kalbach.cpp @@ -125,23 +125,9 @@ void KalbachMann::sample_params( // Sample between the ith and [i+1]th bin int l = r > prn(seed) ? i + 1 : i; - // Interpolation for energy E1 and EK - int n_energy_out = distribution_[i].e_out.size(); - int n_discrete = distribution_[i].n_discrete; - double E_i_1 = distribution_[i].e_out[n_discrete]; - double E_i_K = distribution_[i].e_out[n_energy_out - 1]; - - n_energy_out = distribution_[i + 1].e_out.size(); - n_discrete = distribution_[i + 1].n_discrete; - double E_i1_1 = distribution_[i + 1].e_out[n_discrete]; - double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1]; - - double E_1 = E_i_1 + r * (E_i1_1 - E_i_1); - double E_K = E_i_K + r * (E_i1_K - E_i_K); - // Determine outgoing energy bin - n_energy_out = distribution_[l].e_out.size(); - n_discrete = distribution_[l].n_discrete; + int n_energy_out = distribution_[l].e_out.size(); + int n_discrete = distribution_[l].n_discrete; double r1 = prn(seed); double c_k = distribution_[l].c[0]; int k = 0; @@ -170,8 +156,12 @@ void KalbachMann::sample_params( double E_l_k = distribution_[l].e_out[k]; if (k < n_discrete) { - // Discrete case - E_out = E_l_k; + // Discrete case. Discrete lines correspond by index between adjacent + // incident energies, so interpolate the line energy rather than using the + // value tabulated at grid l alone. + const auto& e_lo = distribution_[i].e_out; + const auto& e_hi = distribution_[i + 1].e_out; + E_out = e_lo[k] + r * (e_hi[k] - e_lo[k]); km_r = distribution_[l].r[k]; km_a = distribution_[l].a[k]; @@ -215,11 +205,31 @@ void KalbachMann::sample_params( (distribution_[l].a[k + 1] - distribution_[l].a[k]); } + // Unit-base interpolation between the continua at i and i + 1. The ACE + // format requires the number of discrete lines to match at every incident + // energy, but does not require a continuum to be present at each one, so + // check for each side and fall back to whichever one has a continuum. + const int n_i = distribution_[i].e_out.size(); + const int nd_i = distribution_[i].n_discrete; + const int n_i1 = distribution_[i + 1].e_out.size(); + const int nd_i1 = distribution_[i + 1].n_discrete; + const bool cont_i = nd_i < n_i; + const bool cont_i1 = nd_i1 < n_i1; + + const double r_c = cont_i ? (cont_i1 ? r : 0.0) : 1.0; + const double E_i_1 = cont_i ? distribution_[i].e_out[nd_i] : 0.0; + const double E_i_K = cont_i ? distribution_[i].e_out[n_i - 1] : 0.0; + const double E_i1_1 = cont_i1 ? distribution_[i + 1].e_out[nd_i1] : 0.0; + const double E_i1_K = cont_i1 ? distribution_[i + 1].e_out[n_i1 - 1] : 0.0; + + const double E_1 = E_i_1 + r_c * (E_i1_1 - E_i_1); + const double E_K = E_i_K + r_c * (E_i1_K - E_i_K); + // Now interpolate between incident energy bins i and i + 1 - if (l == i) { - E_out = E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1); - } else { - E_out = E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1); + const double E_lo = (l == i) ? E_i_1 : E_i1_1; + const double E_hi = (l == i) ? E_i_K : E_i1_K; + if (E_hi != E_lo) { + E_out = E_1 + (E_out - E_lo) * (E_K - E_1) / (E_hi - E_lo); } } } From 12c3f499c0bf28fee9c79d2c38f328be6af9580c Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 2 Sep 2026 01:25:45 +0300 Subject: [PATCH 2/4] fix --- src/distribution_energy.cpp | 52 +++++++++++++++++++++---------------- 1 file changed, 29 insertions(+), 23 deletions(-) diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp index 68784d8555c..d1564432ad8 100644 --- a/src/distribution_energy.cpp +++ b/src/distribution_energy.cpp @@ -256,30 +256,36 @@ double ContinuousTabular::sample(double E, uint64_t* seed) const "distribution."}; } - // Now interpolate between incident energy bins i and i + 1 - if (!histogram_interp && n_energy_out > 1) { - // Interpolation for energy E1 and EK - n_energy_out = distribution_[i].e_out.size(); - n_discrete = distribution_[i].n_discrete; - const double E_i_1 = distribution_[i].e_out[n_discrete]; - const double E_i_K = distribution_[i].e_out[n_energy_out - 1]; - - n_energy_out = distribution_[i + 1].e_out.size(); - n_discrete = distribution_[i + 1].n_discrete; - const double E_i1_1 = distribution_[i + 1].e_out[n_discrete]; - const double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1]; - - const double E_1 = E_i_1 + r * (E_i1_1 - E_i_1); - const double E_K = E_i_K + r * (E_i1_K - E_i_K); - - if (l == i) { - return E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1); - } else { - return E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1); - } - } else { + // Histogram interpolation on the incident energy grid means the table at i + // applies as tabulated, with no interpolation between grids + if (histogram_interp) return E_out; - } + + // Unit-base interpolation between the continua at i and i + 1. The ACE + // format requires the number of discrete lines to match at every incident + // energy, but does not require a continuum to be present at each one, so + // check for each side and fall back to whichever one has a continuum. + const int n_i = distribution_[i].e_out.size(); + const int nd_i = distribution_[i].n_discrete; + const int n_i1 = distribution_[i + 1].e_out.size(); + const int nd_i1 = distribution_[i + 1].n_discrete; + const bool cont_i = nd_i < n_i; + const bool cont_i1 = nd_i1 < n_i1; + + const double r_c = cont_i ? (cont_i1 ? r : 0.0) : 1.0; + const double E_i_1 = cont_i ? distribution_[i].e_out[nd_i] : 0.0; + const double E_i_K = cont_i ? distribution_[i].e_out[n_i - 1] : 0.0; + const double E_i1_1 = cont_i1 ? distribution_[i + 1].e_out[nd_i1] : 0.0; + const double E_i1_K = cont_i1 ? distribution_[i + 1].e_out[n_i1 - 1] : 0.0; + + const double E_1 = E_i_1 + r_c * (E_i1_1 - E_i_1); + const double E_K = E_i_K + r_c * (E_i1_K - E_i_K); + + // Now interpolate between incident energy bins i and i + 1 + const double E_lo = (l == i) ? E_i_1 : E_i1_1; + const double E_hi = (l == i) ? E_i_K : E_i1_K; + return (E_hi != E_lo) ? E_1 + (E_out - E_lo) * (E_K - E_1) / (E_hi - E_lo) + : E_out; } } From deff3ef04217d5d14eb17ee7f06cdc89b7687149 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 2 Sep 2026 01:59:28 +0300 Subject: [PATCH 3/4] update regression test --- tests/regression_tests/weightwindows/shared/results_true.dat | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/regression_tests/weightwindows/shared/results_true.dat b/tests/regression_tests/weightwindows/shared/results_true.dat index 6d1bb9789c5..367ce39b860 100644 --- a/tests/regression_tests/weightwindows/shared/results_true.dat +++ b/tests/regression_tests/weightwindows/shared/results_true.dat @@ -1 +1 @@ -19b676a5e0f2784747d95fe884a44af3b31b4f5193e5ac514378d6d83cfa8ae9f4722da408a98b9b9cd3e369517d0b29bc523fd6f5492856f457e14d15007fa0 \ No newline at end of file +3916cc59d580f8caf1dac2ae30d7a5926f650ab1e314af4ecedd3521c87aa7ab60f3ddfee940ed221dbfd2a40deb0f726eb6252ca3333445c7aa8155e7003094 \ No newline at end of file From 771b54cd4a4cd61bb306dbd9a871481544048c75 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 12 Sep 2026 11:50:22 +0000 Subject: [PATCH 4/4] Guard secondary distribution sampling against a single incident energy get_energy_index() returns i = 0 with an interpolation factor of zero when a distribution is tabulated at a single incident energy, since no interval exists. The samplers nevertheless read distribution_[i + 1] unconditionally, which is an out-of-bounds access on such a table -- the crash reported in openmc-dev/openmc#1098 for a delayed neutron spectrum with one incident energy. Add upper_energy_index(), which clamps the upper index of the incident energy interval to the end of the grid, and use it in ContinuousTabular, CorrelatedAngleEnergy, KalbachMann and AngleDistribution. For a grid with two or more points get_energy_index() keeps i <= n - 2, so the helper always returns i + 1 and results are unchanged. For a single incident energy the interpolation factor is zero and interpolating between the clamped indices is a no-op, so the distribution is sampled as tabulated instead of reading past the end of the array. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01MxFQpQ6zFTaLKsTKtVxy9V --- include/openmc/math_functions.h | 9 +++++++++ src/distribution_angle.cpp | 5 +++-- src/distribution_energy.cpp | 14 ++++++++------ src/secondary_correlated.cpp | 14 ++++++++------ src/secondary_kalbach.cpp | 14 ++++++++------ 5 files changed, 36 insertions(+), 20 deletions(-) diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h index d3bccca7b5d..3feb404f0cc 100644 --- a/include/openmc/math_functions.h +++ b/include/openmc/math_functions.h @@ -237,6 +237,15 @@ double cyl_bessel_j(int n, double x); void get_energy_index( const vector& energies, double E, int& i, double& f); +//! Upper index of the interval located by get_energy_index. Normally i + 1, +//! but a distribution tabulated at a single incident energy has no upper +//! point; the interpolation factor is zero there, so returning i makes +//! interpolating between the two indices a no-op. +inline int upper_energy_index(const vector& energies, int i) +{ + return (i + 1 < static_cast(energies.size())) ? i + 1 : i; +} + //============================================================================== //! Calculate the cumulative distribution function of the standard normal //! distribution at a given value. diff --git a/src/distribution_angle.cpp b/src/distribution_angle.cpp index f3efa27f507..10ab62c9836 100644 --- a/src/distribution_angle.cpp +++ b/src/distribution_angle.cpp @@ -71,7 +71,7 @@ double AngleDistribution::sample(double E, uint64_t* seed) const // Sample between the ith and (i+1)th bin if (r > prn(seed)) - ++i; + i = upper_energy_index(energy_, i); // Sample i-th distribution double mu = distribution_[i]->sample(seed).first; @@ -88,7 +88,8 @@ double AngleDistribution::evaluate(double E, double mu) const int i; double r; get_energy_index(energy_, E, i, r); - return r * distribution_[i + 1]->evaluate(mu) + + const int i1 = upper_energy_index(energy_, i); + return r * distribution_[i1]->evaluate(mu) + (1.0 - r) * distribution_[i]->evaluate(mu); } diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp index e3e840dbb16..55d4894d613 100644 --- a/src/distribution_energy.cpp +++ b/src/distribution_energy.cpp @@ -162,12 +162,14 @@ double ContinuousTabular::sample(double E, uint64_t* seed) const double r; get_energy_index(energy_, E, i, r); + const int i1 = upper_energy_index(energy_, i); + // Sample between the ith and [i+1]th bin int l; if (histogram_interp) { l = i; } else { - l = r > prn(seed) ? i + 1 : i; + l = r > prn(seed) ? i1 : i; } // Determine outgoing energy bin @@ -207,7 +209,7 @@ double ContinuousTabular::sample(double E, uint64_t* seed) const // the value tabulated at grid l alone. Histogram interpolation on the // incident energy grid pins l to i, so no interpolation is applied there. const auto& e_lo = distribution_[i].e_out; - const auto& e_hi = distribution_[i + 1].e_out; + const auto& e_hi = distribution_[i1].e_out; return e_lo[k] + (histogram_interp ? 0.0 : r) * (e_hi[k] - e_lo[k]); } else { // Continuous case @@ -256,16 +258,16 @@ double ContinuousTabular::sample(double E, uint64_t* seed) const // check for each side and fall back to whichever one has a continuum. const int n_i = distribution_[i].e_out.size(); const int nd_i = distribution_[i].n_discrete; - const int n_i1 = distribution_[i + 1].e_out.size(); - const int nd_i1 = distribution_[i + 1].n_discrete; + const int n_i1 = distribution_[i1].e_out.size(); + const int nd_i1 = distribution_[i1].n_discrete; const bool cont_i = nd_i < n_i; const bool cont_i1 = nd_i1 < n_i1; const double r_c = cont_i ? (cont_i1 ? r : 0.0) : 1.0; const double E_i_1 = cont_i ? distribution_[i].e_out[nd_i] : 0.0; const double E_i_K = cont_i ? distribution_[i].e_out[n_i - 1] : 0.0; - const double E_i1_1 = cont_i1 ? distribution_[i + 1].e_out[nd_i1] : 0.0; - const double E_i1_K = cont_i1 ? distribution_[i + 1].e_out[n_i1 - 1] : 0.0; + const double E_i1_1 = cont_i1 ? distribution_[i1].e_out[nd_i1] : 0.0; + const double E_i1_K = cont_i1 ? distribution_[i1].e_out[n_i1 - 1] : 0.0; const double E_1 = E_i_1 + r_c * (E_i1_1 - E_i_1); const double E_K = E_i_K + r_c * (E_i1_K - E_i_K); diff --git a/src/secondary_correlated.cpp b/src/secondary_correlated.cpp index 16424783cea..5bb7d56b0a5 100644 --- a/src/secondary_correlated.cpp +++ b/src/secondary_correlated.cpp @@ -163,8 +163,10 @@ Distribution& CorrelatedAngleEnergy::sample_dist( double r; get_energy_index(energy_, E_in, i, r); + const int i1 = upper_energy_index(energy_, i); + // Sample between the ith and [i+1]th bin - int l = r > prn(seed) ? i + 1 : i; + int l = r > prn(seed) ? i1 : i; // Determine outgoing energy bin int n_energy_out = distribution_[l].e_out.size(); @@ -201,7 +203,7 @@ Distribution& CorrelatedAngleEnergy::sample_dist( // incident energies, so interpolate the line energy rather than using the // value tabulated at grid l alone. const auto& e_lo = distribution_[i].e_out; - const auto& e_hi = distribution_[i + 1].e_out; + const auto& e_hi = distribution_[i1].e_out; E_out = e_lo[k] + r * (e_hi[k] - e_lo[k]); } else { // Continuous case @@ -236,16 +238,16 @@ Distribution& CorrelatedAngleEnergy::sample_dist( // check for each side and fall back to whichever one has a continuum. const int n_i = distribution_[i].e_out.size(); const int nd_i = distribution_[i].n_discrete; - const int n_i1 = distribution_[i + 1].e_out.size(); - const int nd_i1 = distribution_[i + 1].n_discrete; + const int n_i1 = distribution_[i1].e_out.size(); + const int nd_i1 = distribution_[i1].n_discrete; const bool cont_i = nd_i < n_i; const bool cont_i1 = nd_i1 < n_i1; const double r_c = cont_i ? (cont_i1 ? r : 0.0) : 1.0; const double E_i_1 = cont_i ? distribution_[i].e_out[nd_i] : 0.0; const double E_i_K = cont_i ? distribution_[i].e_out[n_i - 1] : 0.0; - const double E_i1_1 = cont_i1 ? distribution_[i + 1].e_out[nd_i1] : 0.0; - const double E_i1_K = cont_i1 ? distribution_[i + 1].e_out[n_i1 - 1] : 0.0; + const double E_i1_1 = cont_i1 ? distribution_[i1].e_out[nd_i1] : 0.0; + const double E_i1_K = cont_i1 ? distribution_[i1].e_out[n_i1 - 1] : 0.0; const double E_1 = E_i_1 + r_c * (E_i1_1 - E_i_1); const double E_K = E_i_K + r_c * (E_i1_K - E_i_K); diff --git a/src/secondary_kalbach.cpp b/src/secondary_kalbach.cpp index 680ca0bc72f..34cd76106da 100644 --- a/src/secondary_kalbach.cpp +++ b/src/secondary_kalbach.cpp @@ -122,8 +122,10 @@ void KalbachMann::sample_params( double r; get_energy_index(energy_, E_in, i, r); + const int i1 = upper_energy_index(energy_, i); + // Sample between the ith and [i+1]th bin - int l = r > prn(seed) ? i + 1 : i; + int l = r > prn(seed) ? i1 : i; // Determine outgoing energy bin int n_energy_out = distribution_[l].e_out.size(); @@ -160,7 +162,7 @@ void KalbachMann::sample_params( // incident energies, so interpolate the line energy rather than using the // value tabulated at grid l alone. const auto& e_lo = distribution_[i].e_out; - const auto& e_hi = distribution_[i + 1].e_out; + const auto& e_hi = distribution_[i1].e_out; E_out = e_lo[k] + r * (e_hi[k] - e_lo[k]); km_r = distribution_[l].r[k]; km_a = distribution_[l].a[k]; @@ -211,16 +213,16 @@ void KalbachMann::sample_params( // check for each side and fall back to whichever one has a continuum. const int n_i = distribution_[i].e_out.size(); const int nd_i = distribution_[i].n_discrete; - const int n_i1 = distribution_[i + 1].e_out.size(); - const int nd_i1 = distribution_[i + 1].n_discrete; + const int n_i1 = distribution_[i1].e_out.size(); + const int nd_i1 = distribution_[i1].n_discrete; const bool cont_i = nd_i < n_i; const bool cont_i1 = nd_i1 < n_i1; const double r_c = cont_i ? (cont_i1 ? r : 0.0) : 1.0; const double E_i_1 = cont_i ? distribution_[i].e_out[nd_i] : 0.0; const double E_i_K = cont_i ? distribution_[i].e_out[n_i - 1] : 0.0; - const double E_i1_1 = cont_i1 ? distribution_[i + 1].e_out[nd_i1] : 0.0; - const double E_i1_K = cont_i1 ? distribution_[i + 1].e_out[n_i1 - 1] : 0.0; + const double E_i1_1 = cont_i1 ? distribution_[i1].e_out[nd_i1] : 0.0; + const double E_i1_K = cont_i1 ? distribution_[i1].e_out[n_i1 - 1] : 0.0; const double E_1 = E_i_1 + r_c * (E_i1_1 - E_i_1); const double E_K = E_i_K + r_c * (E_i1_K - E_i_K);