diff --git a/src/secondary_correlated.cpp b/src/secondary_correlated.cpp index 8e0f23d520..4a1ccad6e4 100644 --- a/src/secondary_correlated.cpp +++ b/src/secondary_correlated.cpp @@ -2,14 +2,14 @@ #include // for copy #include -#include // for size_t -#include // for back_inserter +#include // for size_t +#include // for back_inserter #include "xtensor/xarray.hpp" #include "xtensor/xview.hpp" -#include "openmc/hdf5_interface.h" #include "openmc/endf.h" +#include "openmc/hdf5_interface.h" #include "openmc/random_lcg.h" #include "openmc/search.h" @@ -63,7 +63,7 @@ CorrelatedAngleEnergy::CorrelatedAngleEnergy(hid_t group) int j = offsets[i]; int n; if (i < n_energy - 1) { - n = offsets[i+1] - j; + n = offsets[i + 1] - j; } else { n = eout.shape()[1] - j; } @@ -74,9 +74,9 @@ CorrelatedAngleEnergy::CorrelatedAngleEnergy(hid_t group) d.n_discrete = n_discrete[i]; // Copy data - d.e_out = xt::view(eout, 0, xt::range(j, j+n)); - d.p = xt::view(eout, 1, xt::range(j, j+n)); - d.c = xt::view(eout, 2, xt::range(j, j+n)); + d.e_out = xt::view(eout, 0, xt::range(j, j + n)); + d.p = xt::view(eout, 1, xt::range(j, j + n)); + d.c = xt::view(eout, 2, xt::range(j, j + n)); // To get answers that match ACE data, for now we still use the tabulated // CDF values that were passed through to the HDF5 library. At a later @@ -88,20 +88,20 @@ CorrelatedAngleEnergy::CorrelatedAngleEnergy(hid_t group) if (k == 0) { d.c[k] = d.p[k]; } else { - d.c[k] = d.c[k-1] + d.p[k]; + d.c[k] = d.c[k - 1] + d.p[k]; } } // Continuous portion for (int k = d.n_discrete; k < n; ++k) { if (k == d.n_discrete) { - d.c[k] = d.c[k-1] + d.p[k]; + d.c[k] = d.c[k - 1] + d.p[k]; } else { if (d.interpolation == Interpolation::histogram) { - d.c[k] = d.c[k-1] + d.p[k-1]*(d.e_out[k] - d.e_out[k-1]); + d.c[k] = d.c[k - 1] + d.p[k - 1] * (d.e_out[k] - d.e_out[k - 1]); } else if (d.interpolation == Interpolation::lin_lin) { - d.c[k] = d.c[k-1] + 0.5*(d.p[k-1] + d.p[k]) * - (d.e_out[k] - d.e_out[k-1]); + d.c[k] = d.c[k - 1] + 0.5 * (d.p[k - 1] + d.p[k]) * + (d.e_out[k] - d.e_out[k - 1]); } } } @@ -119,7 +119,7 @@ CorrelatedAngleEnergy::CorrelatedAngleEnergy(hid_t group) int offset_mu = std::lround(eout(4, offsets[i] + j)); int m; if (offsets[i] + j + 1 < eout.shape()[1]) { - m = std::lround(eout(4, offsets[i]+j+1)) - offset_mu; + m = std::lround(eout(4, offsets[i] + j + 1)) - offset_mu; } else { m = mu.shape()[1] - offset_mu; } @@ -128,7 +128,8 @@ CorrelatedAngleEnergy::CorrelatedAngleEnergy(hid_t group) // may be given as discrete mu values. In this case, interpolation values // of zero appear in the HDF5 file. Here we change it to a 1 so that // int2interp doesn't fail. - if (interp_mu == 0) interp_mu = 1; + if (interp_mu == 0) + interp_mu = 1; auto interp = int2interp(interp_mu); auto xs = xt::view(mu, 0, xt::range(offset_mu, offset_mu + m)); @@ -143,7 +144,7 @@ CorrelatedAngleEnergy::CorrelatedAngleEnergy(hid_t group) // CDF values that were passed through to the HDF5 library. At a later // time, we can remove the CDF values from the HDF5 library and // reconstruct them using the PDF - Tabular* mudist = new Tabular{x.data(), p.data(), m, interp, c.data()}; + Tabular* mudist = new Tabular {x.data(), p.data(), m, interp, c.data()}; d.angle.emplace_back(mudist); } // outgoing energies @@ -152,15 +153,15 @@ CorrelatedAngleEnergy::CorrelatedAngleEnergy(hid_t group) } // incoming energies } -void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu, - uint64_t* seed) const +void CorrelatedAngleEnergy::sample( + double E_in, double& E_out, double& mu, uint64_t* seed) const { // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<< // Before the secondary distribution refactor, an isotropic polar cosine was // always sampled but then overwritten with the polar cosine sampled from the // correlated distribution. To preserve the random number stream, we keep // this dummy sampling here but can remove it later (will change answers) - mu = 2.0*prn(seed) - 1.0; + mu = 2.0 * prn(seed) - 1.0; // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<< // Find energy bin and calculate interpolation factor -- if the energy is @@ -176,7 +177,7 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu, r = 1.0; } else { i = lower_bound_index(energy_.begin(), energy_.end(), E_in); - r = (E_in - energy_[i]) / (energy_[i+1] - energy_[i]); + r = (E_in - energy_[i]) / (energy_[i + 1] - energy_[i]); } // Sample between the ith and [i+1]th bin @@ -188,13 +189,13 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu, 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]; + 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); + 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(); @@ -218,8 +219,9 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu, double c_k1; for (int j = n_discrete; j < end; ++j) { k = j; - c_k1 = distribution_[l].c[k+1]; - if (r1 < c_k1) break; + c_k1 = distribution_[l].c[k + 1]; + if (r1 < c_k1) + break; k = j + 1; c_k = c_k1; } @@ -229,37 +231,40 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu, if (distribution_[l].interpolation == Interpolation::histogram) { // Histogram interpolation if (p_l_k > 0.0 && k >= n_discrete) { - E_out = E_l_k + (r1 - c_k)/p_l_k; + E_out = E_l_k + (r1 - c_k) / p_l_k; } else { E_out = E_l_k; } } else if (distribution_[l].interpolation == Interpolation::lin_lin) { // Linear-linear interpolation - double E_l_k1 = distribution_[l].e_out[k+1]; - double p_l_k1 = distribution_[l].p[k+1]; + double E_l_k1 = distribution_[l].e_out[k + 1]; + double p_l_k1 = distribution_[l].p[k + 1]; - double frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k); + double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k); if (frac == 0.0) { - E_out = E_l_k + (r1 - c_k)/p_l_k; + E_out = E_l_k + (r1 - c_k) / p_l_k; } else { - E_out = E_l_k + (std::sqrt(std::max(0.0, p_l_k*p_l_k + - 2.0*frac*(r1 - c_k))) - p_l_k)/frac; + E_out = + E_l_k + + (std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) - + p_l_k) / + frac; } - } // Now interpolate between incident energy bins i and i + 1 - if (k >= n_discrete){ + if (k >= n_discrete) { if (l == i) { - E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1); + 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); + E_out = E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1); } } // Find correlated angular distribution for closest outgoing energy bin - if (r1 - c_k < c_k1 - r1) { + if (r1 - c_k < c_k1 - r1 || + distribution_[l].interpolation == Interpolation::histogram) { mu = distribution_[l].angle[k]->sample(seed); } else { mu = distribution_[l].angle[k + 1]->sample(seed); diff --git a/tests/regression_tests/deplete/last_step_reference_materials.xml b/tests/regression_tests/deplete/last_step_reference_materials.xml index 10815c9a2f..70554af9ad 100644 --- a/tests/regression_tests/deplete/last_step_reference_materials.xml +++ b/tests/regression_tests/deplete/last_step_reference_materials.xml @@ -4,1001 +4,1001 @@ - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + diff --git a/tests/regression_tests/deplete/test_reference.h5 b/tests/regression_tests/deplete/test_reference.h5 index 7331f40596..32e13c74f3 100644 Binary files a/tests/regression_tests/deplete/test_reference.h5 and b/tests/regression_tests/deplete/test_reference.h5 differ diff --git a/tests/unit_tests/test_deplete_resultslist.py b/tests/unit_tests/test_deplete_resultslist.py index be81da1377..94d533d9d3 100644 --- a/tests/unit_tests/test_deplete_resultslist.py +++ b/tests/unit_tests/test_deplete_resultslist.py @@ -22,7 +22,7 @@ def test_get_atoms(res): t_ref = np.array([0.0, 1296000.0, 2592000.0, 3888000.0]) n_ref = np.array( - [6.67473282e+08, 3.76987065e+14, 3.68586723e+14, 3.91338392e+14]) + [6.67473282e+08, 3.76987065e+14, 3.31744559e+14, 3.78705433e+14]) np.testing.assert_allclose(t, t_ref) np.testing.assert_allclose(n, n_ref) @@ -48,8 +48,8 @@ def test_get_reaction_rate(res): t, r = res.get_reaction_rate("1", "Xe135", "(n,gamma)") t_ref = [0.0, 1296000.0, 2592000.0, 3888000.0] - n_ref = [6.67473282e+08, 3.76987065e+14, 3.68586723e+14, 3.91338392e+14] - xs_ref = [3.32282064e-05, 2.76208092e-05, 4.10987995e-05, 3.72454755e-05] + n_ref = [6.67473282e+08, 3.76987065e+14, 3.31744559e+14, 3.78705433e+14] + xs_ref = [3.32282064e-05, 2.76403474e-05, 3.57047568e-05, 3.03844011e-05] np.testing.assert_allclose(t, t_ref) np.testing.assert_allclose(r, np.array(n_ref) * xs_ref) @@ -60,8 +60,8 @@ def test_get_eigenvalue(res): t, k = res.get_eigenvalue() t_ref = [0.0, 1296000.0, 2592000.0, 3888000.0] - k_ref = [1.16984322, 1.19097429, 1.03012517, 1.20045563] - u_ref = [0.0375587, 0.03476389, 0.07215969, 0.02839639] + k_ref = [1.16984322, 1.24064521, 1.1826205 , 1.21236309] + u_ref = [0.0375587018699, 0.0255584778413, 0.0241749230118, 0.0217929393239] np.testing.assert_allclose(t, t_ref) np.testing.assert_allclose(k[:, 0], k_ref)