From 04393200cc79defe5cbfc5484376ea5d753f067c Mon Sep 17 00:00:00 2001 From: John Tramm Date: Tue, 11 Feb 2025 11:23:26 -0600 Subject: [PATCH] Random Ray Source Region Refactor (#3288) --- CMakeLists.txt | 1 + include/openmc/openmp_interface.h | 12 +- .../openmc/random_ray/flat_source_domain.h | 132 +----- .../openmc/random_ray/linear_source_domain.h | 30 +- include/openmc/random_ray/source_region.h | 400 ++++++++++++++++++ src/random_ray/flat_source_domain.cpp | 290 +++++-------- src/random_ray/linear_source_domain.cpp | 129 ++---- src/random_ray/random_ray.cpp | 60 ++- src/random_ray/random_ray_simulation.cpp | 2 +- src/random_ray/source_region.cpp | 236 +++++++++++ 10 files changed, 832 insertions(+), 460 deletions(-) create mode 100644 include/openmc/random_ray/source_region.h create mode 100644 src/random_ray/source_region.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index f7079aebc..94caec757 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -401,6 +401,7 @@ list(APPEND libopenmc_SOURCES src/random_ray/flat_source_domain.cpp src/random_ray/linear_source_domain.cpp src/random_ray/moment_matrix.cpp + src/random_ray/source_region.cpp src/reaction.cpp src/reaction_product.cpp src/scattdata.cpp diff --git a/include/openmc/openmp_interface.h b/include/openmc/openmp_interface.h index fe517415b..30bd6f008 100644 --- a/include/openmc/openmp_interface.h +++ b/include/openmc/openmp_interface.h @@ -36,13 +36,15 @@ inline int thread_num() class OpenMPMutex { public: - OpenMPMutex() + void init() { #ifdef _OPENMP omp_init_lock(&mutex_); #endif } + OpenMPMutex() { init(); } + ~OpenMPMutex() { #ifdef _OPENMP @@ -62,14 +64,10 @@ public: // rather, it produces two different mutexes. // Copy constructor - OpenMPMutex(const OpenMPMutex& other) { OpenMPMutex(); } + OpenMPMutex(const OpenMPMutex& other) { init(); } // Copy assignment operator - OpenMPMutex& operator=(const OpenMPMutex& other) - { - OpenMPMutex(); - return *this; - } + OpenMPMutex& operator=(const OpenMPMutex& other) { return *this; } //! Lock the mutex. // diff --git a/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h index b5b5db706..011ff7ccd 100644 --- a/include/openmc/random_ray/flat_source_domain.h +++ b/include/openmc/random_ray/flat_source_domain.h @@ -4,83 +4,12 @@ #include "openmc/constants.h" #include "openmc/openmp_interface.h" #include "openmc/position.h" +#include "openmc/random_ray/source_region.h" #include "openmc/source.h" +#include namespace openmc { -//---------------------------------------------------------------------------- -// Helper Functions - -// The hash_combine function is the standard hash combine function from boost -// that is typically used for combining multiple hash values into a single hash -// as is needed for larger objects being stored in a hash map. The function is -// taken from: -// https://www.boost.org/doc/libs/1_55_0/doc/html/hash/reference.html#boost.hash_combine -// which carries the following license: -// -// Boost Software License - Version 1.0 - August 17th, 2003 -// Permission is hereby granted, free of charge, to any person or organization -// obtaining a copy of the software and accompanying documentation covered by -// this license (the "Software") to use, reproduce, display, distribute, -// execute, and transmit the Software, and to prepare derivative works of the -// Software, and to permit third-parties to whom the Software is furnished to -// do so, all subject to the following: -// The copyright notices in the Software and this entire statement, including -// the above license grant, this restriction and the following disclaimer, -// must be included in all copies of the Software, in whole or in part, and -// all derivative works of the Software, unless such copies or derivative -// works are solely in the form of machine-executable object code generated by -// a source language processor. -// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -// FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT -// SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE -// FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, -// ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER -// DEALINGS IN THE SOFTWARE. -inline void hash_combine(size_t& seed, const size_t v) -{ - seed ^= (v + 0x9e3779b9 + (seed << 6) + (seed >> 2)); -} - -//---------------------------------------------------------------------------- -// Helper Structs - -// A mapping object that is used to map between a specific random ray -// source region and an OpenMC native tally bin that it should score to -// every iteration. -struct TallyTask { - int tally_idx; - int filter_idx; - int score_idx; - int score_type; - TallyTask(int tally_idx, int filter_idx, int score_idx, int score_type) - : tally_idx(tally_idx), filter_idx(filter_idx), score_idx(score_idx), - score_type(score_type) - {} - TallyTask() = default; - - // Comparison and Hash operators are defined to allow usage of the - // TallyTask struct as a key in an unordered_set - bool operator==(const TallyTask& other) const - { - return tally_idx == other.tally_idx && filter_idx == other.filter_idx && - score_idx == other.score_idx && score_type == other.score_type; - } - - struct HashFunctor { - size_t operator()(const TallyTask& task) const - { - size_t seed = 0; - hash_combine(seed, task.tally_idx); - hash_combine(seed, task.filter_idx); - hash_combine(seed, task.score_idx); - hash_combine(seed, task.score_type); - return seed; - } - }; -}; - /* * The FlatSourceDomain class encompasses data and methods for storing * scalar flux and source region for all flat source regions in a @@ -108,15 +37,16 @@ public: void random_ray_tally(); virtual void accumulate_iteration_flux(); void output_to_vtk() const; - virtual void all_reduce_replicated_source_regions(); + void all_reduce_replicated_source_regions(); void convert_external_sources(); void count_external_source_regions(); void set_adjoint_sources(const vector& forward_flux); - virtual void flux_swap(); + void flux_swap(); virtual double evaluate_flux_at_point(Position r, int64_t sr, int g) const; double compute_fixed_source_normalization_factor() const; void flatten_xs(); void transpose_scattering_matrix(); + void serialize_final_fluxes(vector& flux); //---------------------------------------------------------------------------- // Static Data members @@ -129,7 +59,6 @@ public: //---------------------------------------------------------------------------- // Public Data members - bool mapped_all_tallies_ {false}; // If all source regions have been visited int64_t n_source_regions_ {0}; // Total number of source regions in the model @@ -140,22 +69,6 @@ public: // in model::cells vector source_region_offsets_; - // 1D arrays representing values for all source regions - vector lock_; - vector volume_; - vector volume_t_; - vector position_recorded_; - vector position_; - - // 2D arrays stored in 1D representing values for all source regions x energy - // groups - vector scalar_flux_old_; - vector scalar_flux_new_; - vector source_; - vector external_source_; - vector external_source_present_; - vector scalar_flux_final_; - // 2D arrays stored in 1D representing values for all materials x energy // groups int n_materials_; @@ -168,20 +81,22 @@ public: // groups x energy groups vector sigma_s_; + // The abstract container holding all source region-specific data + SourceRegionContainer source_regions_; + protected: //---------------------------------------------------------------------------- // Methods void apply_external_source_to_source_region( - Discrete* discrete, double strength_factor, int64_t source_region); + Discrete* discrete, double strength_factor, int64_t sr); void apply_external_source_to_cell_instances(int32_t i_cell, Discrete* discrete, double strength_factor, int target_material_id, const vector& instances); void apply_external_source_to_cell_and_children(int32_t i_cell, Discrete* discrete, double strength_factor, int32_t target_material_id); - virtual void set_flux_to_flux_plus_source( - int64_t idx, double volume, int material, int g); - void set_flux_to_source(int64_t idx); - virtual void set_flux_to_old_flux(int64_t idx); + virtual void set_flux_to_flux_plus_source(int64_t sr, double volume, int g); + void set_flux_to_source(int64_t sr, int g); + virtual void set_flux_to_old_flux(int64_t sr, int g); //---------------------------------------------------------------------------- // Private data members @@ -193,20 +108,6 @@ protected: simulation_volume_; // Total physical volume of the simulation domain, as // defined by the 3D box of the random ray source - // 2D array representing values for all source elements x tally - // tasks - vector> tally_task_; - - // 1D array representing values for all source regions, with each region - // containing a set of volume tally tasks. This more complicated data - // structure is convenient for ensuring that volumes are only tallied once per - // source region, regardless of how many energy groups are used for tallying. - vector> volume_task_; - - // 1D arrays representing values for all source regions - vector material_; - vector volume_naive_; - // Volumes for each tally and bin/score combination. This intermediate data // structure is used when tallying quantities that must be normalized by // volume (i.e., flux). The vector is index by tally index, while the inner 2D @@ -250,15 +151,6 @@ T convert_to_big_endian(T in) return out; } -template -void parallel_fill(vector& arr, T value) -{ -#pragma omp parallel for schedule(static) - for (int i = 0; i < arr.size(); i++) { - arr[i] = value; - } -} - } // namespace openmc #endif // OPENMC_RANDOM_RAY_FLAT_SOURCE_DOMAIN_H diff --git a/include/openmc/random_ray/linear_source_domain.h b/include/openmc/random_ray/linear_source_domain.h index 4812d1433..67fdd99f8 100644 --- a/include/openmc/random_ray/linear_source_domain.h +++ b/include/openmc/random_ray/linear_source_domain.h @@ -18,48 +18,22 @@ namespace openmc { class LinearSourceDomain : public FlatSourceDomain { public: - //---------------------------------------------------------------------------- - // Constructors - LinearSourceDomain(); - //---------------------------------------------------------------------------- // Methods void update_neutron_source(double k_eff) override; - double compute_k_eff(double k_eff_old) const; void normalize_scalar_flux_and_volumes( double total_active_distance_per_iteration) override; void batch_reset() override; - void convert_source_regions_to_tallies(); - void reset_tally_volumes(); - void random_ray_tally(); void accumulate_iteration_flux() override; void output_to_vtk() const; - void all_reduce_replicated_source_regions() override; - void convert_external_sources(); - void count_external_source_regions(); - void flux_swap() override; double evaluate_flux_at_point(Position r, int64_t sr, int g) const override; - //---------------------------------------------------------------------------- - // Public Data members - - vector source_gradients_; - vector flux_moments_old_; - vector flux_moments_new_; - vector flux_moments_t_; - vector centroid_; - vector centroid_iteration_; - vector centroid_t_; - vector mom_matrix_; - vector mom_matrix_t_; - protected: //---------------------------------------------------------------------------- // Methods - void set_flux_to_flux_plus_source( - int64_t idx, double volume, int material, int g) override; - void set_flux_to_old_flux(int64_t idx) override; + void set_flux_to_flux_plus_source(int64_t sr, double volume, int g) override; + void set_flux_to_old_flux(int64_t sr, int g) override; }; // class LinearSourceDomain diff --git a/include/openmc/random_ray/source_region.h b/include/openmc/random_ray/source_region.h new file mode 100644 index 000000000..0dc617a64 --- /dev/null +++ b/include/openmc/random_ray/source_region.h @@ -0,0 +1,400 @@ +#ifndef OPENMC_RANDOM_RAY_SOURCE_REGION_H +#define OPENMC_RANDOM_RAY_SOURCE_REGION_H + +#include "openmc/openmp_interface.h" +#include "openmc/position.h" +#include "openmc/random_ray/moment_matrix.h" +#include "openmc/settings.h" + +namespace openmc { + +//---------------------------------------------------------------------------- +// Helper Functions + +// The hash_combine function is the standard hash combine function from boost +// that is typically used for combining multiple hash values into a single hash +// as is needed for larger objects being stored in a hash map. The function is +// taken from: +// https://www.boost.org/doc/libs/1_55_0/doc/html/hash/reference.html#boost.hash_combine +// which carries the following license: +// +// Boost Software License - Version 1.0 - August 17th, 2003 +// Permission is hereby granted, free of charge, to any person or organization +// obtaining a copy of the software and accompanying documentation covered by +// this license (the "Software") to use, reproduce, display, distribute, +// execute, and transmit the Software, and to prepare derivative works of the +// Software, and to permit third-parties to whom the Software is furnished to +// do so, all subject to the following: +// The copyright notices in the Software and this entire statement, including +// the above license grant, this restriction and the following disclaimer, +// must be included in all copies of the Software, in whole or in part, and +// all derivative works of the Software, unless such copies or derivative +// works are solely in the form of machine-executable object code generated by +// a source language processor. +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT +// SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE +// FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, +// ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +// DEALINGS IN THE SOFTWARE. +inline void hash_combine(size_t& seed, const size_t v) +{ + seed ^= (v + 0x9e3779b9 + (seed << 6) + (seed >> 2)); +} + +//---------------------------------------------------------------------------- +// Helper Structs + +// A mapping object that is used to map between a specific random ray +// source region and an OpenMC native tally bin that it should score to +// every iteration. +struct TallyTask { + int tally_idx; + int filter_idx; + int score_idx; + int score_type; + TallyTask(int tally_idx, int filter_idx, int score_idx, int score_type) + : tally_idx(tally_idx), filter_idx(filter_idx), score_idx(score_idx), + score_type(score_type) + {} + TallyTask() = default; + + // Comparison and Hash operators are defined to allow usage of the + // TallyTask struct as a key in an unordered_set + bool operator==(const TallyTask& other) const + { + return tally_idx == other.tally_idx && filter_idx == other.filter_idx && + score_idx == other.score_idx && score_type == other.score_type; + } + + struct HashFunctor { + size_t operator()(const TallyTask& task) const + { + size_t seed = 0; + hash_combine(seed, task.tally_idx); + hash_combine(seed, task.filter_idx); + hash_combine(seed, task.score_idx); + hash_combine(seed, task.score_type); + return seed; + } + }; +}; + +class SourceRegion { +public: + //---------------------------------------------------------------------------- + // Constructors + SourceRegion(int negroups, bool is_linear); + SourceRegion() = default; + + //---------------------------------------------------------------------------- + // Public Data members + + // Scalar fields + int material_ {0}; + OpenMPMutex lock_; + double volume_ {0.0}; + double volume_t_ {0.0}; + double volume_naive_ {0.0}; + int position_recorded_ {0}; + int external_source_present_ {0}; + Position position_ {0.0, 0.0, 0.0}; + Position centroid_ {0.0, 0.0, 0.0}; + Position centroid_iteration_ {0.0, 0.0, 0.0}; + Position centroid_t_ {0.0, 0.0, 0.0}; + MomentMatrix mom_matrix_ {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + MomentMatrix mom_matrix_t_ {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + // A set of volume tally tasks. This more complicated data structure is + // convenient for ensuring that volumes are only tallied once per source + // region, regardless of how many energy groups are used for tallying. + std::unordered_set volume_task_; + + // Energy group-wise 1D arrays + vector scalar_flux_old_; + vector scalar_flux_new_; + vector source_; + vector external_source_; + vector scalar_flux_final_; + + vector source_gradients_; + vector flux_moments_old_; + vector flux_moments_new_; + vector flux_moments_t_; + + // 2D array representing values for all energy groups x tally + // tasks. Each group may have a different number of tally tasks + // associated with it, necessitating the use of a jagged array. + vector> tally_task_; + +}; // class SourceRegion + +class SourceRegionContainer { +public: + //---------------------------------------------------------------------------- + // Constructors + SourceRegionContainer(int negroups, bool is_linear) + : negroups_(negroups), is_linear_(is_linear) + {} + SourceRegionContainer() = default; + + //---------------------------------------------------------------------------- + // Public Accessors + int& material(int64_t sr) { return material_[sr]; } + const int& material(int64_t sr) const { return material_[sr]; } + + OpenMPMutex& lock(int64_t sr) { return lock_[sr]; } + const OpenMPMutex& lock(int64_t sr) const { return lock_[sr]; } + + double& volume(int64_t sr) { return volume_[sr]; } + const double& volume(int64_t sr) const { return volume_[sr]; } + + double& volume_t(int64_t sr) { return volume_t_[sr]; } + const double& volume_t(int64_t sr) const { return volume_t_[sr]; } + + double& volume_naive(int64_t sr) { return volume_naive_[sr]; } + const double& volume_naive(int64_t sr) const { return volume_naive_[sr]; } + + int& position_recorded(int64_t sr) { return position_recorded_[sr]; } + const int& position_recorded(int64_t sr) const + { + return position_recorded_[sr]; + } + + int& external_source_present(int64_t sr) + { + return external_source_present_[sr]; + } + const int& external_source_present(int64_t sr) const + { + return external_source_present_[sr]; + } + + Position& position(int64_t sr) { return position_[sr]; } + const Position& position(int64_t sr) const { return position_[sr]; } + + Position& centroid(int64_t sr) { return centroid_[sr]; } + const Position& centroid(int64_t sr) const { return centroid_[sr]; } + + Position& centroid_iteration(int64_t sr) { return centroid_iteration_[sr]; } + const Position& centroid_iteration(int64_t sr) const + { + return centroid_iteration_[sr]; + } + + Position& centroid_t(int64_t sr) { return centroid_t_[sr]; } + const Position& centroid_t(int64_t sr) const { return centroid_t_[sr]; } + + MomentMatrix& mom_matrix(int64_t sr) { return mom_matrix_[sr]; } + const MomentMatrix& mom_matrix(int64_t sr) const { return mom_matrix_[sr]; } + + MomentMatrix& mom_matrix_t(int64_t sr) { return mom_matrix_t_[sr]; } + const MomentMatrix& mom_matrix_t(int64_t sr) const + { + return mom_matrix_t_[sr]; + } + + MomentArray& source_gradients(int64_t sr, int g) + { + return source_gradients_[index(sr, g)]; + } + const MomentArray& source_gradients(int64_t sr, int g) const + { + return source_gradients_[index(sr, g)]; + } + MomentArray& source_gradients(int64_t se) { return source_gradients_[se]; } + const MomentArray& source_gradients(int64_t se) const + { + return source_gradients_[se]; + } + + MomentArray& flux_moments_old(int64_t sr, int g) + { + return flux_moments_old_[index(sr, g)]; + } + const MomentArray& flux_moments_old(int64_t sr, int g) const + { + return flux_moments_old_[index(sr, g)]; + } + MomentArray& flux_moments_old(int64_t se) { return flux_moments_old_[se]; } + const MomentArray& flux_moments_old(int64_t se) const + { + return flux_moments_old_[se]; + } + + MomentArray& flux_moments_new(int64_t sr, int g) + { + return flux_moments_new_[index(sr, g)]; + } + const MomentArray& flux_moments_new(int64_t sr, int g) const + { + return flux_moments_new_[index(sr, g)]; + } + MomentArray& flux_moments_new(int64_t se) { return flux_moments_new_[se]; } + const MomentArray& flux_moments_new(int64_t se) const + { + return flux_moments_new_[se]; + } + + MomentArray& flux_moments_t(int64_t sr, int g) + { + return flux_moments_t_[index(sr, g)]; + } + const MomentArray& flux_moments_t(int64_t sr, int g) const + { + return flux_moments_t_[index(sr, g)]; + } + MomentArray& flux_moments_t(int64_t se) { return flux_moments_t_[se]; } + const MomentArray& flux_moments_t(int64_t se) const + { + return flux_moments_t_[se]; + } + + double& scalar_flux_old(int64_t sr, int g) + { + return scalar_flux_old_[index(sr, g)]; + } + const double& scalar_flux_old(int64_t sr, int g) const + { + return scalar_flux_old_[index(sr, g)]; + } + double& scalar_flux_old(int64_t se) { return scalar_flux_old_[se]; } + const double& scalar_flux_old(int64_t se) const + { + return scalar_flux_old_[se]; + } + + double& scalar_flux_new(int64_t sr, int g) + { + return scalar_flux_new_[index(sr, g)]; + } + const double& scalar_flux_new(int64_t sr, int g) const + { + return scalar_flux_new_[index(sr, g)]; + } + double& scalar_flux_new(int64_t se) { return scalar_flux_new_[se]; } + const double& scalar_flux_new(int64_t se) const + { + return scalar_flux_new_[se]; + } + + double& scalar_flux_final(int64_t sr, int g) + { + return scalar_flux_final_[index(sr, g)]; + } + const double& scalar_flux_final(int64_t sr, int g) const + { + return scalar_flux_final_[index(sr, g)]; + } + double& scalar_flux_final(int64_t se) { return scalar_flux_final_[se]; } + const double& scalar_flux_final(int64_t se) const + { + return scalar_flux_final_[se]; + } + + float& source(int64_t sr, int g) { return source_[index(sr, g)]; } + const float& source(int64_t sr, int g) const { return source_[index(sr, g)]; } + float& source(int64_t se) { return source_[se]; } + const float& source(int64_t se) const { return source_[se]; } + + float& external_source(int64_t sr, int g) + { + return external_source_[index(sr, g)]; + } + const float& external_source(int64_t sr, int g) const + { + return external_source_[index(sr, g)]; + } + float& external_source(int64_t se) { return external_source_[se]; } + const float& external_source(int64_t se) const + { + return external_source_[se]; + } + + vector& tally_task(int64_t sr, int g) + { + return tally_task_[index(sr, g)]; + } + const vector& tally_task(int64_t sr, int g) const + { + return tally_task_[index(sr, g)]; + } + vector& tally_task(int64_t se) { return tally_task_[se]; } + const vector& tally_task(int64_t se) const + { + return tally_task_[se]; + } + + std::unordered_set& volume_task(int64_t sr) + { + return volume_task_[sr]; + } + const std::unordered_set& volume_task( + int64_t sr) const + { + return volume_task_[sr]; + } + + //---------------------------------------------------------------------------- + // Public Methods + + void push_back(const SourceRegion& sr); + void assign(int n_source_regions, const SourceRegion& source_region); + void flux_swap(); + void mpi_sync_ranks(bool reduce_position); + +private: + //---------------------------------------------------------------------------- + // Private Data Members + int64_t n_source_regions_ {0}; + int negroups_ {0}; + bool is_linear_ {false}; + + // SoA storage for scalar fields (one item per source region) + vector material_; + vector lock_; + vector volume_; + vector volume_t_; + vector volume_naive_; + vector position_recorded_; + vector external_source_present_; + vector position_; + vector centroid_; + vector centroid_iteration_; + vector centroid_t_; + vector mom_matrix_; + vector mom_matrix_t_; + // A set of volume tally tasks. This more complicated data structure is + // convenient for ensuring that volumes are only tallied once per source + // region, regardless of how many energy groups are used for tallying. + vector> volume_task_; + + // SoA energy group-wise 2D arrays flattened to 1D + vector scalar_flux_old_; + vector scalar_flux_new_; + vector scalar_flux_final_; + vector source_; + vector external_source_; + + vector source_gradients_; + vector flux_moments_old_; + vector flux_moments_new_; + vector flux_moments_t_; + + // SoA 3D array representing values for all source regions x energy groups x + // tally tasks. The outer two dimensions (source regions and energy groups) + // are flattened to 1D. Each group may have a different number of tally tasks + // associated with it, necessitating the use of a jagged array for the inner + // dimension. + vector> tally_task_; + + //---------------------------------------------------------------------------- + // Private Methods + + // Helper function for indexing + inline int index(int64_t sr, int g) const { return sr * negroups_ + g; } +}; + +} // namespace openmc + +#endif // OPENMC_RANDOM_RAY_SOURCE_REGION_H \ No newline at end of file diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 19913c03c..39a50ace6 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -46,40 +46,17 @@ FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_) } // Initialize cell-wise arrays - lock_.resize(n_source_regions_); - material_.resize(n_source_regions_); - position_recorded_.assign(n_source_regions_, 0); - position_.resize(n_source_regions_); - volume_.assign(n_source_regions_, 0.0); - volume_t_.assign(n_source_regions_, 0.0); - volume_naive_.assign(n_source_regions_, 0.0); + bool is_linear = RandomRay::source_shape_ != RandomRaySourceShape::FLAT; + source_regions_ = SourceRegionContainer(negroups_, is_linear); + source_regions_.assign(n_source_regions_, SourceRegion(negroups_, is_linear)); - // Initialize element-wise arrays - scalar_flux_new_.assign(n_source_elements_, 0.0); - scalar_flux_final_.assign(n_source_elements_, 0.0); - source_.resize(n_source_elements_); - - tally_task_.resize(n_source_elements_); - volume_task_.resize(n_source_regions_); - - if (settings::run_mode == RunMode::EIGENVALUE) { - // If in eigenvalue mode, set starting flux to guess of unity - scalar_flux_old_.assign(n_source_elements_, 1.0); - } else { - // If in fixed source mode, set starting flux to guess of zero - // and initialize external source arrays - scalar_flux_old_.assign(n_source_elements_, 0.0); - external_source_.assign(n_source_elements_, 0.0); - external_source_present_.assign(n_source_regions_, false); - } - - // Initialize material array + // Initialize materials int64_t source_region_id = 0; for (int i = 0; i < model::cells.size(); i++) { Cell& cell = *model::cells[i]; if (cell.type_ == Fill::MATERIAL) { for (int j = 0; j < cell.n_instances_; j++) { - material_[source_region_id++] = cell.material(j); + source_regions_.material(source_region_id++) = cell.material(j); } } } @@ -113,17 +90,23 @@ FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_) void FlatSourceDomain::batch_reset() { - // Reset scalar fluxes, iteration volume tallies, and region hit flags to - // zero - parallel_fill(scalar_flux_new_, 0.0); - parallel_fill(volume_, 0.0); +// Reset scalar fluxes and iteration volume tallies to zero +#pragma omp parallel for + for (int64_t sr = 0; sr < n_source_regions_; sr++) { + source_regions_.volume(sr) = 0.0; + } +#pragma omp parallel for + for (int64_t se = 0; se < n_source_elements_; se++) { + source_regions_.scalar_flux_new(se) = 0.0; + } } void FlatSourceDomain::accumulate_iteration_flux() { #pragma omp parallel for for (int64_t se = 0; se < n_source_elements_; se++) { - scalar_flux_final_[se] += scalar_flux_new_[se]; + source_regions_.scalar_flux_final(se) += + source_regions_.scalar_flux_new(se); } } @@ -137,40 +120,40 @@ void FlatSourceDomain::update_neutron_source(double k_eff) // Add scattering source #pragma omp parallel for - for (int sr = 0; sr < n_source_regions_; sr++) { - int material = material_[sr]; + for (int64_t sr = 0; sr < n_source_regions_; sr++) { + int material = source_regions_.material(sr); for (int g_out = 0; g_out < negroups_; g_out++) { double sigma_t = sigma_t_[material * negroups_ + g_out]; double scatter_source = 0.0; for (int g_in = 0; g_in < negroups_; g_in++) { - double scalar_flux = scalar_flux_old_[sr * negroups_ + g_in]; + double scalar_flux = source_regions_.scalar_flux_old(sr, g_in); double sigma_s = sigma_s_[material * negroups_ * negroups_ + g_out * negroups_ + g_in]; scatter_source += sigma_s * scalar_flux; } - source_[sr * negroups_ + g_out] = scatter_source / sigma_t; + source_regions_.source(sr, g_out) = scatter_source / sigma_t; } } // Add fission source #pragma omp parallel for - for (int sr = 0; sr < n_source_regions_; sr++) { - int material = material_[sr]; + for (int64_t sr = 0; sr < n_source_regions_; sr++) { + int material = source_regions_.material(sr); for (int g_out = 0; g_out < negroups_; g_out++) { double sigma_t = sigma_t_[material * negroups_ + g_out]; double fission_source = 0.0; for (int g_in = 0; g_in < negroups_; g_in++) { - double scalar_flux = scalar_flux_old_[sr * negroups_ + g_in]; + double scalar_flux = source_regions_.scalar_flux_old(sr, g_in); double nu_sigma_f = nu_sigma_f_[material * negroups_ + g_in]; double chi = chi_[material * negroups_ + g_out]; fission_source += nu_sigma_f * scalar_flux * chi; } - source_[sr * negroups_ + g_out] += + source_regions_.source(sr, g_out) += fission_source * inverse_k_eff / sigma_t; } } @@ -179,7 +162,7 @@ void FlatSourceDomain::update_neutron_source(double k_eff) if (settings::run_mode == RunMode::FIXED_SOURCE) { #pragma omp parallel for for (int64_t se = 0; se < n_source_elements_; se++) { - source_[se] += external_source_[se]; + source_regions_.source(se) += source_regions_.external_source(se); } } @@ -194,38 +177,42 @@ void FlatSourceDomain::normalize_scalar_flux_and_volumes( double volume_normalization_factor = 1.0 / (total_active_distance_per_iteration * simulation::current_batch); -// Normalize scalar flux to total distance travelled by all rays this iteration +// Normalize scalar flux to total distance travelled by all rays this +// iteration #pragma omp parallel for - for (int64_t se = 0; se < scalar_flux_new_.size(); se++) { - scalar_flux_new_[se] *= normalization_factor; + for (int64_t se = 0; se < n_source_elements_; se++) { + source_regions_.scalar_flux_new(se) *= normalization_factor; } // Accumulate cell-wise ray length tallies collected this iteration, then // update the simulation-averaged cell-wise volume estimates #pragma omp parallel for for (int64_t sr = 0; sr < n_source_regions_; sr++) { - volume_t_[sr] += volume_[sr]; - volume_naive_[sr] = volume_[sr] * normalization_factor; - volume_[sr] = volume_t_[sr] * volume_normalization_factor; + source_regions_.volume_t(sr) += source_regions_.volume(sr); + source_regions_.volume_naive(sr) = + source_regions_.volume(sr) * normalization_factor; + source_regions_.volume(sr) = + source_regions_.volume_t(sr) * volume_normalization_factor; } } void FlatSourceDomain::set_flux_to_flux_plus_source( - int64_t idx, double volume, int material, int g) + int64_t sr, double volume, int g) { - double sigma_t = sigma_t_[material * negroups_ + g]; - scalar_flux_new_[idx] /= (sigma_t * volume); - scalar_flux_new_[idx] += source_[idx]; + double sigma_t = sigma_t_[source_regions_.material(sr) * negroups_ + g]; + source_regions_.scalar_flux_new(sr, g) /= (sigma_t * volume); + source_regions_.scalar_flux_new(sr, g) += source_regions_.source(sr, g); } -void FlatSourceDomain::set_flux_to_old_flux(int64_t idx) +void FlatSourceDomain::set_flux_to_old_flux(int64_t sr, int g) { - scalar_flux_new_[idx] = scalar_flux_old_[idx]; + source_regions_.scalar_flux_new(sr, g) = + source_regions_.scalar_flux_old(sr, g); } -void FlatSourceDomain::set_flux_to_source(int64_t idx) +void FlatSourceDomain::set_flux_to_source(int64_t sr, int g) { - scalar_flux_new_[idx] = source_[idx]; + source_regions_.scalar_flux_new(sr, g) = source_regions_.source(sr, g); } // Combine transport flux contributions and flat source contributions from the @@ -235,20 +222,16 @@ int64_t FlatSourceDomain::add_source_to_scalar_flux() int64_t n_hits = 0; #pragma omp parallel for reduction(+ : n_hits) - for (int sr = 0; sr < n_source_regions_; sr++) { + for (int64_t sr = 0; sr < n_source_regions_; sr++) { - double volume_simulation_avg = volume_[sr]; - double volume_iteration = volume_naive_[sr]; + double volume_simulation_avg = source_regions_.volume(sr); + double volume_iteration = source_regions_.volume_naive(sr); // Increment the number of hits if cell was hit this iteration if (volume_iteration) { n_hits++; } - // Check if an external source is present in this source region - bool external_source_present = - external_source_present_.size() && external_source_present_[sr]; - // The volume treatment depends on the volume estimator type // and whether or not an external source is present in the cell. double volume; @@ -260,7 +243,7 @@ int64_t FlatSourceDomain::add_source_to_scalar_flux() volume = volume_simulation_avg; break; case RandomRayVolumeEstimator::HYBRID: - if (external_source_present) { + if (source_regions_.external_source_present(sr)) { volume = volume_iteration; } else { volume = volume_simulation_avg; @@ -270,17 +253,14 @@ int64_t FlatSourceDomain::add_source_to_scalar_flux() fatal_error("Invalid volume estimator type"); } - int material = material_[sr]; for (int g = 0; g < negroups_; g++) { - int64_t idx = (sr * negroups_) + g; - // There are three scenarios we need to consider: if (volume_iteration > 0.0) { // 1. If the FSR was hit this iteration, then the new flux is equal to // the flat source from the previous iteration plus the contributions // from rays passing through the source region (computed during the // transport sweep) - set_flux_to_flux_plus_source(idx, volume, material, g); + set_flux_to_flux_plus_source(sr, volume, g); } else if (volume_simulation_avg > 0.0) { // 2. If the FSR was not hit this iteration, but has been hit some // previous iteration, then we need to make a choice about what @@ -294,10 +274,10 @@ int64_t FlatSourceDomain::add_source_to_scalar_flux() // in the cell we will use the previous iteration's flux estimate. This // injects a small degree of correlation into the simulation, but this // is going to be trivial when the miss rate is a few percent or less. - if (external_source_present) { - set_flux_to_old_flux(idx); + if (source_regions_.external_source_present(sr)) { + set_flux_to_old_flux(sr, g); } else { - set_flux_to_source(idx); + set_flux_to_source(sr, g); } } // If the FSR was not hit this iteration, and it has never been hit in @@ -322,24 +302,25 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const vector p(n_source_regions_, 0.0f); #pragma omp parallel for reduction(+ : fission_rate_old, fission_rate_new) - for (int sr = 0; sr < n_source_regions_; sr++) { + for (int64_t sr = 0; sr < n_source_regions_; sr++) { // If simulation averaged volume is zero, don't include this cell - double volume = volume_[sr]; + double volume = source_regions_.volume(sr); if (volume == 0.0) { continue; } - int material = material_[sr]; + int material = source_regions_.material(sr); double sr_fission_source_old = 0; double sr_fission_source_new = 0; for (int g = 0; g < negroups_; g++) { - int64_t idx = (sr * negroups_) + g; double nu_sigma_f = nu_sigma_f_[material * negroups_ + g]; - sr_fission_source_old += nu_sigma_f * scalar_flux_old_[idx]; - sr_fission_source_new += nu_sigma_f * scalar_flux_new_[idx]; + sr_fission_source_old += + nu_sigma_f * source_regions_.scalar_flux_old(sr, g); + sr_fission_source_new += + nu_sigma_f * source_regions_.scalar_flux_new(sr, g); } // Compute total fission rates in FSR @@ -361,7 +342,7 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const double inverse_sum = 1 / fission_rate_new; #pragma omp parallel for reduction(+ : H) - for (int sr = 0; sr < n_source_regions_; sr++) { + for (int64_t sr = 0; sr < n_source_regions_; sr++) { // Only if FSR has non-negative and non-zero fission source if (p[sr] > 0.0f) { // Normalize to total weight of bank sites. p_i for better performance @@ -419,11 +400,11 @@ void FlatSourceDomain::convert_source_regions_to_tallies() // Attempt to generate mapping for all source regions #pragma omp parallel for - for (int sr = 0; sr < n_source_regions_; sr++) { + for (int64_t sr = 0; sr < n_source_regions_; sr++) { // If this source region has not been hit by a ray yet, then // we aren't going to be able to map it, so skip it. - if (!position_recorded_[sr]) { + if (!source_regions_.position_recorded(sr)) { all_source_regions_mapped = false; continue; } @@ -432,8 +413,8 @@ void FlatSourceDomain::convert_source_regions_to_tallies() // crossing through this source region is used to estabilish // the spatial location of the source region Particle p; - p.r() = position_[sr]; - p.r_last() = position_[sr]; + p.r() = source_regions_.position(sr); + p.r_last() = source_regions_.position(sr); bool found = exhaustive_find_cell(p); // Loop over energy groups (so as to support energy filters) @@ -449,7 +430,7 @@ void FlatSourceDomain::convert_source_regions_to_tallies() // If this task has already been populated, we don't need to do // it again. - if (tally_task_[source_element].size() > 0) { + if (source_regions_.tally_task(sr, g).size() > 0) { continue; } @@ -479,11 +460,11 @@ void FlatSourceDomain::convert_source_regions_to_tallies() // If a valid tally, filter, and score combination has been found, // then add it to the list of tally tasks for this source element. TallyTask task(i_tally, filter_index, score_index, score_bin); - tally_task_[source_element].push_back(task); + source_regions_.tally_task(sr, g).push_back(task); // Also add this task to the list of volume tasks for this source // region. - volume_task_[sr].insert(task); + source_regions_.volume_task(sr).insert(task); } } } @@ -529,13 +510,13 @@ double FlatSourceDomain::compute_fixed_source_normalization_factor() const // total external source strength in the simulation. double simulation_external_source_strength = 0.0; #pragma omp parallel for reduction(+ : simulation_external_source_strength) - for (int sr = 0; sr < n_source_regions_; sr++) { - int material = material_[sr]; - double volume = volume_[sr] * simulation_volume_; + for (int64_t sr = 0; sr < n_source_regions_; sr++) { + int material = source_regions_.material(sr); + double volume = source_regions_.volume(sr) * simulation_volume_; for (int g = 0; g < negroups_; g++) { double sigma_t = sigma_t_[material * negroups_ + g]; simulation_external_source_strength += - external_source_[sr * negroups_ + g] * sigma_t * volume; + source_regions_.external_source(sr, g) * sigma_t * volume; } } @@ -576,7 +557,7 @@ void FlatSourceDomain::random_ray_tally() // element, we check if there are any scores needed and apply // them. #pragma omp parallel for - for (int sr = 0; sr < n_source_regions_; sr++) { + for (int64_t sr = 0; sr < n_source_regions_; sr++) { // The fsr.volume_ is the unitless fractional simulation averaged volume // (i.e., it is the FSR's fraction of the overall simulation volume). The // simulation_volume_ is the total 3D physical volume in cm^3 of the @@ -586,15 +567,15 @@ void FlatSourceDomain::random_ray_tally() // eigenvalue solves, but useful in fixed source solves for returning the // flux shape with a magnitude that makes sense relative to the fixed // source strength. - double volume = volume_[sr] * simulation_volume_; + double volume = source_regions_.volume(sr) * simulation_volume_; - double material = material_[sr]; + double material = source_regions_.material(sr); for (int g = 0; g < negroups_; g++) { - int idx = sr * negroups_ + g; - double flux = scalar_flux_new_[idx] * source_normalization_factor; + double flux = + source_regions_.scalar_flux_new(sr, g) * source_normalization_factor; // Determine numerical score value - for (auto& task : tally_task_[idx]) { + for (auto& task : source_regions_.tally_task(sr, g)) { double score; switch (task.score_type) { @@ -637,7 +618,7 @@ void FlatSourceDomain::random_ray_tally() // for normalizing the flux. We store this volume in a separate tensor. // We only contribute to each volume tally bin once per FSR. if (volume_normalized_flux_tallies_) { - for (const auto& task : volume_task_[sr]) { + for (const auto& task : source_regions_.volume_task(sr)) { if (task.score_type == SCORE_FLUX) { #pragma omp atomic tally_volumes_[task.tally_idx](task.filter_idx, task.score_idx) += @@ -676,7 +657,6 @@ void FlatSourceDomain::random_ray_tally() void FlatSourceDomain::all_reduce_replicated_source_regions() { #ifdef OPENMC_MPI - // If we only have 1 MPI rank, no need // to reduce anything. if (mpi::n_procs <= 1) @@ -684,80 +664,15 @@ void FlatSourceDomain::all_reduce_replicated_source_regions() simulation::time_bank_sendrecv.start(); - // The "position_recorded" variable needs to be allreduced (and maxed), - // as whether or not a cell was hit will affect some decisions in how the - // source is calculated in the next iteration so as to avoid dividing - // by zero. We take the max rather than the sum as the hit values are - // expected to be zero or 1. - MPI_Allreduce(MPI_IN_PLACE, position_recorded_.data(), n_source_regions_, - MPI_INT, MPI_MAX, mpi::intracomm); - - // The position variable is more complicated to reduce than the others, - // as we do not want the sum of all positions in each cell, rather, we - // want to just pick any single valid position. Thus, we perform a gather - // and then pick the first valid position we find for all source regions - // that have had a position recorded. This operation does not need to - // be broadcast back to other ranks, as this value is only used for the - // tally conversion operation, which is only performed on the master rank. - // While this is expensive, it only needs to be done for active batches, - // and only if we have not mapped all the tallies yet. Once tallies are - // fully mapped, then the position vector is fully populated, so this - // operation can be skipped. - // First, we broadcast the fully mapped tally status variable so that // all ranks are on the same page int mapped_all_tallies_i = static_cast(mapped_all_tallies_); MPI_Bcast(&mapped_all_tallies_i, 1, MPI_INT, 0, mpi::intracomm); - // Then, we perform the gather of position data, if needed - if (simulation::current_batch > settings::n_inactive && - !mapped_all_tallies_i) { + bool reduce_position = + simulation::current_batch > settings::n_inactive && !mapped_all_tallies_i; - // Master rank will gather results and pick valid positions - if (mpi::master) { - // Initialize temporary vector for receiving positions - vector> all_position; - all_position.resize(mpi::n_procs); - for (int i = 0; i < mpi::n_procs; i++) { - all_position[i].resize(n_source_regions_); - } - - // Copy master rank data into gathered vector for convenience - all_position[0] = position_; - - // Receive all data into gather vector - for (int i = 1; i < mpi::n_procs; i++) { - MPI_Recv(all_position[i].data(), n_source_regions_ * 3, MPI_DOUBLE, i, - 0, mpi::intracomm, MPI_STATUS_IGNORE); - } - - // Scan through gathered data and pick first valid cell posiiton - for (int sr = 0; sr < n_source_regions_; sr++) { - if (position_recorded_[sr] == 1) { - for (int i = 0; i < mpi::n_procs; i++) { - if (all_position[i][sr].x != 0.0 || all_position[i][sr].y != 0.0 || - all_position[i][sr].z != 0.0) { - position_[sr] = all_position[i][sr]; - break; - } - } - } - } - } else { - // Other ranks just send in their data - MPI_Send(position_.data(), n_source_regions_ * 3, MPI_DOUBLE, 0, 0, - mpi::intracomm); - } - } - - // For the rest of the source region data, we simply perform an all reduce, - // as these values will be needed on all ranks for transport during the - // next iteration. - MPI_Allreduce(MPI_IN_PLACE, volume_.data(), n_source_regions_, MPI_DOUBLE, - MPI_SUM, mpi::intracomm); - - MPI_Allreduce(MPI_IN_PLACE, scalar_flux_new_.data(), n_source_elements_, - MPI_DOUBLE, MPI_SUM, mpi::intracomm); + source_regions_.mpi_sync_ranks(reduce_position); simulation::time_bank_sendrecv.stop(); #endif @@ -766,7 +681,7 @@ void FlatSourceDomain::all_reduce_replicated_source_regions() double FlatSourceDomain::evaluate_flux_at_point( Position r, int64_t sr, int g) const { - return scalar_flux_final_[sr * negroups_ + g] / + return source_regions_.scalar_flux_final(sr, g) / (settings::n_batches - settings::n_inactive); } @@ -894,7 +809,7 @@ void FlatSourceDomain::output_to_vtk() const std::fprintf(plot, "SCALARS Materials int\n"); std::fprintf(plot, "LOOKUP_TABLE default\n"); for (int fsr : voxel_indices) { - int mat = material_[fsr]; + int mat = source_regions_.material(fsr); mat = convert_to_big_endian(mat); std::fwrite(&mat, sizeof(int), 1, plot); } @@ -906,7 +821,7 @@ void FlatSourceDomain::output_to_vtk() const int64_t fsr = voxel_indices[i]; float total_fission = 0.0; - int mat = material_[fsr]; + int mat = source_regions_.material(fsr); for (int g = 0; g < negroups_; g++) { int64_t source_element = fsr * negroups_ + g; float flux = evaluate_flux_at_point(voxel_positions[i], fsr, g); @@ -922,16 +837,16 @@ void FlatSourceDomain::output_to_vtk() const } void FlatSourceDomain::apply_external_source_to_source_region( - Discrete* discrete, double strength_factor, int64_t source_region) + Discrete* discrete, double strength_factor, int64_t sr) { - external_source_present_[source_region] = true; + source_regions_.external_source_present(sr) = 1; const auto& discrete_energies = discrete->x(); const auto& discrete_probs = discrete->prob(); for (int i = 0; i < discrete_energies.size(); i++) { int g = data::mg.get_group_index(discrete_energies[i]); - external_source_[source_region * negroups_ + g] += + source_regions_.external_source(sr, g) += discrete_probs[i] * strength_factor; } } @@ -983,8 +898,8 @@ void FlatSourceDomain::count_external_source_regions() { n_external_source_regions_ = 0; #pragma omp parallel for reduction(+ : n_external_source_regions_) - for (int sr = 0; sr < n_source_regions_; sr++) { - if (external_source_present_[sr]) { + for (int64_t sr = 0; sr < n_source_regions_; sr++) { + if (source_regions_.external_source_present(sr)) { n_external_source_regions_++; } } @@ -1029,18 +944,17 @@ void FlatSourceDomain::convert_external_sources() // Divide the fixed source term by sigma t (to save time when applying each // iteration) #pragma omp parallel for - for (int sr = 0; sr < n_source_regions_; sr++) { - int material = material_[sr]; + for (int64_t sr = 0; sr < n_source_regions_; sr++) { for (int g = 0; g < negroups_; g++) { - double sigma_t = sigma_t_[material * negroups_ + g]; - external_source_[sr * negroups_ + g] /= sigma_t; + double sigma_t = sigma_t_[source_regions_.material(sr) * negroups_ + g]; + source_regions_.external_source(sr, g) /= sigma_t; } } } void FlatSourceDomain::flux_swap() { - scalar_flux_old_.swap(scalar_flux_new_); + source_regions_.flux_swap(); } void FlatSourceDomain::flatten_xs() @@ -1096,17 +1010,16 @@ void FlatSourceDomain::set_adjoint_sources(const vector& forward_flux) // so we must convert it to a "per batch" quantity #pragma omp parallel for for (int64_t se = 0; se < n_source_elements_; se++) { - external_source_[se] = 1.0 / forward_flux[se]; + source_regions_.external_source(se) = 1.0 / forward_flux[se]; } // Divide the fixed source term by sigma t (to save time when applying each // iteration) #pragma omp parallel for - for (int sr = 0; sr < n_source_regions_; sr++) { - int material = material_[sr]; + for (int64_t sr = 0; sr < n_source_regions_; sr++) { for (int g = 0; g < negroups_; g++) { - double sigma_t = sigma_t_[material * negroups_ + g]; - external_source_[sr * negroups_ + g] /= sigma_t; + double sigma_t = sigma_t_[source_regions_.material(sr) * negroups_ + g]; + source_regions_.external_source(sr, g) /= sigma_t; } } } @@ -1129,4 +1042,15 @@ void FlatSourceDomain::transpose_scattering_matrix() } } +void FlatSourceDomain::serialize_final_fluxes(vector& flux) +{ + // Ensure array is correct size + flux.resize(n_source_regions_ * negroups_); +// Serialize the final fluxes for output +#pragma omp parallel for + for (int64_t se = 0; se < n_source_elements_; se++) { + flux[se] = source_regions_.scalar_flux_final(se); + } +} + } // namespace openmc diff --git a/src/random_ray/linear_source_domain.cpp b/src/random_ray/linear_source_domain.cpp index f3f6fa068..fb6502a3e 100644 --- a/src/random_ray/linear_source_domain.cpp +++ b/src/random_ray/linear_source_domain.cpp @@ -20,33 +20,17 @@ namespace openmc { // LinearSourceDomain implementation //============================================================================== -LinearSourceDomain::LinearSourceDomain() : FlatSourceDomain() -{ - // First order spatial moment of the scalar flux - flux_moments_old_.assign(n_source_elements_, {0.0, 0.0, 0.0}); - flux_moments_new_.assign(n_source_elements_, {0.0, 0.0, 0.0}); - flux_moments_t_.assign(n_source_elements_, {0.0, 0.0, 0.0}); - // Source gradients given by M inverse multiplied by source moments - source_gradients_.assign(n_source_elements_, {0.0, 0.0, 0.0}); - - centroid_.assign(n_source_regions_, {0.0, 0.0, 0.0}); - centroid_iteration_.assign(n_source_regions_, {0.0, 0.0, 0.0}); - centroid_t_.assign(n_source_regions_, {0.0, 0.0, 0.0}); - mom_matrix_.assign(n_source_regions_, {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}); - mom_matrix_t_.assign(n_source_regions_, {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}); -} - void LinearSourceDomain::batch_reset() { FlatSourceDomain::batch_reset(); #pragma omp parallel for - for (int64_t se = 0; se < n_source_elements_; se++) { - flux_moments_new_[se] = {0.0, 0.0, 0.0}; + for (int64_t sr = 0; sr < n_source_regions_; sr++) { + source_regions_.centroid_iteration(sr) = {0.0, 0.0, 0.0}; + source_regions_.mom_matrix(sr) = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; } #pragma omp parallel for - for (int64_t sr = 0; sr < n_source_regions_; sr++) { - centroid_iteration_[sr] = {0.0, 0.0, 0.0}; - mom_matrix_[sr] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + for (int64_t se = 0; se < n_source_elements_; se++) { + source_regions_.flux_moments_new(se) = {0.0, 0.0, 0.0}; } } @@ -57,10 +41,9 @@ void LinearSourceDomain::update_neutron_source(double k_eff) double inverse_k_eff = 1.0 / k_eff; #pragma omp parallel for - for (int sr = 0; sr < n_source_regions_; sr++) { - - int material = material_[sr]; - MomentMatrix invM = mom_matrix_[sr].inverse(); + for (int64_t sr = 0; sr < n_source_regions_; sr++) { + int material = source_regions_.material(sr); + MomentMatrix invM = source_regions_.mom_matrix(sr).inverse(); for (int g_out = 0; g_out < negroups_; g_out++) { double sigma_t = sigma_t_[material * negroups_ + g_out]; @@ -72,8 +55,8 @@ void LinearSourceDomain::update_neutron_source(double k_eff) for (int g_in = 0; g_in < negroups_; g_in++) { // Handles for the flat and linear components of the flux - double flux_flat = scalar_flux_old_[sr * negroups_ + g_in]; - MomentArray flux_linear = flux_moments_old_[sr * negroups_ + g_in]; + double flux_flat = source_regions_.scalar_flux_old(sr, g_in); + MomentArray flux_linear = source_regions_.flux_moments_old(sr, g_in); // Handles for cross sections double sigma_s = @@ -89,7 +72,7 @@ void LinearSourceDomain::update_neutron_source(double k_eff) } // Compute the flat source term - source_[sr * negroups_ + g_out] = + source_regions_.source(sr, g_out) = (scatter_flat + fission_flat * inverse_k_eff) / sigma_t; // Compute the linear source terms @@ -97,7 +80,7 @@ void LinearSourceDomain::update_neutron_source(double k_eff) // are not well known, we will leave the source gradients as zero // so as to avoid causing any numerical instability. if (simulation::current_batch > 10) { - source_gradients_[sr * negroups_ + g_out] = + source_regions_.source_gradients(sr, g_out) = invM * ((scatter_linear + fission_linear * inverse_k_eff) / sigma_t); } } @@ -107,7 +90,7 @@ void LinearSourceDomain::update_neutron_source(double k_eff) // Add external source to flat source term if in fixed source mode #pragma omp parallel for for (int64_t se = 0; se < n_source_elements_; se++) { - source_[se] += external_source_[se]; + source_regions_.source(se) += source_regions_.external_source(se); } } @@ -123,48 +106,44 @@ void LinearSourceDomain::normalize_scalar_flux_and_volumes( // Normalize flux to total distance travelled by all rays this iteration #pragma omp parallel for - for (int64_t se = 0; se < scalar_flux_new_.size(); se++) { - scalar_flux_new_[se] *= normalization_factor; - flux_moments_new_[se] *= normalization_factor; + for (int64_t se = 0; se < n_source_elements_; se++) { + source_regions_.scalar_flux_new(se) *= normalization_factor; + source_regions_.flux_moments_new(se) *= normalization_factor; } // Accumulate cell-wise ray length tallies collected this iteration, then // update the simulation-averaged cell-wise volume estimates #pragma omp parallel for for (int64_t sr = 0; sr < n_source_regions_; sr++) { - centroid_t_[sr] += centroid_iteration_[sr]; - mom_matrix_t_[sr] += mom_matrix_[sr]; - volume_t_[sr] += volume_[sr]; - volume_naive_[sr] = volume_[sr] * normalization_factor; - volume_[sr] = volume_t_[sr] * volume_normalization_factor; - if (volume_t_[sr] > 0.0) { - double inv_volume = 1.0 / volume_t_[sr]; - centroid_[sr] = centroid_t_[sr]; - centroid_[sr] *= inv_volume; - mom_matrix_[sr] = mom_matrix_t_[sr]; - mom_matrix_[sr] *= inv_volume; + source_regions_.centroid_t(sr) += source_regions_.centroid_iteration(sr); + source_regions_.mom_matrix_t(sr) += source_regions_.mom_matrix(sr); + source_regions_.volume_t(sr) += source_regions_.volume(sr); + source_regions_.volume_naive(sr) = + source_regions_.volume(sr) * normalization_factor; + source_regions_.volume(sr) = + source_regions_.volume_t(sr) * volume_normalization_factor; + if (source_regions_.volume_t(sr) > 0.0) { + double inv_volume = 1.0 / source_regions_.volume_t(sr); + source_regions_.centroid(sr) = source_regions_.centroid_t(sr); + source_regions_.centroid(sr) *= inv_volume; + source_regions_.mom_matrix(sr) = source_regions_.mom_matrix_t(sr); + source_regions_.mom_matrix(sr) *= inv_volume; } } } void LinearSourceDomain::set_flux_to_flux_plus_source( - int64_t idx, double volume, int material, int g) + int64_t sr, double volume, int g) { - scalar_flux_new_[idx] /= volume; - scalar_flux_new_[idx] += source_[idx]; - flux_moments_new_[idx] *= (1.0 / volume); + source_regions_.scalar_flux_new(sr, g) /= volume; + source_regions_.scalar_flux_new(sr, g) += source_regions_.source(sr, g); + source_regions_.flux_moments_new(sr, g) *= (1.0 / volume); } -void LinearSourceDomain::set_flux_to_old_flux(int64_t idx) +void LinearSourceDomain::set_flux_to_old_flux(int64_t sr, int g) { - scalar_flux_new_[idx] = scalar_flux_old_[idx]; - flux_moments_new_[idx] = flux_moments_old_[idx]; -} - -void LinearSourceDomain::flux_swap() -{ - FlatSourceDomain::flux_swap(); - flux_moments_old_.swap(flux_moments_new_); + source_regions_.scalar_flux_new(g) = source_regions_.scalar_flux_old(g); + source_regions_.flux_moments_new(g) = source_regions_.flux_moments_old(g); } void LinearSourceDomain::accumulate_iteration_flux() @@ -175,48 +154,20 @@ void LinearSourceDomain::accumulate_iteration_flux() // Accumulate scalar flux moments #pragma omp parallel for for (int64_t se = 0; se < n_source_elements_; se++) { - flux_moments_t_[se] += flux_moments_new_[se]; + source_regions_.flux_moments_t(se) += source_regions_.flux_moments_new(se); } } -void LinearSourceDomain::all_reduce_replicated_source_regions() -{ -#ifdef OPENMC_MPI - FlatSourceDomain::all_reduce_replicated_source_regions(); - simulation::time_bank_sendrecv.start(); - - // We are going to assume we can safely cast Position, MomentArray, - // and MomentMatrix to contiguous arrays of doubles for the MPI - // allreduce operation. This is a safe assumption as typically - // compilers will at most pad to 8 byte boundaries. If a new FP32 MomentArray - // type is introduced, then there will likely be padding, in which case this - // function will need to become more complex. - if (sizeof(MomentArray) != 3 * sizeof(double) || - sizeof(MomentMatrix) != 6 * sizeof(double)) { - fatal_error("Unexpected buffer padding in linear source domain reduction."); - } - - MPI_Allreduce(MPI_IN_PLACE, static_cast(flux_moments_new_.data()), - n_source_elements_ * 3, MPI_DOUBLE, MPI_SUM, mpi::intracomm); - MPI_Allreduce(MPI_IN_PLACE, static_cast(mom_matrix_.data()), - n_source_regions_ * 6, MPI_DOUBLE, MPI_SUM, mpi::intracomm); - MPI_Allreduce(MPI_IN_PLACE, static_cast(centroid_iteration_.data()), - n_source_regions_ * 3, MPI_DOUBLE, MPI_SUM, mpi::intracomm); - - simulation::time_bank_sendrecv.stop(); -#endif -} - double LinearSourceDomain::evaluate_flux_at_point( Position r, int64_t sr, int g) const { double phi_flat = FlatSourceDomain::evaluate_flux_at_point(r, sr, g); - Position local_r = r - centroid_[sr]; - MomentArray phi_linear = flux_moments_t_[sr * negroups_ + g]; + Position local_r = r - source_regions_.centroid(sr); + MomentArray phi_linear = source_regions_.flux_moments_t(sr, g); phi_linear *= 1.0 / (settings::n_batches - settings::n_inactive); - MomentMatrix invM = mom_matrix_[sr].inverse(); + MomentMatrix invM = source_regions_.mom_matrix(sr).inverse(); MomentArray phi_solved = invM * phi_linear; return phi_flat + phi_solved.dot(local_r); diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp index 7a359f356..9c2a890ec 100644 --- a/src/random_ray/random_ray.cpp +++ b/src/random_ray/random_ray.cpp @@ -309,11 +309,9 @@ void RandomRay::attenuate_flux_flat_source(double distance, bool is_active) int i_cell = lowest_coord().cell; // The source region is the spatial region index - int64_t source_region = - domain_->source_region_offsets_[i_cell] + cell_instance(); + int64_t sr = domain_->source_region_offsets_[i_cell] + cell_instance(); // The source element is the energy-specific region index - int64_t source_element = source_region * negroups_; int material = this->material(); // MOC incoming flux attenuation + source contribution/attenuation equation @@ -322,7 +320,7 @@ void RandomRay::attenuate_flux_flat_source(double distance, bool is_active) float tau = sigma_t * distance; float exponential = cjosey_exponential(tau); // exponential = 1 - exp(-tau) float new_delta_psi = - (angular_flux_[g] - domain_->source_[source_element + g]) * exponential; + (angular_flux_[g] - domain_->source_regions_.source(sr, g)) * exponential; delta_psi_[g] = new_delta_psi; angular_flux_[g] -= new_delta_psi; } @@ -332,28 +330,28 @@ void RandomRay::attenuate_flux_flat_source(double distance, bool is_active) if (is_active) { // Aquire lock for source region - domain_->lock_[source_region].lock(); + domain_->source_regions_.lock(sr).lock(); // Accumulate delta psi into new estimate of source region flux for // this iteration for (int g = 0; g < negroups_; g++) { - domain_->scalar_flux_new_[source_element + g] += delta_psi_[g]; + domain_->source_regions_.scalar_flux_new(sr, g) += delta_psi_[g]; } // Accomulate volume (ray distance) into this iteration's estimate // of the source region's volume - domain_->volume_[source_region] += distance; + domain_->source_regions_.volume(sr) += distance; // Tally valid position inside the source region (e.g., midpoint of // the ray) if not done already - if (!domain_->position_recorded_[source_region]) { + if (!domain_->source_regions_.position_recorded(sr)) { Position midpoint = r() + u() * (distance / 2.0); - domain_->position_[source_region] = midpoint; - domain_->position_recorded_[source_region] = 1; + domain_->source_regions_.position(sr) = midpoint; + domain_->source_regions_.position_recorded(sr) = 1; } // Release lock - domain_->lock_[source_region].unlock(); + domain_->source_regions_.lock(sr).unlock(); } } @@ -373,14 +371,12 @@ void RandomRay::attenuate_flux_linear_source(double distance, bool is_active) int i_cell = lowest_coord().cell; // The source region is the spatial region index - int64_t source_region = - domain_->source_region_offsets_[i_cell] + cell_instance(); + int64_t sr = domain_->source_region_offsets_[i_cell] + cell_instance(); // The source element is the energy-specific region index - int64_t source_element = source_region * negroups_; int material = this->material(); - Position& centroid = domain->centroid_[source_region]; + Position& centroid = domain_->source_regions_.centroid(sr); Position midpoint = r() + u() * (distance / 2.0); // Determine the local position of the midpoint and the ray origin @@ -393,7 +389,7 @@ void RandomRay::attenuate_flux_linear_source(double distance, bool is_active) // be no estimate of its centroid. We detect this by checking if it has // any accumulated volume. If its volume is zero, just use the midpoint // of the ray as the region's centroid. - if (domain->volume_t_[source_region]) { + if (domain_->source_regions_.volume_t(sr)) { rm_local = midpoint - centroid; r0_local = r() - centroid; } else { @@ -420,9 +416,10 @@ void RandomRay::attenuate_flux_linear_source(double distance, bool is_active) // calculated from the source gradients dot product with local centroid // and direction, respectively. float spatial_source = - domain_->source_[source_element + g] + - rm_local.dot(domain->source_gradients_[source_element + g]); - float dir_source = u().dot(domain->source_gradients_[source_element + g]); + domain_->source_regions_.source(sr, g) + + rm_local.dot(domain_->source_regions_.source_gradients(sr, g)); + float dir_source = + u().dot(domain_->source_regions_.source_gradients(sr, g)); float gn = exponentialG(tau); float f1 = 1.0f - tau * gn; @@ -465,33 +462,33 @@ void RandomRay::attenuate_flux_linear_source(double distance, bool is_active) rm_local, u(), distance); // Aquire lock for source region - domain_->lock_[source_region].lock(); + domain_->source_regions_.lock(sr).lock(); // Accumulate deltas into the new estimate of source region flux for this // iteration for (int g = 0; g < negroups_; g++) { - domain_->scalar_flux_new_[source_element + g] += delta_psi_[g]; - domain->flux_moments_new_[source_element + g] += delta_moments_[g]; + domain_->source_regions_.scalar_flux_new(sr, g) += delta_psi_[g]; + domain_->source_regions_.flux_moments_new(sr, g) += delta_moments_[g]; } // Accumulate the volume (ray segment distance), centroid, and spatial // momement estimates into the running totals for the iteration for this // source region. The centroid and spatial momements estimates are scaled by // the ray segment length as part of length averaging of the estimates. - domain_->volume_[source_region] += distance; - domain->centroid_iteration_[source_region] += midpoint * distance; + domain_->source_regions_.volume(sr) += distance; + domain_->source_regions_.centroid_iteration(sr) += midpoint * distance; moment_matrix_estimate *= distance; - domain->mom_matrix_[source_region] += moment_matrix_estimate; + domain_->source_regions_.mom_matrix(sr) += moment_matrix_estimate; // Tally valid position inside the source region (e.g., midpoint of // the ray) if not done already - if (!domain_->position_recorded_[source_region]) { - domain_->position_[source_region] = midpoint; - domain_->position_recorded_[source_region] = 1; + if (!domain_->source_regions_.position_recorded(sr)) { + domain_->source_regions_.position(sr) = midpoint; + domain_->source_regions_.position_recorded(sr) = 1; } // Release lock - domain_->lock_[source_region].unlock(); + domain_->source_regions_.lock(sr).unlock(); } } @@ -537,11 +534,10 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain) // Initialize ray's starting angular flux to starting location's isotropic // source int i_cell = lowest_coord().cell; - int64_t source_region_idx = - domain_->source_region_offsets_[i_cell] + cell_instance(); + int64_t sr = domain_->source_region_offsets_[i_cell] + cell_instance(); for (int g = 0; g < negroups_; g++) { - angular_flux_[g] = domain_->source_[source_region_idx * negroups_ + g]; + angular_flux_[g] = domain_->source_regions_.source(sr, g); } } diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index a8e6d1b82..c4f719b23 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -66,7 +66,7 @@ void openmc_run_random_ray() simulation::time_total.stop(); // Normalize and save the final forward flux - forward_flux = sim.domain()->scalar_flux_final_; + sim.domain()->serialize_final_fluxes(forward_flux); double source_normalization_factor = sim.domain()->compute_fixed_source_normalization_factor() / diff --git a/src/random_ray/source_region.cpp b/src/random_ray/source_region.cpp new file mode 100644 index 000000000..0a20af8a9 --- /dev/null +++ b/src/random_ray/source_region.cpp @@ -0,0 +1,236 @@ +#include "openmc/random_ray/source_region.h" + +#include "openmc/error.h" +#include "openmc/message_passing.h" + +namespace openmc { + +//============================================================================== +// SourceRegion implementation +//============================================================================== +SourceRegion::SourceRegion(int negroups, bool is_linear) +{ + if (settings::run_mode == RunMode::EIGENVALUE) { + // If in eigenvalue mode, set starting flux to guess of 1 + scalar_flux_old_.assign(negroups, 1.0); + } else { + // If in fixed source mode, set starting flux to guess of zero + // and initialize external source arrays + scalar_flux_old_.assign(negroups, 0.0); + external_source_.assign(negroups, 0.0); + } + + scalar_flux_new_.assign(negroups, 0.0); + source_.resize(negroups); + scalar_flux_final_.assign(negroups, 0.0); + + tally_task_.resize(negroups); + + if (is_linear) { + source_gradients_.resize(negroups); + flux_moments_old_.resize(negroups); + flux_moments_new_.resize(negroups); + flux_moments_t_.resize(negroups); + } +} + +//============================================================================== +// SourceRegionContainer implementation +//============================================================================== +void SourceRegionContainer::push_back(const SourceRegion& sr) +{ + n_source_regions_++; + + // Scalar fields + material_.push_back(sr.material_); + lock_.push_back(sr.lock_); + volume_.push_back(sr.volume_); + volume_t_.push_back(sr.volume_t_); + volume_naive_.push_back(sr.volume_naive_); + position_recorded_.push_back(sr.position_recorded_); + external_source_present_.push_back(sr.external_source_present_); + position_.push_back(sr.position_); + volume_task_.push_back(sr.volume_task_); + + // Only store these fields if is_linear_ is true + if (is_linear_) { + centroid_.push_back(sr.centroid_); + centroid_iteration_.push_back(sr.centroid_iteration_); + centroid_t_.push_back(sr.centroid_t_); + mom_matrix_.push_back(sr.mom_matrix_); + mom_matrix_t_.push_back(sr.mom_matrix_t_); + } + + // Energy-dependent fields + for (int g = 0; g < negroups_; ++g) { + scalar_flux_old_.push_back(sr.scalar_flux_old_[g]); + scalar_flux_new_.push_back(sr.scalar_flux_new_[g]); + scalar_flux_final_.push_back(sr.scalar_flux_final_[g]); + source_.push_back(sr.source_[g]); + if (settings::run_mode == RunMode::FIXED_SOURCE) { + external_source_.push_back(sr.external_source_[g]); + } + + // Only store these fields if is_linear_ is true + if (is_linear_) { + source_gradients_.push_back(sr.source_gradients_[g]); + flux_moments_old_.push_back(sr.flux_moments_old_[g]); + flux_moments_new_.push_back(sr.flux_moments_new_[g]); + flux_moments_t_.push_back(sr.flux_moments_t_[g]); + } + + // Tally tasks + tally_task_.emplace_back(sr.tally_task_[g]); + } +} + +void SourceRegionContainer::assign( + int n_source_regions, const SourceRegion& source_region) +{ + // Clear existing data + n_source_regions_ = 0; + material_.clear(); + lock_.clear(); + volume_.clear(); + volume_t_.clear(); + volume_naive_.clear(); + position_recorded_.clear(); + external_source_present_.clear(); + position_.clear(); + + if (is_linear_) { + centroid_.clear(); + centroid_iteration_.clear(); + centroid_t_.clear(); + mom_matrix_.clear(); + mom_matrix_t_.clear(); + } + + scalar_flux_old_.clear(); + scalar_flux_new_.clear(); + scalar_flux_final_.clear(); + source_.clear(); + external_source_.clear(); + + if (is_linear_) { + source_gradients_.clear(); + flux_moments_old_.clear(); + flux_moments_new_.clear(); + flux_moments_t_.clear(); + } + + tally_task_.clear(); + volume_task_.clear(); + + // Fill with copies of source_region + for (int i = 0; i < n_source_regions; ++i) { + push_back(source_region); + } +} + +void SourceRegionContainer::flux_swap() +{ + scalar_flux_old_.swap(scalar_flux_new_); + if (is_linear_) { + flux_moments_old_.swap(flux_moments_new_); + } +} + +void SourceRegionContainer::mpi_sync_ranks(bool reduce_position) +{ +#ifdef OPENMC_MPI + + // The "position_recorded" variable needs to be allreduced (and maxed), + // as whether or not a cell was hit will affect some decisions in how the + // source is calculated in the next iteration so as to avoid dividing + // by zero. We take the max rather than the sum as the hit values are + // expected to be zero or 1. + MPI_Allreduce(MPI_IN_PLACE, position_recorded_.data(), 1, MPI_INT, MPI_MAX, + mpi::intracomm); + + // The position variable is more complicated to reduce than the others, + // as we do not want the sum of all positions in each cell, rather, we + // want to just pick any single valid position. Thus, we perform a gather + // and then pick the first valid position we find for all source regions + // that have had a position recorded. This operation does not need to + // be broadcast back to other ranks, as this value is only used for the + // tally conversion operation, which is only performed on the master rank. + // While this is expensive, it only needs to be done for active batches, + // and only if we have not mapped all the tallies yet. Once tallies are + // fully mapped, then the position vector is fully populated, so this + // operation can be skipped. + + // Then, we perform the gather of position data, if needed + if (reduce_position) { + + // Master rank will gather results and pick valid positions + if (mpi::master) { + // Initialize temporary vector for receiving positions + vector> all_position; + all_position.resize(mpi::n_procs); + for (int i = 0; i < mpi::n_procs; i++) { + all_position[i].resize(n_source_regions_); + } + + // Copy master rank data into gathered vector for convenience + all_position[0] = position_; + + // Receive all data into gather vector + for (int i = 1; i < mpi::n_procs; i++) { + MPI_Recv(all_position[i].data(), n_source_regions_ * 3, MPI_DOUBLE, i, + 0, mpi::intracomm, MPI_STATUS_IGNORE); + } + + // Scan through gathered data and pick first valid cell posiiton + for (int64_t sr = 0; sr < n_source_regions_; sr++) { + if (position_recorded_[sr] == 1) { + for (int i = 0; i < mpi::n_procs; i++) { + if (all_position[i][sr].x != 0.0 || all_position[i][sr].y != 0.0 || + all_position[i][sr].z != 0.0) { + position_[sr] = all_position[i][sr]; + break; + } + } + } + } + } else { + // Other ranks just send in their data + MPI_Send(position_.data(), n_source_regions_ * 3, MPI_DOUBLE, 0, 0, + mpi::intracomm); + } + } + + // For the rest of the source region data, we simply perform an all reduce, + // as these values will be needed on all ranks for transport during the + // next iteration. + MPI_Allreduce(MPI_IN_PLACE, volume_.data(), n_source_regions_, MPI_DOUBLE, + MPI_SUM, mpi::intracomm); + + MPI_Allreduce(MPI_IN_PLACE, scalar_flux_new_.data(), + n_source_regions_ * negroups_, MPI_DOUBLE, MPI_SUM, mpi::intracomm); + + if (is_linear_) { + // We are going to assume we can safely cast Position, MomentArray, + // and MomentMatrix to contiguous arrays of doubles for the MPI + // allreduce operation. This is a safe assumption as typically + // compilers will at most pad to 8 byte boundaries. If a new FP32 + // MomentArray type is introduced, then there will likely be padding, in + // which case this function will need to become more complex. + if (sizeof(MomentArray) != 3 * sizeof(double) || + sizeof(MomentMatrix) != 6 * sizeof(double)) { + fatal_error( + "Unexpected buffer padding in linear source domain reduction."); + } + + MPI_Allreduce(MPI_IN_PLACE, static_cast(flux_moments_new_.data()), + n_source_regions_ * negroups_ * 3, MPI_DOUBLE, MPI_SUM, mpi::intracomm); + MPI_Allreduce(MPI_IN_PLACE, static_cast(mom_matrix_.data()), + n_source_regions_ * 6, MPI_DOUBLE, MPI_SUM, mpi::intracomm); + MPI_Allreduce(MPI_IN_PLACE, static_cast(centroid_iteration_.data()), + n_source_regions_ * 3, MPI_DOUBLE, MPI_SUM, mpi::intracomm); + } + +#endif +} + +} // namespace openmc \ No newline at end of file