diff --git a/CMakeLists.txt b/CMakeLists.txt index b8cda33ba..bbeb774ad 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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 diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 5ec852d2d..ef84b78f4 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -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; diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index c7053c608..7b30bb792 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -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 k_generation; diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index a4f68bcbe..75fa3d10a 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -5,6 +5,7 @@ #include "openmc/tallies/trigger.h" #include "pugixml.hpp" +#include "xtensor/xfixed.hpp" #include "xtensor/xtensor.hpp" #include // 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 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 results_; + //---------------------------------------------------------------------------- // Miscellaneous public members. @@ -102,17 +116,22 @@ private: //============================================================================== namespace model { + extern std::vector> tallies; + extern std::vector active_tallies; + extern std::vector active_analog_tallies; + extern std::vector active_tracklength_tallies; + extern std::vector active_collision_tallies; + extern std::vector active_meshsurf_tallies; + extern std::vector active_surface_tallies; +} -extern std::vector> tallies; +namespace simulation { + //! Global tallies (such as k-effective estimators) + extern xt::xtensor_fixed> global_tallies; -extern std::vector active_tallies; -extern std::vector active_analog_tallies; -extern std::vector active_tracklength_tallies; -extern std::vector active_collision_tallies; -extern std::vector active_meshsurf_tallies; -extern std::vector 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 using adaptor_type = xt::xtensor_adaptor, 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(); diff --git a/src/constants.F90 b/src/constants.F90 index 1fd5d00e9..ff74cd68a 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -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, & diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index 0f52649d8..b3678a021 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -46,7 +46,7 @@ xt::xtensor 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 kv {}; xt::xtensor cov = xt::zeros({3, 3}); diff --git a/src/finalize.cpp b/src/finalize.cpp index e67708ee1..c38577716 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -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; diff --git a/src/output.cpp b/src/output.cpp index 400e56d32..29d873173 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -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"; diff --git a/src/photon.cpp b/src/photon.cpp index d3fa9f465..7a6bc2452 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -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 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); } diff --git a/src/simulation.cpp b/src/simulation.cpp index 89bd93421..35d8d14a6 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -21,6 +21,7 @@ #include "openmc/tallies/trigger.h" #include +#include "xtensor/xview.hpp" #include #include @@ -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 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 diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index 8976c7047..675306346 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -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 diff --git a/src/source.cpp b/src/source.cpp index 231a95e21..187b2e6b4 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -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()) { diff --git a/src/state_point.F90 b/src/state_point.F90 deleted file mode 100644 index 33daef4c3..000000000 --- a/src/state_point.F90 +++ /dev/null @@ -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 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 diff --git a/src/state_point.cpp b/src/state_point.cpp index f86a7b131..5e409322f 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -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 results_copy = xt::zeros_like(results); + xt::xtensor 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; diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 1304b02d7..ba4ee7f2f 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -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 !=============================================================================== diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 93ba72a19..31806b5b4 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -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> tallies; - std::vector active_tallies; std::vector active_analog_tallies; std::vector active_tracklength_tallies; @@ -47,6 +47,11 @@ namespace model { std::vector active_surface_tallies; } +namespace simulation { + xt::xtensor_fixed> 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({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 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 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 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 //============================================================================== diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index f61c99a98..0167a23fb 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -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 diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 690f61f81..9edac5ec2 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -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& 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& 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; } } diff --git a/src/tallies/trigger.cpp b/src/tallies/trigger.cpp index a53cef817..2f7f36920 100644 --- a/src/tallies/trigger.cpp +++ b/src/tallies/trigger.cpp @@ -29,14 +29,12 @@ namespace settings { std::pair 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];