crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge pull request #1821 from HunterBelanger/correlated-angle-energy-fix

Browse source 
Fixes angle sampling in CorrelatedAngleEnergy
This commit is contained in:
Paul Romano 2021-04-26 14:44:22 -05:00 • committed by GitHub
commit 7f7977091d
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
4 changed files with 697 additions and 692 deletions
Download patch file Download diff file Expand all files Collapse all files

83
src/secondary_correlated.cpp
View file

@ -2,14 +2,14 @@
#include <algorithm> // for copy
#include <cmath>
#include <cstddef> // for size_t
#include <iterator> // for back_inserter
#include <cstddef> // for size_t
#include <iterator> // 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);

1296
tests/regression_tests/deplete/last_step_reference_materials.xml
View file

File diff suppressed because it is too large Load diff

BIN
tests/regression_tests/deplete/test_reference.h5
View file

Binary file not shown.

10
tests/unit_tests/test_deplete_resultslist.py
View file

@ -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)