diff --git a/include/openmc/constants.h b/include/openmc/constants.h index fa2251b84d..8bb2cb1bfa 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -5,6 +5,7 @@ #define OPENMC_CONSTANTS_H #include +#include #include #include "openmc/array.h" @@ -339,6 +340,11 @@ enum class RunMode { enum class GeometryType { CSG, DAG }; +//============================================================================== +// Volume Calculation Constants + +constexpr uint64_t UINT64_T_MAX {std::numeric_limits::max()}; + } // namespace openmc #endif // OPENMC_CONSTANTS_H diff --git a/include/openmc/volume_calc.h b/include/openmc/volume_calc.h index db96f250f9..2efd955f04 100644 --- a/include/openmc/volume_calc.h +++ b/include/openmc/volume_calc.h @@ -1,6 +1,9 @@ #ifndef OPENMC_VOLUME_CALC_H #define OPENMC_VOLUME_CALC_H +#include +#include + #include "openmc/array.h" #include "openmc/position.h" #include "openmc/tallies/trigger.h" @@ -10,7 +13,6 @@ #include "xtensor/xtensor.hpp" #include -#include namespace openmc { @@ -69,7 +71,8 @@ private: //! \param[in] i_material Index in global materials vector //! \param[in,out] indices Vector of material indices //! \param[in,out] hits Number of hits corresponding to each material - void check_hit(int i_material, vector& indices, vector& hits) const; + void check_hit( + int i_material, vector& indices, vector& hits) const; }; //============================================================================== diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index 3ae98fdc38..e36b889041 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -99,16 +99,16 @@ vector VolumeCalculation::execute() const { // Shared data that is collected from all threads int n = domain_ids_.size(); - vector> master_indices( + vector> master_indices( n); // List of material indices for each domain - vector> master_hits( + vector> master_hits( n); // Number of hits for each material in each domain int iterations = 0; // Divide work over MPI processes - size_t min_samples = n_samples_ / mpi::n_procs; - size_t remainder = n_samples_ % mpi::n_procs; - size_t i_start, i_end; + uint64_t min_samples = n_samples_ / mpi::n_procs; + uint64_t remainder = n_samples_ % mpi::n_procs; + uint64_t i_start, i_end; if (mpi::rank < remainder) { i_start = (min_samples + 1) * mpi::rank; i_end = i_start + min_samples + 1; @@ -123,14 +123,14 @@ vector VolumeCalculation::execute() const #pragma omp parallel { // Variables that are private to each thread - vector> indices(n); - vector> hits(n); + vector> indices(n); + vector> hits(n); Particle p; // Sample locations and count hits #pragma omp for for (size_t i = i_start; i < i_end; i++) { - int64_t id = iterations * n_samples_ + i; + uint64_t id = iterations * n_samples_ + i; uint64_t seed = init_seed(id, STREAM_VOLUME); p.n_coord() = 1; @@ -223,7 +223,15 @@ vector VolumeCalculation::execute() const // bump iteration counter and get total number // of samples at this point iterations++; - size_t total_samples = iterations * n_samples_; + uint64_t total_samples = iterations * n_samples_; + + // warn user if total sample size is greater than what the size_t type can + // represent + if (total_samples == UINT64_T_MAX) { + warning("The number of samples has exceeded the type used to track hits. " + "Volume " + "results may be inaccurate."); + } // reset double trigger_val = -INFTY; @@ -246,10 +254,11 @@ vector VolumeCalculation::execute() const if (mpi::master) { for (int j = 1; j < mpi::n_procs; j++) { int q; + // retrieve results MPI_Recv( - &q, 1, MPI_INTEGER, j, 2 * j, mpi::intracomm, MPI_STATUS_IGNORE); - vector buffer(2 * q); - MPI_Recv(buffer.data(), 2 * q, MPI_INTEGER, j, 2 * j + 1, + &q, 1, MPI_UINT64_T, j, 2 * j, mpi::intracomm, MPI_STATUS_IGNORE); + vector buffer(2 * q); + MPI_Recv(buffer.data(), 2 * q, MPI_UINT64_T, j, 2 * j + 1, mpi::intracomm, MPI_STATUS_IGNORE); for (int k = 0; k < q; ++k) { bool already_added = false; @@ -268,20 +277,20 @@ vector VolumeCalculation::execute() const } } else { int q = master_indices[i_domain].size(); - vector buffer(2 * q); + vector buffer(2 * q); for (int k = 0; k < q; ++k) { buffer[2 * k] = master_indices[i_domain][k]; buffer[2 * k + 1] = master_hits[i_domain][k]; } - MPI_Send(&q, 1, MPI_INTEGER, 0, 2 * mpi::rank, mpi::intracomm); - MPI_Send(buffer.data(), 2 * q, MPI_INTEGER, 0, 2 * mpi::rank + 1, + MPI_Send(&q, 1, MPI_UINT64_T, 0, 2 * mpi::rank, mpi::intracomm); + MPI_Send(buffer.data(), 2 * q, MPI_UINT64_T, 0, 2 * mpi::rank + 1, mpi::intracomm); } #endif if (mpi::master) { - int total_hits = 0; + size_t total_hits = 0; for (int j = 0; j < master_indices[i_domain].size(); ++j) { total_hits += master_hits[i_domain][j]; double f = @@ -464,7 +473,7 @@ void VolumeCalculation::to_hdf5( } void VolumeCalculation::check_hit( - int i_material, vector& indices, vector& hits) const + int i_material, vector& indices, vector& hits) const { // Check if this material was previously hit and if so, increment count