Xtensor for photoelectric cross sections matches the flux spectrum for uranium of the previous iteration of OpenMC

This commit is contained in:
myerspat 2022-06-09 09:52:50 -05:00
parent 6df0286b2e
commit 61db44f40e
2 changed files with 20 additions and 17 deletions

View file

@ -29,6 +29,7 @@
#include <algorithm> // for max, min, max_element
#include <cmath> // for sqrt, exp, log, abs, copysign
#include <xtensor/xview.hpp>
namespace openmc {
@ -344,25 +345,24 @@ void sample_photon_reaction(Particle& p)
// Photoelectric effect
double prob_after = prob + micro.photoelectric;
int i_grid = micro.index_grid;
double f = micro.interp_factor;
const auto& xs_lower = xt::row(element.cross_sections_, i_grid);
const auto& xs_upper = xt::row(element.cross_sections_, i_grid + 1);
if (prob_after > cutoff) {
for (int i_shell = 0; i_shell < element.shells_.size(); ++i_shell) {
const auto& shell {element.shells_[i_shell]};
// Get grid index and interpolation factor
int i_grid = micro.index_grid;
double f = micro.interp_factor;
// Check threshold of reaction
int i_start = shell.threshold;
if (i_grid < i_start)
if (xs_lower(i_shell) == 0)
continue;
// Evaluation subshell photoionization cross section
double xs = std::exp(element.cross_sections_(i_grid - i_start) +
f * (element.cross_sections_(i_grid + 1 - i_start) -
element.cross_sections_(i_grid - i_start)));
// Evaluation subshell photoionization cross section
prob += std::exp(xs_lower(i_shell) +
f * (xs_upper(i_shell) -
xs_lower(i_shell)));
prob += xs;
if (prob > cutoff) {
double E_electron = p.E() - shell.binding_energy;