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

Off by 1 error on Nuclide::cs_cdf

Browse source 
Error occured in regression_tests/resonance_scattering/test.py
This commit is contained in:
Olaf Schumann 2022-01-06 14:51:38 +01:00
parent b150a58a11
commit 9c7dcbba63
2 changed files with 15 additions and 14 deletions
Download patch file Download diff file Expand all files Collapse all files

27
src/physics.cpp
View file

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