Hey everyone,
So, the reason we use “nearest” mode temperature interpolation is so OpenMC works as you’d expect if you have only one temperature available for a nuclear data library. The reason being that obviously you can’t interpolate if you only have one temperature.
I think we could just write the interpolation code to work like nearest mode if only one temperature is available. It would mean changing the interpolation code from looking like this (nuclide.cpp):
for (i_temp = 0; i_temp < kTs_.size() - 1; ++i_temp) {
if (kTs_[i_temp] <= kT && kT < kTs_[i_temp + 1]) break;
}
// Randomly sample between temperature i and i+1
f = (kT - kTs_[i_temp]) / (kTs_[i_temp + 1] - kTs_[i_temp]);
if (f > prn(p.current_seed())) ++i_temp;
break;
}
to
for (i_temp = 0; i_temp < kTs_.size() - 1; ++i_temp) {
if (kTs_[i_temp] <= kT && kT < kTs_[i_temp + 1]) break;
}
i_temp_upper = i_temp+static_cast<int>(kTs_.size() > 1);
// Randomly sample between temperature i and i+1
if (kT - kTs_[i_temp] > (kTs_[i_temp_upper] - kTs_[i_temp]) * prn(p.current_seed())) i_temp = i_temp_upper;
break;
}
Where you can see in the latter case, out-of-bounds accesses and division by zero are avoided if only one temperature is available. Now, the bigger problem is that in “nearest” mode, a random number is not consumed by the LCG, so switching to the above in regression tests that previously used “nearest” mode would cause a loss of reproducibility compared to the old test.
Interested to get thoughts on this.