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

Move results and global_tallies over to C++

Browse source
This commit is contained in:
Paul Romano 2019-02-20 07:33:06 -06:00
parent 1542cf7bfd
commit 3eea8311d3
19 changed files with 315 additions and 613 deletions
Download patch file Download diff file Expand all files Collapse all files

1
CMakeLists.txt
View file

@ -319,7 +319,6 @@ add_library(libopenmc SHARED
src/relaxng
src/settings.F90
src/simulation_header.F90
src/state_point.F90
src/stl_vector.F90
src/string.F90
src/vector_header.F90

1
include/openmc/capi.h
View file

@ -164,7 +164,6 @@ extern "C" {
// Global variables
extern char openmc_err_msg[256];
extern int32_t n_realizations;
extern int32_t n_sab_tables;
extern int32_t n_tallies;

2
include/openmc/simulation.h
View file

@ -36,7 +36,7 @@ extern "C" bool need_depletion_rx; //!< need to calculate depletion rx?
extern "C" int restart_batch; //!< batch at which a restart job resumed
extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied?
extern "C" int total_gen; //!< total number of generations simulated
extern "C" double total_weight; //!< Total source weight in a batch
extern double total_weight; //!< Total source weight in a batch
extern "C" int64_t work; //!< number of particles per process
extern std::vector<double> k_generation;

55
include/openmc/tallies/tally.h
View file

@ -5,6 +5,7 @@
#include "openmc/tallies/trigger.h"
#include "pugixml.hpp"
#include "xtensor/xfixed.hpp"
#include "xtensor/xtensor.hpp"
#include <memory> // for unique_ptr
@ -47,6 +48,10 @@ public:
void init_triggers(pugi::xml_node node, int i_tally);
void init_results();
void accumulate();
//----------------------------------------------------------------------------
// Major public data members.
@ -62,6 +67,9 @@ public:
//! Whether this tally is currently being updated
bool active_ {false};
//! Number of realizations
int n_realizations_ {0};
std::vector<int> scores_; //!< Filter integrands (e.g. flux, fission)
//! Index of each nuclide to be tallied. -1 indicates total material.
@ -70,6 +78,12 @@ public:
//! True if this tally has a bin for every nuclide in the problem
bool all_nuclides_ {false};
//! Results for each bin -- the first dimension of the array is for scores
//! (e.g. flux, total reaction rate, fission reaction rate, etc.) and the
//! second dimension of the array is for the combination of filters
//! (e.g. specific cell, specific energy group, etc.)
xt::xtensor<double, 3> results_;
//----------------------------------------------------------------------------
// Miscellaneous public members.
@ -102,17 +116,22 @@ private:
//==============================================================================
namespace model {
extern std::vector<std::unique_ptr<Tally>> tallies;
extern std::vector<int> active_tallies;
extern std::vector<int> active_analog_tallies;
extern std::vector<int> active_tracklength_tallies;
extern std::vector<int> active_collision_tallies;
extern std::vector<int> active_meshsurf_tallies;
extern std::vector<int> active_surface_tallies;
}
extern std::vector<std::unique_ptr<Tally>> tallies;
namespace simulation {
//! Global tallies (such as k-effective estimators)
extern xt::xtensor_fixed<double, xt::xshape<N_GLOBAL_TALLIES, 3>> global_tallies;
extern std::vector<int> active_tallies;
extern std::vector<int> active_analog_tallies;
extern std::vector<int> active_tracklength_tallies;
extern std::vector<int> active_collision_tallies;
extern std::vector<int> active_meshsurf_tallies;
extern std::vector<int> active_surface_tallies;
} // namespace model
//! Number of realizations for global tallies
extern int32_t n_realizations;
}
// It is possible to protect accumulate operations on global tallies by using an
// atomic update. However, when multiple threads accumulate to the same global
@ -130,23 +149,21 @@ extern double global_tally_leakage;
// Non-member functions
//==============================================================================
//! \brief Accumulate the sum of the contributions from each history within the
//! batch to a new random variable
void accumulate_tallies();
//! Determine which tallies should be active
void setup_active_tallies();
// Alias for the type returned by xt::adapt(...). N is the dimension of the
// multidimensional array
template <std::size_t N>
using adaptor_type = xt::xtensor_adaptor<xt::xbuffer_adaptor<double*&, xt::no_ownership>, N>;
//! Get the global tallies as a multidimensional array
//! \return Global tallies array
adaptor_type<2> global_tallies();
//! Get tally results as a multidimensional array
//! \param idx Index in tallies array
//! \return Tally results array
adaptor_type<3> tally_results(int idx);
#ifdef OPENMC_MPI
//! Collect all tally results onto master process
extern "C" void reduce_tally_results();
void reduce_tally_results();
#endif
extern "C" void free_memory_tally_c();

8
src/constants.F90
View file

@ -127,14 +127,6 @@ module constants
FILTER_MATERIAL = 2, &
FILTER_CELL = 3
! Global tally parameters
integer, parameter :: N_GLOBAL_TALLIES = 4
integer, parameter :: &
K_COLLISION = 1, &
K_ABSORPTION = 2, &
K_TRACKLENGTH = 3, &
LEAKAGE = 4
! Differential tally independent variables
integer, parameter :: &
DIFF_DENSITY = 1, &

10
src/eigenvalue.cpp
View file

@ -46,7 +46,7 @@ xt::xtensor<double, 1> source_frac;
void calculate_generation_keff()
{
auto gt = global_tallies();
const auto& gt = simulation::global_tallies;
// Get keff for this generation by subtracting off the starting value
simulation::keff_generation = gt(K_TRACKLENGTH, RESULT_VALUE) - simulation::keff_generation;
@ -309,7 +309,7 @@ void calculate_average_keff()
int i = overall_generation() - 1;
int n;
if (simulation::current_batch > settings::n_inactive) {
n = settings::gen_per_batch*n_realizations + simulation::current_gen;
n = settings::gen_per_batch*simulation::n_realizations + simulation::current_gen;
} else {
n = 0;
}
@ -390,15 +390,15 @@ int openmc_get_keff(double* k_combined)
// Make sure we have at least four realizations. Notice that at the end,
// there is a N-3 term in a denominator.
if (n_realizations <= 3) {
if (simulation::n_realizations <= 3) {
return -1;
}
// Initialize variables
int64_t n = n_realizations;
int64_t n = simulation::n_realizations;
// Copy estimates of k-effective and its variance (not variance of the mean)
auto gt = global_tallies();
const auto& gt = simulation::global_tallies;
std::array<double, 3> kv {};
xt::xtensor<double, 2> cov = xt::zeros<double>({3, 3});

16
src/finalize.cpp
View file

@ -12,6 +12,8 @@
#include "openmc/timer.h"
#include "openmc/tallies/tally.h"
#include "xtensor/xview.hpp"
using namespace openmc;
// Functions defined in Fortran
@ -100,17 +102,9 @@ int openmc_reset()
openmc_tally_reset(i);
}
// Reset global tallies (can't really use global_tallies() right now because
// it doesn't have any information about whether the underlying buffer was
// allocated)
n_realizations = 0;
double* buffer = nullptr;
openmc_global_tallies(&buffer);
if (buffer) {
for (int i = 0; i < 3*N_GLOBAL_TALLIES; ++i) {
buffer[i] = 0.0;
}
}
// Reset global tallies
simulation::n_realizations = 0;
xt::view(simulation::global_tallies, xt::all()) = 0.0;
simulation::k_col_abs = 0.0;
simulation::k_col_tra = 0.0;

16
src/output.cpp
View file

@ -370,7 +370,7 @@ void print_generation()
// Determine overall generation and number of active generations
int i = overall_generation() - 1;
int n = simulation::current_batch > settings::n_inactive ?
settings::gen_per_batch*n_realizations + simulation::current_gen : 0;
settings::gen_per_batch*simulation::n_realizations + simulation::current_gen : 0;
// Set format for values
std::cout << std::fixed << std::setprecision(5);
@ -404,7 +404,7 @@ void print_batch_keff()
// Determine overall generation and number of active generations
int i = simulation::current_batch*settings::gen_per_batch - 1;
int n = n_realizations*settings::gen_per_batch;
int n = simulation::n_realizations*settings::gen_per_batch;
// Set format for values
std::cout << std::fixed << std::setprecision(5);
@ -538,7 +538,7 @@ void print_results()
header("Results", 4);
// Calculate t-value for confidence intervals
int n = n_realizations;
int n = simulation::n_realizations;
double alpha, t_n1, t_n3;
if (settings::confidence_intervals) {
alpha = 1.0 - CONFIDENCE_LEVEL;
@ -553,7 +553,7 @@ void print_results()
std::cout << std::fixed << std::setprecision(5);
// write global tallies
auto gt = global_tallies();
const auto& gt = simulation::global_tallies;
double mean, stdev;
if (n > 1) {
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
@ -633,16 +633,12 @@ write_tallies()
// Loop over each tally.
for (auto i_tally = 0; i_tally < model::tallies.size(); ++i_tally) {
const auto& tally {*model::tallies[i_tally]};
auto results = tally_results(i_tally);
// TODO: get this directly from the tally object when it's been translated
int32_t n_realizations;
auto err = openmc_tally_get_n_realizations(i_tally+1, &n_realizations);
// Calculate t-value for confidence intervals
double t_value = 1;
if (settings::confidence_intervals) {
auto alpha = 1 - CONFIDENCE_LEVEL;
t_value = t_percentile(1 - alpha*0.5, n_realizations - 1);
t_value = t_percentile(1 - alpha*0.5, tally.n_realizations_ - 1);
}
// Write header block.
@ -719,7 +715,7 @@ write_tallies()
double mean, stdev;
//TODO: off-by-one
std::tie(mean, stdev) = mean_stdev(
&results(filter_index-1, score_index, 0), n_realizations);
&tally.results_(filter_index-1, score_index, 0), tally.n_realizations_);
tallies_out << std::string(indent+1, ' ') << std::left
<< std::setw(36) << score_name << " " << mean << " +/- "
<< t_value * stdev << "\n";

3
src/photon.cpp
View file

@ -176,7 +176,7 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element)
read_dataset(rgroup, "J", profile_pdf_);
// Get Compton profile momentum grid. By deafult, an xtensor has a size of 1.
// TODO: Update version of xtensor and change to 0
// TODO: Change to zero when xtensor is updated
if (data::compton_profile_pz.size() == 1) {
read_dataset(rgroup, "pz", data::compton_profile_pz);
}
@ -208,6 +208,7 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element)
// Get energy grids used for bremsstrahlung DCS and for stopping powers
xt::xtensor<double, 1> electron_energy;
read_dataset(rgroup, "electron_energy", electron_energy);
// TODO: Change to zero when xtensor is updated
if (data::ttb_k_grid.size() == 1) {
read_dataset(rgroup, "photon_energy", data::ttb_k_grid);
}

23
src/simulation.cpp
View file

@ -21,6 +21,7 @@
#include "openmc/tallies/trigger.h"
#include <omp.h>
#include "xtensor/xview.hpp"
#include <algorithm>
#include <string>
@ -28,9 +29,6 @@
namespace openmc {
// data/functions from Fortran side
extern "C" void accumulate_tallies();
extern "C" void allocate_tally_results();
extern "C" void setup_active_tallies();
extern "C" void write_tallies();
} // namespace openmc
@ -78,7 +76,9 @@ int openmc_simulation_init()
allocate_banks();
// Allocate tally results arrays if they're not allocated yet
allocate_tally_results();
for (auto& t : model::tallies) {
t->init_results();
}
// Set up material nuclide index mapping
for (auto& mat : model::materials) {
@ -263,6 +263,7 @@ bool need_depletion_rx {false};
int restart_batch;
bool satisfy_triggers {false};
int total_gen {0};
double total_weight;
int64_t work;
std::vector<double> k_generation;
@ -360,9 +361,8 @@ void finalize_batch()
// Reset global tally results
if (simulation::current_batch <= settings::n_inactive) {
auto gt = global_tallies();
std::fill(gt.begin(), gt.end(), 0.0);
n_realizations = 0;
xt::view(simulation::global_tallies, xt::all()) = 0.0;
simulation::n_realizations = 0;
}
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
@ -415,13 +415,14 @@ void initialize_generation()
if (settings::ufs_on) ufs_count_sites();
// Store current value of tracklength k
simulation::keff_generation = global_tallies()(K_TRACKLENGTH, RESULT_VALUE);
simulation::keff_generation = simulation::global_tallies(
K_TRACKLENGTH, RESULT_VALUE);
}
}
void finalize_generation()
{
auto gt = global_tallies();
auto& gt = simulation::global_tallies;
// Update global tallies with the omp private accumulation variables
#pragma omp parallel
@ -547,7 +548,7 @@ void broadcast_results() {
// Create a new datatype that consists of all values for a given filter
// bin and then use that to broadcast. This is done to minimize the
// chance of the 'count' argument of MPI_BCAST exceeding 2**31
auto results = tally_results(i);
auto& results = model::tallies[i]->results_;
auto shape = results.shape();
int count_per_filter = shape[1] * shape[2];
@ -559,7 +560,7 @@ void broadcast_results() {
}
// Also broadcast global tally results
auto gt = global_tallies();
auto& gt = simulation::global_tallies;
MPI_Bcast(gt.data(), gt.size(), MPI_DOUBLE, 0, mpi::intracomm);
// These guys are needed so that non-master processes can calculate the

3
src/simulation_header.F90
View file

@ -30,9 +30,6 @@ module simulation_header
! Temporary k-effective values
real(C_DOUBLE), bind(C) :: keff ! average k over active batches
real(C_DOUBLE), bind(C) :: keff_std ! standard deviation of average k
real(C_DOUBLE), bind(C) :: k_col_abs ! sum over batches of k_collision * k_absorption
real(C_DOUBLE), bind(C) :: k_col_tra ! sum over batches of k_collision * k_tracklength
real(C_DOUBLE), bind(C) :: k_abs_tra ! sum over batches of k_absorption * k_tracklength
! ============================================================================
! PARALLEL PROCESSING VARIABLES

9
src/source.cpp
View file

@ -336,15 +336,6 @@ extern "C" void free_memory_source()
model::external_sources.clear();
}
extern "C" double total_source_strength()
{
double strength = 0.0;
for (const auto& s : model::external_sources) {
strength += s.strength();
}
return strength;
}
void fill_source_bank_fixedsource()
{
if (settings::path_source.empty()) {

105
src/state_point.F90
View file

@ -1,105 +0,0 @@
module state_point
!===============================================================================
! STATE_POINT -- This module handles writing and reading state point
! files. State points are contain complete tally results, source sites, and
! various other data. They can be used to restart a run or to reconstruct
! confidence intervals for tallies (this requires post-processing via Python
! scripts).
!
! State points can be written at any batch during a simulation, or at specified
! intervals, using the <state_point ... /> tag.
!===============================================================================
use, intrinsic :: ISO_C_BINDING
use hdf5_interface
use settings
use string, only: to_str
use tally_header
implicit none
contains
!===============================================================================
! OPENMC_STATEPOINT_WRITE writes an HDF5 statepoint file to disk
!===============================================================================
subroutine statepoint_write_f(file_id) bind(C)
integer(HID_T), value :: file_id
integer :: i
integer(HID_T) :: tallies_group, tally_group
! Open tallies group
tallies_group = open_group(file_id, "tallies")
if (reduce_tallies) then
! Write global tallies
call write_dataset(file_id, "global_tallies", global_tallies)
! Write tallies
if (active_tallies_size() > 0) then
! Indicate that tallies are on
call write_attribute(file_id, "tallies_present", 1)
! Write all tally results
TALLY_RESULTS: do i = 1, n_tallies
associate (tally => tallies(i) % obj)
! Write sum and sum_sq for each bin
tally_group = open_group(tallies_group, "tally " &
// to_str(tally % id()))
call tally % write_results_hdf5(tally_group)
call close_group(tally_group)
end associate
end do TALLY_RESULTS
else
! Indicate tallies are off
call write_attribute(file_id, "tallies_present", 0)
end if
end if
call close_group(tallies_group)
end subroutine
!===============================================================================
! LOAD_STATE_POINT
!===============================================================================
subroutine load_state_point_f(file_id) bind(C)
integer(HID_T), value :: file_id
integer :: i
integer :: temp
integer(HID_T) :: tallies_group
integer(HID_T) :: tally_group
! Read global tally data
call read_dataset(global_tallies, file_id, "global_tallies")
! Check if tally results are present
call read_attribute(temp, file_id, "tallies_present")
! Read in sum and sum squared
if (temp == 1) then
tallies_group = open_group(file_id, "tallies")
TALLY_RESULTS: do i = 1, n_tallies
associate (t => tallies(i) % obj)
! Read sum, sum_sq, and N for each bin
tally_group = open_group(tallies_group, "tally " // &
trim(to_str(t % id())))
call t % read_results_hdf5(tally_group)
call read_dataset(t % n_realizations, tally_group, &
"n_realizations")
call close_group(tally_group)
end associate
end do TALLY_RESULTS
call close_group(tallies_group)
end if
end subroutine load_state_point_f
end module state_point

95
src/state_point.cpp
View file

@ -169,11 +169,8 @@ openmc_statepoint_write(const char* filename, bool* write_source)
write_attribute(tallies_group, "ids", tally_ids);
// Write all tally information except results
//TODO: use these two lines when openmc_get_n_realizations isn't needed
//for (const auto& tally_ptr : model::tallies) {
// const auto& tally {*tally_ptr};
for (auto i_tally = 0; i_tally < model::tallies.size(); ++i_tally) {
const auto& tally {*model::tallies[i_tally]};
for (const auto& tally_ptr : model::tallies) {
const auto& tally {*tally_ptr};
hid_t tally_group = create_group(tallies_group,
"tally " + std::to_string(tally.id_));
@ -187,10 +184,7 @@ openmc_statepoint_write(const char* filename, bool* write_source)
write_dataset(tally_group, "estimator", "collision");
}
// TODO: get this directly from the tally object when it's moved to C++
int32_t n_realizations;
auto err = openmc_tally_get_n_realizations(i_tally+1, &n_realizations);
write_dataset(tally_group, "n_realizations", n_realizations);
write_dataset(tally_group, "n_realizations", tally.n_realizations_);
// Write the ID of each filter attached to this tally
write_dataset(tally_group, "n_filters", tally.filters().size());
@ -231,8 +225,33 @@ openmc_statepoint_write(const char* filename, bool* write_source)
}
if (settings::reduce_tallies) {
// Write global tallies
write_dataset(file_id, "global_tallies", simulation::global_tallies);
// Write tallies
if (model::active_tallies.size() > 0) {
// Indicate that tallies are on
write_attribute(file_id, "tallies_present", 1);
// Write all tally results
for (int i = 0; i < model::tallies.size(); ++i) {
// Write sum and sum_sq for each bin
const auto& tally {model::tallies[i]};
std::string name = "tally " + std::to_string(tally->id_);
hid_t tally_group = open_group(tallies_group, name.c_str());
auto& results = tally->results_;
write_tally_results(tally_group, results.shape()[0],
results.shape()[1], results.data());
close_group(tally_group);
}
} else {
// Indicate tallies are off
write_attribute(file_id, "tallies_present", 0);
}
}
close_group(tallies_group);
statepoint_write_f(file_id);
}
// Check for the no-tally-reduction method
@ -243,7 +262,7 @@ openmc_statepoint_write(const char* filename, bool* write_source)
} else if (mpi::master) {
// Write number of global realizations
write_dataset(file_id, "n_realizations", n_realizations);
write_dataset(file_id, "n_realizations", simulation::n_realizations);
// Write out the runtime metrics.
using namespace simulation;
@ -292,7 +311,7 @@ void restart_set_keff()
simulation::k_sum[0] += simulation::k_generation[i];
simulation::k_sum[1] += std::pow(simulation::k_generation[i], 2);
}
int n = settings::gen_per_batch*n_realizations;
int n = settings::gen_per_batch*simulation::n_realizations;
simulation::keff = simulation::k_sum[0] / n;
} else {
simulation::keff = simulation::k_generation.back();
@ -375,7 +394,7 @@ void load_state_point()
}
// Read number of realizations for global tallies
read_dataset(file_id, "n_realizations", n_realizations);
read_dataset(file_id, "n_realizations", simulation::n_realizations);
// Set k_sum, keff, and current_batch based on whether restart file is part
// of active cycle or inactive cycle
@ -395,8 +414,32 @@ void load_state_point()
#else
if (mpi::master) {
#endif
// Read tally results
load_state_point_f(file_id);
// Read global tally data
read_dataset(file_id, "global_tallies", H5T_NATIVE_DOUBLE,
simulation::global_tallies.data());
// Check if tally results are present
bool present;
read_attribute(file_id, "tallies_present", present);
// Read in sum and sum squared
if (present) {
hid_t tallies_group = open_group(file_id, "tallies");
for (int i = 0; i < model::tallies.size(); ++i) {
// Read sum, sum_sq, and N for each bin
auto& tally {model::tallies[i]};
std::string name = "tally " + std::to_string(tally->id_);
hid_t tally_group = open_group(tallies_group, name.c_str());
auto& results = tally->results_;
read_tally_results(tally_group, results.shape()[0],
results.shape()[1], results.data());
read_dataset(tally_group, "n_realizations", tally->n_realizations_);
close_group(tally_group);
}
close_group(tallies_group);
}
}
// Read source if in eigenvalue mode
@ -623,7 +666,7 @@ void write_tally_results_nr(hid_t file_id)
hid_t tallies_group;
if (mpi::master) {
// Write number of realizations
write_dataset(file_id, "n_realizations", n_realizations);
write_dataset(file_id, "n_realizations", simulation::n_realizations);
// Write number of global tallies
write_dataset(file_id, "n_global_tallies", N_GLOBAL_TALLIES);
@ -631,8 +674,8 @@ void write_tally_results_nr(hid_t file_id)
tallies_group = open_group(file_id, "tallies");
}
// Get pointer to global tallies
auto gt = global_tallies();
// Get global tallies
auto& gt = simulation::global_tallies;
#ifdef OPENMC_MPI
// Reduce global tallies
@ -656,18 +699,15 @@ void write_tally_results_nr(hid_t file_id)
for (int i = 0; i < n_tallies; ++i) {
// Skip any tallies that are not active
bool active;
// TODO: off-by-one
openmc_tally_get_active(i+1, &active);
if (!active) continue;
const auto& t {model::tallies[i]};
if (!t->active_) continue;
if (mpi::master && !object_exists(file_id, "tallies_present")) {
write_attribute(file_id, "tallies_present", 1);
}
// Get view of accumulated tally values
auto results = tally_results(i);
auto values_view = xt::view(results, xt::all(), xt::all(),
auto values_view = xt::view(t->results_, xt::all(), xt::all(),
xt::range(RESULT_SUM, RESULT_SUM_SQ + 1));
// Make copy of tally values in contiguous array
@ -675,10 +715,7 @@ void write_tally_results_nr(hid_t file_id)
if (mpi::master) {
// Open group for tally
int id;
// TODO: off-by-one
openmc_tally_get_id(i+1, &id);
std::string groupname {"tally " + std::to_string(id)};
std::string groupname {"tally " + std::to_string(t->id_)};
hid_t tally_group = open_group(tallies_group, groupname.c_str());
// The MPI_IN_PLACE specifier allows the master to copy values into
@ -696,7 +733,7 @@ void write_tally_results_nr(hid_t file_id)
}
// Put in temporary tally result
xt::xtensor<double, 3> results_copy = xt::zeros_like(results);
xt::xtensor<double, 3> results_copy = xt::zeros_like(t->results_);
auto copy_view = xt::view(results_copy, xt::all(), xt::all(),
xt::range(RESULT_SUM, RESULT_SUM_SQ + 1));
copy_view = values;

91
src/tallies/tally.F90
View file

@ -16,97 +16,6 @@ module tally
contains
!===============================================================================
! ACCUMULATE_TALLIES accumulates the sum of the contributions from each history
! within the batch to a new random variable
!===============================================================================
subroutine accumulate_tallies() bind(C)
integer :: i
real(C_DOUBLE) :: k_col ! Copy of batch collision estimate of keff
real(C_DOUBLE) :: k_abs ! Copy of batch absorption estimate of keff
real(C_DOUBLE) :: k_tra ! Copy of batch tracklength estimate of keff
real(C_DOUBLE) :: val
#ifdef OPENMC_MPI
interface
subroutine reduce_tally_results() bind(C)
end subroutine
end interface
! Combine tally results onto master process
if (reduce_tallies) call reduce_tally_results()
#endif
! Increase number of realizations (only used for global tallies)
if (reduce_tallies) then
n_realizations = n_realizations + 1
else
n_realizations = n_realizations + n_procs
end if
! Accumulate on master only unless run is not reduced then do it on all
if (master .or. (.not. reduce_tallies)) then
if (run_mode == MODE_EIGENVALUE) then
if (current_batch > n_inactive) then
! Accumulate products of different estimators of k
k_col = global_tallies(RESULT_VALUE, K_COLLISION) / total_weight
k_abs = global_tallies(RESULT_VALUE, K_ABSORPTION) / total_weight
k_tra = global_tallies(RESULT_VALUE, K_TRACKLENGTH) / total_weight
k_col_abs = k_col_abs + k_col * k_abs
k_col_tra = k_col_tra + k_col * k_tra
k_abs_tra = k_abs_tra + k_abs * k_tra
end if
end if
! Accumulate results for global tallies
do i = 1, size(global_tallies, 2)
val = global_tallies(RESULT_VALUE, i)/total_weight
global_tallies(RESULT_VALUE, i) = ZERO
global_tallies(RESULT_SUM, i) = global_tallies(RESULT_SUM, i) + val
global_tallies(RESULT_SUM_SQ, i) = &
global_tallies(RESULT_SUM_SQ, i) + val*val
end do
end if
! Accumulate results for each tally
do i = 1, active_tallies_size()
call tallies(active_tallies_data(i)+1) % obj % accumulate()
end do
end subroutine accumulate_tallies
!===============================================================================
! SETUP_ACTIVE_TALLIES
!===============================================================================
subroutine setup_active_tallies() bind(C)
integer :: i
integer(C_INT) :: err
logical(C_BOOL) :: active
interface
subroutine setup_active_tallies_c() bind(C)
end subroutine
end interface
call setup_active_tallies_c()
do i = 1, n_tallies
associate (t => tallies(i) % obj)
err = openmc_tally_get_active(i, active)
if (active) then
! Check if tally contains depletion reactions and if so, set flag
if (t % depletion_rx()) need_depletion_rx = .true.
end if
end associate
end do
end subroutine setup_active_tallies
!===============================================================================
! C API FUNCTIONS
!===============================================================================

201
src/tallies/tally.cpp
View file

@ -9,6 +9,7 @@
#include "openmc/reaction_product.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/source.h"
#include "openmc/tallies/derivative.h"
#include "openmc/tallies/filter.h"
#include "openmc/tallies/filter_cell.h"
@ -38,7 +39,6 @@ namespace openmc {
namespace model {
std::vector<std::unique_ptr<Tally>> tallies;
std::vector<int> active_tallies;
std::vector<int> active_analog_tallies;
std::vector<int> active_tracklength_tallies;
@ -47,6 +47,11 @@ namespace model {
std::vector<int> active_surface_tallies;
}
namespace simulation {
xt::xtensor_fixed<double, xt::xshape<N_GLOBAL_TALLIES, 3>> global_tallies;
int32_t n_realizations {0};
}
double global_tally_absorption;
double global_tally_collision;
double global_tally_tracklength;
@ -503,49 +508,54 @@ Tally::init_triggers(pugi::xml_node node, int i_tally)
}
}
void Tally::init_results()
{
int n_scores = scores_.size() * nuclides_.size();
results_ = xt::empty<double>({n_filter_bins_, n_scores, 3});
}
void Tally::accumulate()
{
// Increment number of realizations
n_realizations_ += settings::reduce_tallies ? 1 : mpi::n_procs;
if (mpi::master || !settings::reduce_tallies) {
// Calculate total source strength for normalization
double total_source = 0.0;
if (settings::run_mode == RUN_MODE_FIXEDSOURCE) {
for (const auto& s : model::external_sources) {
total_source += s.strength();
}
} else {
total_source = 1.0;
}
// Accumulate each result
for (int i = 0; i < results_.shape()[0]; ++i) {
for (int j = 0; j < results_.shape()[1]; ++j) {
double val = results_(i, j, RESULT_VALUE) /
simulation::total_weight * total_source;
results_(i, j, RESULT_VALUE) = 0.0;
results_(i, j, RESULT_SUM) += val;
results_(i, j, RESULT_SUM_SQ) += val*val;
}
}
}
}
//==============================================================================
// Non-member functions
//==============================================================================
adaptor_type<2> global_tallies()
{
// Get pointer to global tallies
double* buffer;
openmc_global_tallies(&buffer);
// Adapt into xtensor
std::array<size_t, 2> shape = {N_GLOBAL_TALLIES, 3};
std::size_t size {3*N_GLOBAL_TALLIES};
return xt::adapt(buffer, size, xt::no_ownership(), shape);
}
adaptor_type<3> tally_results(int idx)
{
// Get pointer to tally results
double* results;
std::array<std::size_t, 3> shape;
// TODO: off-by-one
openmc_tally_results(idx+1, &results, shape.data());
// Adapt array into xtensor with no ownership
std::size_t size {shape[0] * shape[1] * shape[2]};
return xt::adapt(results, size, xt::no_ownership(), shape);
}
#ifdef OPENMC_MPI
void reduce_tally_results()
{
for (int i = 0; i < n_tallies; ++i) {
for (int i_tally : model::active_tallies) {
// Skip any tallies that are not active
bool active;
// TODO: off-by-one
openmc_tally_get_active(i+1, &active);
if (!active) continue;
auto& tally {model::tallies[i_tally]};
// Get view of accumulated tally values
auto results = tally_results(i);
auto values_view = xt::view(results, xt::all(), xt::all(), RESULT_VALUE);
auto values_view = xt::view(tally->results_, xt::all(), xt::all(), RESULT_VALUE);
// Make copy of tally values in contiguous array
xt::xtensor<double, 2> values = values_view;
@ -564,7 +574,7 @@ void reduce_tally_results()
}
// Get view of global tally values
auto gt = global_tallies();
auto& gt = simulation::global_tallies;
auto gt_values_view = xt::view(gt, xt::all(), RESULT_VALUE);
// Make copy of values in contiguous array
@ -591,8 +601,51 @@ void reduce_tally_results()
}
#endif
extern "C" void
setup_active_tallies_c()
void
accumulate_tallies()
{
#ifdef OPENMC_MPI
// Combine tally results onto master process
if (settings::reduce_tallies) reduce_tally_results();
#endif
// Increase number of realizations (only used for global tallies)
simulation::n_realizations += settings::reduce_tallies ? 1 : mpi::n_procs;
// Accumulate on master only unless run is not reduced then do it on all
if (mpi::master || !settings::reduce_tallies) {
auto& gt = simulation::global_tallies;
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
if (simulation::current_batch > settings::n_inactive) {
// Accumulate products of different estimators of k
double k_col = gt(K_COLLISION, RESULT_VALUE) / simulation::total_weight;
double k_abs = gt(K_ABSORPTION, RESULT_VALUE) / simulation::total_weight;
double k_tra = gt(K_TRACKLENGTH, RESULT_VALUE) / simulation::total_weight;
simulation::k_col_abs += k_col * k_abs;
simulation::k_col_tra += k_col * k_tra;
simulation::k_abs_tra += k_abs * k_tra;
}
}
// Accumulate results for global tallies
for (int i = 0; i < N_GLOBAL_TALLIES; ++i) {
double val = gt(i, RESULT_VALUE)/simulation::total_weight;
gt(i, RESULT_VALUE) = 0.0;
gt(i, RESULT_SUM) += val;
gt(i, RESULT_SUM_SQ) += val*val;
}
}
// Accumulate results for each tally
for (int i_tally : model::active_tallies) {
auto& tally {model::tallies[i_tally]};
tally->accumulate();
}
}
void
setup_active_tallies()
{
model::active_tallies.clear();
model::active_analog_tallies.clear();
@ -628,6 +681,9 @@ setup_active_tallies_c()
case TALLY_SURFACE:
model::active_surface_tallies.push_back(i);
}
// Check if tally contains depletion reactions and if so, set flag
if (tally.depletion_rx_) simulation::need_depletion_rx = true;
}
}
}
@ -833,6 +889,77 @@ openmc_tally_set_filters(int32_t index, int n, const int32_t* indices)
return 0;
}
//! Reset tally results and number of realizations
extern "C" int
openmc_tally_reset(int32_t index)
{
// Make sure the index fits in the array bounds.
if (index < 1 || index > model::tallies.size()) {
set_errmsg("Index in tallies array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
// TODO: off-by-one
auto& t {model::tallies[index-1]};
t->n_realizations_ = 0;
// TODO: Change to zero when xtensor is updated
if (t->results_.size() != 1) {
xt::view(t->results_, xt::all()) = 0.0;
}
return 0;
}
extern "C" int
openmc_tally_get_n_realizations(int32_t index, int32_t* n)
{
// Make sure the index fits in the array bounds.
if (index < 1 || index > model::tallies.size()) {
set_errmsg("Index in tallies array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
// TODO: off-by-one
*n = model::tallies[index - 1]->n_realizations_;
return 0;
}
//! \brief Returns a pointer to a tally results array along with its shape. This
//! allows a user to obtain in-memory tally results from Python directly.
extern "C" int
openmc_tally_results(int32_t index, double** results, size_t* shape)
{
// Make sure the index fits in the array bounds.
if (index < 1 || index > model::tallies.size()) {
set_errmsg("Index in tallies array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
// TODO: off-by-one
const auto& t {model::tallies[index - 1]};
// TODO: Change to zero when xtensor is updated
if (t->results_.size() == 1) {
set_errmsg("Tally results have not been allocated yet.");
return OPENMC_E_ALLOCATE;
}
// Set pointer to results and copy shape
*results = t->results_.data();
auto s = t->results_.shape();
shape[0] = s[0];
shape[1] = s[1];
shape[2] = s[2];
return 0;
}
extern "C" int
openmc_global_tallies(double** ptr)
{
*ptr = simulation::global_tallies.data();
return 0;
}
extern "C" size_t tallies_size() { return model::tallies.size(); }
//==============================================================================
// Fortran compatibility functions
//==============================================================================

243
src/tallies/tally_header.F90
View file

@ -97,22 +97,9 @@ module tally_header
character(len=104) :: name = "" ! user-defined name
real(8) :: volume ! volume of region
! Results for each bin -- the first dimension of the array is for scores
! (e.g. flux, total reaction rate, fission reaction rate, etc.) and the
! second dimension of the array is for the combination of filters
! (e.g. specific cell, specific energy group, etc.)
integer :: total_score_bins
real(C_DOUBLE), allocatable :: results(:,:,:)
! Number of realizations of tally random variables
integer :: n_realizations = 0
contains
procedure :: accumulate => tally_accumulate
procedure :: allocate_results => tally_allocate_results
procedure :: read_results_hdf5 => tally_read_results_hdf5
procedure :: write_results_hdf5 => tally_write_results_hdf5
procedure :: id => tally_get_id
procedure :: set_id => tally_set_id
procedure :: type => tally_get_type
@ -147,135 +134,8 @@ module tally_header
! in case the code needs to find an ID for a new tally.
integer :: largest_tally_id
! Global tallies
! 1) collision estimate of k-eff
! 2) absorption estimate of k-eff
! 3) track-length estimate of k-eff
! 4) leakage fraction
real(C_DOUBLE), public, allocatable, target :: global_tallies(:,:)
! Normalization for statistics
integer(C_INT32_T), public, bind(C) :: n_realizations = 0 ! # of independent realizations
real(C_DOUBLE), public, bind(C) :: total_weight ! total starting particle weight in realization
contains
!===============================================================================
! ACCUMULATE_TALLY
!===============================================================================
subroutine tally_accumulate(this)
class(TallyObject), intent(inout) :: this
integer :: i, j
real(C_DOUBLE) :: val
real(C_DOUBLE) :: total_source
interface
function total_source_strength() result(strength) bind(C)
import C_DOUBLE
real(C_DOUBLE) :: strength
end function
end interface
! Increment number of realizations
if (reduce_tallies) then
this % n_realizations = this % n_realizations + 1
else
this % n_realizations = this % n_realizations + n_procs
end if
if (master .or. (.not. reduce_tallies)) then
! Calculate total source strength for normalization
if (run_mode == MODE_FIXEDSOURCE) then
total_source = total_source_strength()
else
total_source = ONE
end if
! Accumulate each result
do j = 1, size(this % results, 3)
do i = 1, size(this % results, 2)
val = this % results(RESULT_VALUE, i, j)/total_weight * total_source
this % results(RESULT_VALUE, i, j) = ZERO
this % results(RESULT_SUM, i, j) = &
this % results(RESULT_SUM, i, j) + val
this % results(RESULT_SUM_SQ, i, j) = &
this % results(RESULT_SUM_SQ, i, j) + val*val
end do
end do
end if
end subroutine tally_accumulate
subroutine tally_write_results_hdf5(this, group_id)
class(TallyObject), intent(in) :: this
integer(HID_T), intent(in) :: group_id
integer(HSIZE_T) :: n_filter, n_score
interface
subroutine write_tally_results(group_id, n_filter, n_score, results) &
bind(C)
import HID_T, HSIZE_T, C_DOUBLE
integer(HID_T), value :: group_id
integer(HSIZE_T), value :: n_filter
integer(HSIZE_T), value :: n_score
real(C_DOUBLE), intent(in) :: results(*)
end subroutine write_tally_results
end interface
n_filter = size(this % results, 3)
n_score = size(this % results, 2)
call write_tally_results(group_id, n_filter, n_score, this % results)
end subroutine tally_write_results_hdf5
subroutine tally_read_results_hdf5(this, group_id)
class(TallyObject), intent(inout) :: this
integer(HID_T), intent(in) :: group_id
integer(HSIZE_T) :: n_filter, n_score
interface
subroutine read_tally_results(group_id, n_filter, n_score, results) &
bind(C)
import HID_T, HSIZE_T, C_DOUBLE
integer(HID_T), value :: group_id
integer(HSIZE_T), value :: n_filter
integer(HSIZE_T), value :: n_score
real(C_DOUBLE), intent(out) :: results(*)
end subroutine read_tally_results
end interface
n_filter = size(this % results, 3)
n_score = size(this % results, 2)
call read_tally_results(group_id, n_filter, n_score, this % results)
end subroutine tally_read_results_hdf5
!===============================================================================
! ALLOCATE_RESULTS allocates and initializes the results component of the
! TallyObject derived type
!===============================================================================
subroutine tally_allocate_results(this)
class(TallyObject), intent(inout) :: this
! Set total number of filter and scoring bins
this % total_score_bins = this % n_score_bins() * this % n_nuclide_bins()
if (allocated(this % results)) then
! If results was already allocated but shape is wrong, then reallocate it
! to the correct shape
if (this % total_score_bins /= size(this % results, 2) .or. &
this % n_filter_bins() /= size(this % results, 3)) then
deallocate(this % results)
allocate(this % results(3, this % total_score_bins, this % n_filter_bins()))
end if
else
allocate(this % results(3, this % total_score_bins, this % n_filter_bins()))
end if
end subroutine tally_allocate_results
function tally_get_id(this) result(t)
class(TallyObject) :: this
integer(C_INT) :: t
@ -502,28 +362,6 @@ contains
call tally_set_deriv_c(this % ptr, deriv)
end subroutine
!===============================================================================
! CONFIGURE_TALLIES initializes several data structures related to tallies. This
! is called after the basic tally data has already been read from the
! tallies.xml file.
!===============================================================================
subroutine allocate_tally_results() bind(C)
integer :: i
! Allocate global tallies
if (.not. allocated(global_tallies)) then
allocate(global_tallies(3, N_GLOBAL_TALLIES))
end if
! Allocate results arrays for tallies
do i = 1, n_tallies
call tallies(i) % obj % allocate_results()
end do
end subroutine allocate_tally_results
!===============================================================================
! FREE_MEMORY_TALLY deallocates global arrays defined in this module
!===============================================================================
@ -540,8 +378,6 @@ contains
if (allocated(tallies)) deallocate(tallies)
call tally_dict % clear()
largest_tally_id = 0
if (allocated(global_tallies)) deallocate(global_tallies)
end subroutine free_memory_tally
!===============================================================================
@ -612,20 +448,6 @@ contains
end function openmc_get_tally_index
function openmc_global_tallies(ptr) result(err) bind(C)
type(C_PTR), intent(out) :: ptr
integer(C_INT) :: err
if (.not. allocated(global_tallies)) then
err = E_ALLOCATE
call set_errmsg("Global tallies have not been allocated yet.")
else
err = 0
ptr = C_LOC(global_tallies)
end if
end function openmc_global_tallies
function openmc_tally_get_estimator(index, estimator) result(err) bind(C)
! Return the type of estimator of a tally
integer(C_INT32_T), value :: index
@ -658,71 +480,6 @@ contains
end function openmc_tally_get_id
function openmc_tally_get_n_realizations(index, n) result(err) bind(C)
! Return the number of realizations for a tally
integer(C_INT32_T), value :: index
integer(C_INT32_T), intent(out) :: n
integer(C_INT) :: err
if (index >= 1 .and. index <= size(tallies)) then
n = tallies(index) % obj % n_realizations
err = 0
else
err = E_OUT_OF_BOUNDS
call set_errmsg('Index in tallies array is out of bounds.')
end if
end function openmc_tally_get_n_realizations
function openmc_tally_reset(index) result(err) bind(C)
! Reset tally results and number of realizations
integer(C_INT32_T), intent(in), value :: index
integer(C_INT) :: err
if (index >= 1 .and. index <= size(tallies)) then
associate (t => tallies(index) % obj)
t % n_realizations = 0
if (allocated(t % results)) t % results(:, :, :) = ZERO
err = 0
end associate
else
err = E_OUT_OF_BOUNDS
call set_errmsg('Index in tallies array is out of bounds.')
end if
end function openmc_tally_reset
function openmc_tally_results(index, ptr, shape_) result(err) bind(C)
! Returns a pointer to a tally results array along with its shape. This
! allows a user to obtain in-memory tally results from Python directly.
integer(C_INT32_T), intent(in), value :: index
type(C_PTR), intent(out) :: ptr
integer(C_SIZE_T), intent(out) :: shape_(3)
integer(C_INT) :: err
if (index >= 1 .and. index <= size(tallies)) then
associate (t => tallies(index) % obj)
if (allocated(t % results)) then
ptr = C_LOC(t % results(1,1,1))
! Note that shape is reversed since it is assumed to be used from
! C/C++ code
shape_(1) = size(t % results, 3)
shape_(2) = size(t % results, 2)
shape_(3) = size(t % results, 1)
err = 0
else
err = E_ALLOCATE
call set_errmsg("Tally results have not been allocated yet.")
end if
end associate
else
err = E_OUT_OF_BOUNDS
call set_errmsg('Index in tallies array is out of bounds.')
end if
end function openmc_tally_results
function openmc_tally_set_estimator(index, estimator) result(err) bind(C)
! Set the type of estimator a tally
integer(C_INT32_T), value, intent(in) :: index

29
src/tallies/tally_scoring.cpp
View file

@ -146,7 +146,7 @@ void
score_fission_delayed_dg(int i_tally, int d_bin, double score, int score_index)
{
// Save the original delayed group bin
const Tally& tally {*model::tallies[i_tally]};
auto& tally {*model::tallies[i_tally]};
auto i_filt = tally.filters(tally.delayedgroup_filter_);
auto& dg_match {simulation::filter_matches[i_filt]};
auto i_bin = dg_match.i_bin_;
@ -164,10 +164,9 @@ score_fission_delayed_dg(int i_tally, int d_bin, double score, int score_index)
}
// Update the tally result
auto results = tally_results(i_tally);
//TODO: off-by-one
// TODO: off-by-one
#pragma omp atomic
results(filter_index-1, score_index, RESULT_VALUE) += score;
tally.results_(filter_index-1, score_index, RESULT_VALUE) += score;
// Reset the original delayed group bin
dg_match.bins_[i_bin] = original_bin;
@ -181,8 +180,7 @@ score_fission_delayed_dg(int i_tally, int d_bin, double score, int score_index)
void
score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin)
{
const Tally& tally {*model::tallies[i_tally]};
auto results = tally_results(i_tally);
auto& tally {*model::tallies[i_tally]};
auto i_eout_filt = tally.filters()[tally.energyout_filter_];
auto i_bin = simulation::filter_matches[i_eout_filt].i_bin_;
auto bin_energyout = simulation::filter_matches[i_eout_filt].bins_[i_bin];
@ -261,7 +259,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin)
// Update tally results
#pragma omp atomic
results(filter_index-1, i_score, RESULT_VALUE) += score;
tally.results_(filter_index-1, i_score, RESULT_VALUE) += score;
} else if (score_bin == SCORE_DELAYED_NU_FISSION && g != 0) {
@ -309,7 +307,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin)
// Update tally results
#pragma omp atomic
results(filter_index-1, i_score, RESULT_VALUE) += score*filter_weight;
tally.results_(filter_index-1, i_score, RESULT_VALUE) += score*filter_weight;
}
}
}
@ -328,8 +326,7 @@ void
score_general_ce(const Particle* p, int i_tally, int start_index,
int filter_index, int i_nuclide, double atom_density, double flux)
{
const Tally& tally {*model::tallies[i_tally]};
auto results = tally_results(i_tally);
auto& tally {*model::tallies[i_tally]};
// Get the pre-collision energy of the particle.
auto E = p->last_E;
@ -1222,7 +1219,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
// Update tally results
#pragma omp atomic
results(filter_index-1, score_index, RESULT_VALUE) += score;
tally.results_(filter_index-1, score_index, RESULT_VALUE) += score;
}
}
@ -1235,8 +1232,7 @@ void
score_general_mg(const Particle* p, int i_tally, int start_index,
int filter_index, int i_nuclide, double atom_density, double flux)
{
const Tally& tally {*model::tallies[i_tally]};
auto results = tally_results(i_tally);
auto& tally {*model::tallies[i_tally]};
//TODO: off-by-one throughout on p->material
@ -1930,7 +1926,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
// Update tally results
#pragma omp atomic
results(filter_index-1, score_index, RESULT_VALUE) += score;
tally.results_(filter_index-1, score_index, RESULT_VALUE) += score;
}
}
@ -2223,8 +2219,7 @@ score_surface_tally(const Particle* p, const std::vector<int>& tallies)
double flux = p->wgt;
for (auto i_tally : tallies) {
const Tally& tally {*model::tallies[i_tally]};
auto results = tally_results(i_tally);
auto& tally {*model::tallies[i_tally]};
// Initialize an iterator over valid filter bin combinations. If there are
// no valid combinations, use a continue statement to ensure we skip the
@ -2246,7 +2241,7 @@ score_surface_tally(const Particle* p, const std::vector<int>& tallies)
++score_index) {
//TODO: off-by-one
#pragma omp atomic
results(filter_index-1, score_index, RESULT_VALUE) += score;
tally.results_(filter_index-1, score_index, RESULT_VALUE) += score;
}
}

