From 687d43b0ab6e3e3ddfe867767c210f0a1cfb38a7 Mon Sep 17 00:00:00 2001 From: Patrick Myers Date: Wed, 25 May 2022 09:06:53 -0500 Subject: [PATCH 1/8] added external xtensor within PhotonInteraction for all shell photoelectric xs to utilize contiguous memory in calculate_xs --- include/openmc/photon.h | 2 +- src/photon.cpp | 26 +++++++++++--------------- src/physics.cpp | 6 +++--- 3 files changed, 15 insertions(+), 19 deletions(-) diff --git a/include/openmc/photon.h b/include/openmc/photon.h index e8ec82d107..50ee228c44 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -36,7 +36,6 @@ public: int threshold; double n_electrons; double binding_energy; - xt::xtensor cross_section; vector transitions; }; @@ -82,6 +81,7 @@ public: // Photoionization and atomic relaxation data vector shells_; + xt::xtensor cross_sections_; // Compton profile data xt::xtensor profile_pdf_; diff --git a/src/photon.cpp b/src/photon.cpp index b0cd8aa87d..3f9ff6dbf3 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -16,6 +16,9 @@ #include "xtensor/xbuilder.hpp" #include "xtensor/xoperation.hpp" #include "xtensor/xview.hpp" +#include "xtensor/xmath.hpp" +#include "xtensor/xslice.hpp" +#include "xtensor/xtensor_forward.hpp" #include #include @@ -125,6 +128,7 @@ PhotonInteraction::PhotonInteraction(hid_t group) } shells_.resize(n_shell); + cross_sections_.resize({energy_.size(), n_shell}); // Create mapping from designator to index std::unordered_map shell_map; @@ -155,13 +159,14 @@ PhotonInteraction::PhotonInteraction(hid_t group) read_attribute(tgroup, "num_electrons", shell.n_electrons); // Read subshell cross section + xt::xtensor xs; dset = open_dataset(tgroup, "xs"); read_attribute(dset, "threshold_idx", shell.threshold); close_dataset(dset); - read_dataset(tgroup, "xs", shell.cross_section); + read_dataset(tgroup, "xs", xs); - auto& xs = shell.cross_section; - xs = xt::where(xs > 0.0, xt::log(xs), -500.0); + auto cross_section = xt::view(cross_sections_, xt::range(shell.threshold, shell.threshold + xs.size()), i); + cross_section = xt::where(xs > 0.0, xt::log(xs), -500.0); if (object_exists(tgroup, "transitions")) { // Determine dimensions of transitions @@ -565,19 +570,10 @@ void PhotonInteraction::calculate_xs(Particle& p) const incoherent_(i_grid) + f * (incoherent_(i_grid + 1) - incoherent_(i_grid))); // Calculate microscopic photoelectric cross section - xs.photoelectric = 0.0; - for (const auto& shell : shells_) { - // Check threshold of reaction - int i_start = shell.threshold; - if (i_grid < i_start) - continue; + const auto& xs_upper = xt::row(cross_sections_, i_grid); + const auto& xs_lower = xt::row(cross_sections_, i_grid + 1); - // Evaluation subshell photoionization cross section - xs.photoelectric += - std::exp(shell.cross_section(i_grid - i_start) + - f * (shell.cross_section(i_grid + 1 - i_start) - - shell.cross_section(i_grid - i_start))); - } + xs.photoelectric = xt::sum(xt::exp(xs_upper + f * (xs_upper - xs_lower)))[0]; // Calculate microscopic pair production cross section xs.pair_production = std::exp( diff --git a/src/physics.cpp b/src/physics.cpp index 13dd063de5..38d74dbede 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -358,9 +358,9 @@ void sample_photon_reaction(Particle& p) continue; // Evaluation subshell photoionization cross section - double xs = std::exp(shell.cross_section(i_grid - i_start) + - f * (shell.cross_section(i_grid + 1 - i_start) - - shell.cross_section(i_grid - i_start))); + 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))); prob += xs; if (prob > cutoff) { From 6df0286b2e20e39b4eff8d93b262821086a9304f Mon Sep 17 00:00:00 2001 From: Patrick Myers Date: Wed, 25 May 2022 09:14:10 -0500 Subject: [PATCH 2/8] added clang-format --- src/photon.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/photon.cpp b/src/photon.cpp index 3f9ff6dbf3..1be515a621 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -14,11 +14,11 @@ #include "openmc/settings.h" #include "xtensor/xbuilder.hpp" -#include "xtensor/xoperation.hpp" -#include "xtensor/xview.hpp" #include "xtensor/xmath.hpp" +#include "xtensor/xoperation.hpp" #include "xtensor/xslice.hpp" #include "xtensor/xtensor_forward.hpp" +#include "xtensor/xview.hpp" #include #include @@ -165,7 +165,8 @@ PhotonInteraction::PhotonInteraction(hid_t group) close_dataset(dset); read_dataset(tgroup, "xs", xs); - auto cross_section = xt::view(cross_sections_, xt::range(shell.threshold, shell.threshold + xs.size()), i); + auto cross_section = xt::view(cross_sections_, + xt::range(shell.threshold, shell.threshold + xs.size()), i); cross_section = xt::where(xs > 0.0, xt::log(xs), -500.0); if (object_exists(tgroup, "transitions")) { From 61db44f40e33d132d895974852c448beeacc6770 Mon Sep 17 00:00:00 2001 From: myerspat Date: Thu, 9 Jun 2022 09:52:50 -0500 Subject: [PATCH 3/8] Xtensor for photoelectric cross sections matches the flux spectrum for uranium of the previous iteration of OpenMC --- src/photon.cpp | 15 +++++++++------ src/physics.cpp | 22 +++++++++++----------- 2 files changed, 20 insertions(+), 17 deletions(-) diff --git a/src/photon.cpp b/src/photon.cpp index 1be515a621..bb715cdeaf 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -128,7 +128,7 @@ PhotonInteraction::PhotonInteraction(hid_t group) } shells_.resize(n_shell); - cross_sections_.resize({energy_.size(), n_shell}); + cross_sections_ = xt::zeros({energy_.size(), n_shell}); // Create mapping from designator to index std::unordered_map shell_map; @@ -166,8 +166,8 @@ PhotonInteraction::PhotonInteraction(hid_t group) read_dataset(tgroup, "xs", xs); auto cross_section = xt::view(cross_sections_, - xt::range(shell.threshold, shell.threshold + xs.size()), i); - cross_section = xt::where(xs > 0.0, xt::log(xs), -500.0); + xt::range(shell.threshold, xt::placeholders::_), i); + cross_section = xt::where(xs > 0, xt::log(xs), 0); if (object_exists(tgroup, "transitions")) { // Determine dimensions of transitions @@ -571,10 +571,13 @@ void PhotonInteraction::calculate_xs(Particle& p) const incoherent_(i_grid) + f * (incoherent_(i_grid + 1) - incoherent_(i_grid))); // Calculate microscopic photoelectric cross section - const auto& xs_upper = xt::row(cross_sections_, i_grid); - const auto& xs_lower = xt::row(cross_sections_, i_grid + 1); + xs.photoelectric = 0.0; + const auto& xs_lower = xt::row(cross_sections_, i_grid); + const auto& xs_upper = xt::row(cross_sections_, i_grid + 1); - xs.photoelectric = xt::sum(xt::exp(xs_upper + f * (xs_upper - xs_lower)))[0]; + for (int i = 0; i < xs_upper.size(); ++i) + if (xs_lower(i) != 0) + xs.photoelectric += std::exp(xs_lower(i) + f * (xs_upper(i) - xs_lower(i))); // Calculate microscopic pair production cross section xs.pair_production = std::exp( diff --git a/src/physics.cpp b/src/physics.cpp index 38d74dbede..26b3a3520a 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -29,6 +29,7 @@ #include // for max, min, max_element #include // for sqrt, exp, log, abs, copysign +#include 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; From 70b1706460166ff6d0a05744ad804ef9fd0671ef Mon Sep 17 00:00:00 2001 From: myerspat Date: Thu, 9 Jun 2022 10:39:46 -0500 Subject: [PATCH 4/8] Moved initialized variables into if statement --- src/physics.cpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/physics.cpp b/src/physics.cpp index 26b3a3520a..0df8e63c0a 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -345,12 +345,14 @@ 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) { + + 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); + for (int i_shell = 0; i_shell < element.shells_.size(); ++i_shell) { const auto& shell {element.shells_[i_shell]}; From 98a7d5eabc972a986a64307bc709dda7ecadaf3d Mon Sep 17 00:00:00 2001 From: myerspat Date: Fri, 10 Jun 2022 11:57:35 -0500 Subject: [PATCH 5/8] added clang format --- src/photon.cpp | 7 ++++--- src/physics.cpp | 5 ++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/photon.cpp b/src/photon.cpp index bb715cdeaf..85bd435c99 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -165,8 +165,8 @@ PhotonInteraction::PhotonInteraction(hid_t group) close_dataset(dset); read_dataset(tgroup, "xs", xs); - auto cross_section = xt::view(cross_sections_, - xt::range(shell.threshold, xt::placeholders::_), i); + auto cross_section = xt::view( + cross_sections_, xt::range(shell.threshold, xt::placeholders::_), i); cross_section = xt::where(xs > 0, xt::log(xs), 0); if (object_exists(tgroup, "transitions")) { @@ -577,7 +577,8 @@ void PhotonInteraction::calculate_xs(Particle& p) const for (int i = 0; i < xs_upper.size(); ++i) if (xs_lower(i) != 0) - xs.photoelectric += std::exp(xs_lower(i) + f * (xs_upper(i) - xs_lower(i))); + xs.photoelectric += + std::exp(xs_lower(i) + f * (xs_upper(i) - xs_lower(i))); // Calculate microscopic pair production cross section xs.pair_production = std::exp( diff --git a/src/physics.cpp b/src/physics.cpp index 0df8e63c0a..cc712affd7 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -361,9 +361,8 @@ void sample_photon_reaction(Particle& p) continue; // Evaluation subshell photoionization cross section - prob += std::exp(xs_lower(i_shell) + - f * (xs_upper(i_shell) - - xs_lower(i_shell))); + prob += std::exp( + xs_lower(i_shell) + f * (xs_upper(i_shell) - xs_lower(i_shell))); if (prob > cutoff) { double E_electron = p.E() - shell.binding_energy; From d80bf7bdb596a7169f0e28d42b363d1886c43448 Mon Sep 17 00:00:00 2001 From: Patrick Myers <90068356+myerspat@users.noreply.github.com> Date: Wed, 29 Jun 2022 09:25:59 -0500 Subject: [PATCH 6/8] Added comments explaining reasoning for 2D xtensor xs in photon.h Co-authored-by: Paul Romano --- include/openmc/photon.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/include/openmc/photon.h b/include/openmc/photon.h index 50ee228c44..09fb3ba01b 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -79,7 +79,9 @@ public: Tabulated1D coherent_anomalous_real_; Tabulated1D coherent_anomalous_imag_; - // Photoionization and atomic relaxation data + // Photoionization and atomic relaxation data. Subshell cross sections are + // stored separately to improve memory access pattern when calculating the + // total cross section vector shells_; xt::xtensor cross_sections_; From 2c5a845897aeb973448f31cea2d9256754262e3e Mon Sep 17 00:00:00 2001 From: Patrick Myers <90068356+myerspat@users.noreply.github.com> Date: Wed, 29 Jun 2022 09:26:36 -0500 Subject: [PATCH 7/8] Removed an include from photon.cpp Co-authored-by: Paul Romano --- src/photon.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/photon.cpp b/src/photon.cpp index 85bd435c99..a127ac2fee 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -17,7 +17,6 @@ #include "xtensor/xmath.hpp" #include "xtensor/xoperation.hpp" #include "xtensor/xslice.hpp" -#include "xtensor/xtensor_forward.hpp" #include "xtensor/xview.hpp" #include From c38ed4d8759704e829454abfc12dfa2189dd4bf1 Mon Sep 17 00:00:00 2001 From: myerspat Date: Wed, 29 Jun 2022 09:37:18 -0500 Subject: [PATCH 8/8] used placeholders namespace and added comment about variables initialized --- src/photon.cpp | 6 ++++-- src/physics.cpp | 3 ++- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/photon.cpp b/src/photon.cpp index a127ac2fee..2f9bc5e3dd 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -46,6 +46,8 @@ vector> elements; PhotonInteraction::PhotonInteraction(hid_t group) { + using namespace xt::placeholders; + // Set index of element in global vector index_ = data::elements.size(); @@ -164,8 +166,8 @@ PhotonInteraction::PhotonInteraction(hid_t group) close_dataset(dset); read_dataset(tgroup, "xs", xs); - auto cross_section = xt::view( - cross_sections_, xt::range(shell.threshold, xt::placeholders::_), i); + auto cross_section = + xt::view(cross_sections_, xt::range(shell.threshold, _), i); cross_section = xt::where(xs > 0, xt::log(xs), 0); if (object_exists(tgroup, "transitions")) { diff --git a/src/physics.cpp b/src/physics.cpp index cc712affd7..bcdcf7ecf9 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -347,7 +347,8 @@ void sample_photon_reaction(Particle& p) double prob_after = prob + micro.photoelectric; if (prob_after > cutoff) { - + // Get grid index, interpolation factor, and bounding subshell + // cross sections int i_grid = micro.index_grid; double f = micro.interp_factor; const auto& xs_lower = xt::row(element.cross_sections_, i_grid);