17
src/tallies/trigger.cpp
View file

@ -29,14 +29,12 @@ namespace settings {
std::pair<double, double>
get_tally_uncertainty(int i_tally, int score_index, int filter_index)
{
int n;
//TODO: off-by-one
int err = openmc_tally_get_n_realizations(i_tally+1, &n);
const auto& tally {model::tallies[i_tally]};
auto results = tally_results(i_tally);
auto sum = results(filter_index, score_index, RESULT_SUM);
auto sum_sq = results(filter_index, score_index, RESULT_SUM_SQ);
auto sum = tally->results_(filter_index, score_index, RESULT_SUM);
auto sum_sq = tally->results_(filter_index, score_index, RESULT_SUM_SQ);
int n = tally->n_realizations_;
auto mean = sum / n;
double std_dev = std::sqrt((sum_sq/n - mean*mean) / (n-1));
double rel_err = (mean != 0.) ? std_dev / std::abs(mean) : 0.;
@ -59,13 +57,10 @@ check_tally_triggers(double& ratio, int& tally_id, int& score)
const Tally& t {*model::tallies[i_tally]};
// Ignore tallies with less than two realizations.
int n_reals;
//TODO: off-by-one
int err = openmc_tally_get_n_realizations(i_tally+1, &n_reals);
if (n_reals < 2) continue;
if (t.n_realizations_ < 2) continue;
for (const auto& trigger : t.triggers_) {
const auto& results = tally_results(i_tally);
const auto& results = t.results_;
for (auto filter_index = 0; filter_index < results.shape()[0];
++filter_index) {
for (auto score_index = 0; score_index < results.shape()[1];