diff --git a/CMakeLists.txt b/CMakeLists.txt index 61d2cc6c590..3fe310f99a1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -406,6 +406,8 @@ list(APPEND libopenmc_SOURCES src/progress_bar.cpp src/random_dist.cpp src/random_lcg.cpp + src/random_ray/decomposition_map.cpp + src/random_ray/ray_bank.cpp src/random_ray/random_ray_simulation.cpp src/random_ray/random_ray.cpp src/random_ray/flat_source_domain.cpp diff --git a/docs/source/methods/random_ray.rst b/docs/source/methods/random_ray.rst index 8bc2a0a1bf5..d8c886ea182 100644 --- a/docs/source/methods/random_ray.rst +++ b/docs/source/methods/random_ray.rst @@ -1100,14 +1100,176 @@ The adjoint external source will be computed for each source region in the simulation mesh, independent of any tallies. The adjoint external source is always flat, even when a linear scattering and fission source shape is used. -When in adjoint mode, all reported results (e.g., tallies, eigenvalues, etc.) -are derived from the adjoint flux, even when the physical meaning is not -necessarily obvious. These values are still reported, though we emphasize that -the primary use case for adjoint mode is for producing adjoint flux tallies to -support subsequent perturbation studies and weight window generation. Note -however that the adjoint :math:`k_{eff}` is statistically the same as the +When in adjoint mode, all reported results (e.g., tallies, eigenvalues, etc.) +are derived from the adjoint flux, even when the physical meaning is not +necessarily obvious. These values are still reported, though we emphasize that +the primary use case for adjoint mode is for producing adjoint flux tallies to +support subsequent perturbation studies and weight window generation. Note +however that the adjoint :math:`k_{eff}` is statistically the same as the forward :math:`k_{eff}`, despite the flux distributions taking different shapes. +-------------------- +Domain Decomposition +-------------------- + +To enable parallelisation and scalability beyond the resources of a single +computational node, a domain decomposition capability is available for the +random ray solver. + +~~~~~~~~~~~~~~~~~~~~ +Voronoi Tessellation +~~~~~~~~~~~~~~~~~~~~ + +The domain decomposition scheme distributes the source regions across multiple +MPI processes (ranks) according to a `capacity-constrained Voronoi tessellation +`_. Each MPI rank is responsible for the transport sweeps and +result tallying in one Voronoi region of the problem. Source regions are +assigned to MPI ranks (and thus Voronoi regions) using a formula that evaluates +the distance between the first intersection point :math:`\mathbf{x}` of a ray +with that source region and the Voronoi region centroid +:math:`\mathbf{c}_{\mathrm{rank}}`, combined with an additive weight +:math:`\omega_{\mathrm{rank}}`: + +.. math:: + :label: mpi_ownership + + \mathrm{rank}(\mathbf{x}) = \arg\min_{\text{rank}} + \left(\|\mathbf{c}_{\mathrm{rank}} - \mathbf{x}\|^2 - + \omega_{\mathrm{rank}} \right)\;\mathrm{.} + +The initial weight is zero and subsequently changed for load balancing. This +approach yields compact MPI rank subdomains. + +.. figure:: ../_images/c5g7_geometry.png + :width: 48% + :align: center + :figclass: align-center + + C5G7 geometry. + +.. figure:: ../_images/c5g7_voronoi.png + :width: 48% + :align: center + :figclass: align-center + + Voronoi decomposition of C5G7 source regions. + +In the OpenMC random ray implementation, the algorithm is not aware of the +source regions in the geometry a priori. Instead, source regions are discovered +dynamically as rays travel through the geometry and, once discovered, they get +added to a list of known source regions. Whenever a ray enters a previously +unknown source region, the responsible MPI rank is determined using the formula +given above. Ideally, the source region centroid would be used to assign +ownership unambiguously. However, centroid positions are not precalculated, and +ownership is instead decided based on the ray entry point. If a source region +happens to sit on a boundary between two Voronoi regions, it may be hit by rays +from both MPI ranks at different locations, and both MPI ranks may therefore +claim the same source region during the transport sweep. To resolve these +conflicts, after each transport sweep, ownership of contested source regions is +decided based on the estimated load of the ranks involved. + +Once each newly discovered source region has a unique owner rank, the ownership +information is shared across all MPI ranks and saved in a decomposition map. +This decomposition map is used to look up the responsible MPI rank every time a +ray enters a source region that has already been recorded, thereby avoiding the +need to evaluate Equation :eq:`mpi_ownership` again, which can be time-intensive +if many MPI ranks are present. + +~~~~~~~~~~~~~~~~~ +Ray Communication +~~~~~~~~~~~~~~~~~ + +When rays exit an MPI rank subdomain, they must be transmitted to their new +owner rank so that the transport can continue until the rays reach their +termination distance. Every time an MPI rank detects that a ray is leaving its +subdomain, the transport of that ray stops, and the ray attributes (angular flux +values, position, direction, distance traveled, etc.) are stored in a buffer. +Once each MPI rank has processed all rays in its subdomain, i.e. the rays have +either terminated or been moved into the buffer, all MPI ranks send their +buffered ray data in a bulk synchronous communication pattern to the new MPI +owner ranks. After communicating the ray data, each MPI rank reinitializes the +received rays with the transmitted data, and the rays continue traveling. This +communication pattern continues until all rays of a given batch have terminated. + +~~~~~~~~~~~~~~ +Load Balancing +~~~~~~~~~~~~~~ + +To ensure high parallel efficiency, it is crucial to assign each MPI rank +approximately the same amount of computational workload, such that the +individual MPI processes do not spend excessive time waiting at synchronization +steps (like the synchronous communication phase described above), while others +are still performing their calculations. In many computational fields, the +amount of work that is performed is fixed per source region. In such cases, the +overall load is simply a function of how many source regions are contained +within a given subdomain, with each source region requiring the same set of +operations to be performed. + +However, in random ray, the load for a given source region is much +more complex to determine. When a ray crosses a cell, the angular flux increment +:math:`\Delta \psi_{r,g}` for each energy group is calculated according to +Equation :eq:`delta_psi`. Additionally, to determine the length of the ray +through the cell (and which cell comes next), ray tracing operations are +performed. The frequency of these calculations and the cost of the ray trace +operations depend on the size, aspect ratio and definition of a given cell, +which can vary strongly across the simulation problem. To estimate the workload +associated with a given source region, an empirical formula has been set up that +accounts for 1) the number of ray crossings :math:`n_{\mathrm{hits}, i}` +in a source region :math:`i`, and 2) the number of surface ray trace operations +:math:`n_{\mathrm{RT}, i}` associated with the definition of that source +region: + +.. math:: + \mathrm{load}_{\mathrm{estimate}, i} = F_{r} \cdot \left(C_1 \cdot + n_{\mathrm{hits}, i} \cdot N_{G} + C_2 \cdot n_{\mathrm{RT}, i} \right)\; + \mathrm{.} + +The quantities :math:`n_{\mathrm{hits}, i}` and :math:`n_{\mathrm{RT}, i}` are +recorded throughout the simulation. Both contributions are weighted with factors +:math:`C_1` and :math:`C_2`, which represent the relative computational cost +of these operations. The values of these factors are set to :math:`C_1=1.0` and +:math:`C_2=0.1`, according to empirical tests. These estimates per source +region are then scaled by the additional prefactor :math:`F_{\mathrm{rank}}` +for the respective MPI rank, which is calculated from the ratio between measured +and estimated MPI rank load in the current batch. The measured load is +determined based on the transport sweep times, which are recorded by default for +diagnostics. + +Based on these load estimates, a load balancing routine tries to equalize the +work per MPI rank. To do so, the weights +:math:`\omega_{\mathrm{rank}}` in Equation :eq:`mpi_ownership` are +adjusted according to the deviation of the estimated rank load + +.. math:: + \mathrm{load}_{\mathrm{estimate}, \mathrm{rank}} = \frac{\sum\limits_{i \, + \in \, \mathrm{rank}} \mathrm{load}_{\mathrm{estimate},i}}{\sum\limits_{i=1} + ^{M} \mathrm{load}_{\mathrm{estimate}, i}} + +from the target load + +.. math:: + \mathrm{load}_{\mathrm{target}} = \frac{1}{N_{\mathrm{rank}}}\;\mathrm{,} + +where :math:`N_{\mathrm{rank}}` is the total number of MPI ranks. + +Changes to the weights :math:`\mathbf{\omega_{\mathrm{rank}}}` increase or +decrease the reach of a specific MPI rank, and thus the number of source regions +that belong to it. After each weight change, the rank load estimates are updated +according to the anticipated changes in the ownership of source regions, and the +new load estimates are then used again to calculate new weights. These load +optimization iterations continue until the estimated load imbalance is smaller +than 1% or until a maximum of 200 iterations is reached. + +After the load balancing, numerous source regions will belong to new MPI ranks. +The corresponding cell data is transferred to the new owner ranks and erased +from the previous owner ranks. Since both the iterative load optimization and +source region exchange can be computationally expensive, load balancing is +restricted to the first 5 simulation batches. Because random ray simulations +should use an appropriately large ray population, it is expected that sufficient +load estimate data has been recorded for the vast majority of cells after 5 +batches, and the load per MPI rank will not change significantly beyond +stochastic fluctuations associated with the changing quadrature. + --------------------------- Fundamental Sources of Bias --------------------------- @@ -1162,6 +1324,7 @@ in random ray particle transport are: .. _Cosgrove-2023: https://doi.org/10.1080/00295639.2023.2270618 .. _Ferrer-2016: https://doi.org/10.13182/NSE15-6 .. _Gunow-2018: https://dspace.mit.edu/handle/1721.1/119030 +.. _Balzer-2009: https://doi.org/10.1109/ISVD.2009.28 .. only:: html diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst index 2a0d301d4bc..b60fc5d9780 100644 --- a/docs/source/usersguide/install.rst +++ b/docs/source/usersguide/install.rst @@ -327,6 +327,9 @@ Prerequisites cmake -DOPENMC_USE_DAGMC=on -DCMAKE_PREFIX_PATH=/path/to/dagmc/installation .. + Distributed memory calculations with the random ray solver require MOAB + version 5.2.0 or later. + * MCPL_ library for reading and writing .mcpl files This option allows OpenMC to read and write MCPL (Monte Carlo Particle diff --git a/docs/source/usersguide/parallel.rst b/docs/source/usersguide/parallel.rst index ecbdd20b626..a43e99e8261 100644 --- a/docs/source/usersguide/parallel.rst +++ b/docs/source/usersguide/parallel.rst @@ -69,6 +69,9 @@ argument to :func:`openmc.run`:: openmc.run(mpi_args=['mpiexec', '-n', '32']) +Distributed memory calculations with the random ray solver for DAGMC geometries +require MOAB version 5.2.0 or later. + ---------------------- Maximizing Performance ---------------------- diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 1a27ad51602..60ea6197dda 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -93,6 +93,9 @@ class Region { //! Get a vector containing all the surfaces in the region expression vector surfaces() const; + //! Get size of surfaces + int n_surfaces() const { return expression_.size(); } + //---------------------------------------------------------------------------- // Accessors @@ -228,6 +231,9 @@ class Cell { //! Get a vector of surfaces in the cell virtual vector surfaces() const { return vector(); } + //! Get the number of surfaces in the cell + virtual int n_surfaces() const { return 0; } + //! Check if the cell region expression is simple virtual bool is_simple() const { return true; } @@ -420,6 +426,8 @@ class CSGCell : public Cell { // Methods vector surfaces() const override { return region_.surfaces(); } + int n_surfaces() const override { return region_.n_surfaces(); } + std::pair distance(Position r, Direction u, int32_t on_surface, GeometryState* p) const override { diff --git a/include/openmc/constants.h b/include/openmc/constants.h index a1d94e5819e..0770bcc5f0a 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -88,6 +88,15 @@ constexpr double MINIMUM_MACRO_XS {1e-6}; // window games are unbiased regardless of where the thresholds sit. constexpr double WEIGHT_WINDOW_REL_TOL {1e-9}; +// Maximum number of DAGMC entity handles to send when exchanging rays +// between MPI ranks. This caps the RayHistory length to avoid sending +// variable-length vectors. +constexpr int MAX_N_HANDLES {5}; + +// Maximum number of load optimization iterations to perform to balance +// the load between MPI ranks during random ray transport. +constexpr int ITER_LOAD_BALANCE {5}; + // ============================================================================ // MATH AND PHYSICAL CONSTANTS @@ -380,6 +389,7 @@ enum class SolverType { MONTE_CARLO, RANDOM_RAY }; enum class RandomRayVolumeEstimator { NAIVE, SIMULATION_AVERAGED, HYBRID }; enum class RandomRaySourceShape { FLAT, LINEAR, LINEAR_XY }; +enum class RandomRayGeomDim { TWO_DIM, THREE_DIM }; enum class RandomRaySampleMethod { PRNG, HALTON, S2 }; enum class RandomRaySolve { FORWARD, FORWARD_FOR_ADJOINT, ADJOINT }; diff --git a/include/openmc/random_ray/decomposition_map.h b/include/openmc/random_ray/decomposition_map.h new file mode 100644 index 00000000000..93e2fab1a7a --- /dev/null +++ b/include/openmc/random_ray/decomposition_map.h @@ -0,0 +1,123 @@ +#ifndef OPENMC_DECOMPOSITION_MAP_H +#define OPENMC_DECOMPOSITION_MAP_H + +#include "openmc/random_ray/flat_source_domain.h" +#include "openmc/random_ray/source_region.h" +#include "openmc/vector.h" + +namespace openmc { + +class DecompositionMap; + +namespace mpi { +extern DecompositionMap decomp_map; +} // namespace mpi + +class DecompositionMap { +public: + //---------------------------------------------------------------------------- + // Constructors + DecompositionMap(); + + //---------------------------------------------------------------------------- + // Methods + + // Methods to find rank centres that divide spatial domain up into equal + // Voronoi volumes + void initialize(); + void generate_rank_centers(); + void calculate_grid_points(int grid_points_total); + void initialize_voronoi_centers(); + void calculate_voronoi( + vector& position_sum_per_rank, vector& num_points_per_rank); + Position calculate_centroids( + const Position position_sum, const int num_points, int rank); + + // Methods to create and update subdomain list and exchange source region data + void exchange_sr_info( + ParallelMap& + discovered_source_regions); + bool any_discovered_source_regions( + ParallelMap& + discovered_source_regions); + void send_sr_data(int receiver, SourceRegion& sr_send); + void receive_sr_data(int sender, SourceRegion& sr_recv); + + // Methods for balancing the load between ranks + void balance_load(FlatSourceDomain* domain); + void update_load(FlatSourceDomain* domain, bool check_all_ranks, + vector& rank_load_combined, vector& load_ratio); + void redistribute_source_regions(FlatSourceDomain* domain); + + // Methods to find owner of source region + int find_owner(SourceRegionKey sr_key, Position r, + ParallelMap& + discovered_source_regions); + int find_closest_rank(Position r, bool test_all_ranks); + + // Method to calculate the load per rank based on the total number of hits in + // all source regions of a rank + void calculate_rank_load( + FlatSourceDomain* domain, double batch_transport_time); + double calculate_load_ratio(int rank); + + //---------------------------------------------------------------------------- + // Public data members + + // Map that relates a SourceRegionKey to the index of the MPI rank that + // contains that source region in its subdomain. + std::unordered_map + subdomain_map_; + + // Neighbors of each rank's Voronoi cell + std::unordered_set my_neighbors_; + + // Data to estimate rank loads + vector + num_base_source_region_RT_; // number of base source region ray trace + // operations per base source region + vector num_mesh_bin_RT_; // number of mesh bin ray trace operations + // per base source region + vector ray_tracing_cost_; + vector volume_base_sr_; + vector measured_rank_load_fractions_; + + // Load optimization + vector rank_weights_; + double target_load_; + +private: + //---------------------------------------------------------------------------- + // Private data members + SpatialBox* spatial_box_ = nullptr; + + // Voronoi cell calculation + vector grid_points_; + int grid_points_per_rank_ {125}; // default 5x5x5 grid points per rank + vector rank_centers_; // centers of each rank's Voronoi cell + + // Load calculation + vector estimated_rank_load_fractions_; + vector estimated_rank_load_totals_; + double estimated_load_sum_; + + // Coefficients for load calculation + double C1_ = 1.0; + double C2_ = 0.1; + double C3_ = 0.1; + + // Load optimization + double imbalance_tolerance_ = 0.01; // 1% imbalance tolerance + double optimization_history_factor_ = 1.0; + + // Miscellaneous + uint64_t n_base_sr_; + int negroups_; + double max_domain_length_; + bool is_linear_; + +}; // class DecompositionMap + +} // namespace openmc + +#endif // OPENMC_DECOMPOSITION_MAP_H diff --git a/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h index 09414fd4465..d257782687e 100644 --- a/include/openmc/random_ray/flat_source_domain.h +++ b/include/openmc/random_ray/flat_source_domain.h @@ -32,6 +32,7 @@ class FlatSourceDomain { void compute_k_eff(); virtual void normalize_scalar_flux_and_volumes( double total_active_distance_per_iteration); + bool is_geometry_3D(); int64_t add_source_to_scalar_flux(); virtual void batch_reset(); @@ -40,6 +41,7 @@ class FlatSourceDomain { void random_ray_tally(); virtual void accumulate_iteration_flux(); void output_to_vtk() const; + void output_to_vtk_decomp() const; void convert_external_sources(bool use_adjoint_sources); void count_external_source_regions(); void set_fw_adjoint_sources(); diff --git a/include/openmc/random_ray/parallel_map.h b/include/openmc/random_ray/parallel_map.h index 7f4f06d9996..cc77c787ad8 100644 --- a/include/openmc/random_ray/parallel_map.h +++ b/include/openmc/random_ray/parallel_map.h @@ -1,6 +1,7 @@ #ifndef OPENMC_RANDOM_RAY_PARALLEL_HASH_MAP_H #define OPENMC_RANDOM_RAY_PARALLEL_HASH_MAP_H +#include "openmc/error.h" #include "openmc/openmp_interface.h" #include @@ -142,7 +143,11 @@ class ParallelMap { ValueType& operator[](const KeyType& key) { Bucket& bucket = get_bucket(key); - return *bucket.map_[key].get(); + auto it = bucket.map_.find(key); + if (it == bucket.map_.end()) { + fatal_error("ParallelMap::operator[]: key not present in map."); + } + return *it->second; } ValueType* emplace(KeyType key, const ValueType& value) @@ -174,6 +179,23 @@ class ParallelMap { HashFunctor>::iterator()); } + // Return summed length of all buckets. + uint64_t size() + { + uint64_t total_size = 0; + for (const auto& bucket : buckets_) { + total_size += bucket.map_.size(); + } + return total_size; + } + + // Erase element by key. + void erase(const KeyType& key) + { + Bucket& bucket = get_bucket(key); + bucket.map_.erase(key); + } + private: //---------------------------------------------------------------------------- // Private Methods diff --git a/include/openmc/random_ray/random_ray.h b/include/openmc/random_ray/random_ray.h index b61d2d67aa8..5699ff63756 100644 --- a/include/openmc/random_ray/random_ray.h +++ b/include/openmc/random_ray/random_ray.h @@ -7,8 +7,79 @@ #include "openmc/random_ray/moment_matrix.h" #include "openmc/source.h" +#ifdef OPENMC_DAGMC_ENABLED +#include "DagMC.hpp" +#endif + namespace openmc { +// Container for MPI exchange +struct RayBufferContainer { + Position position; + Direction direction; + double distance_travelled; + vector angular_flux; + int surface; + bool is_active; + uint64_t ray_id; + int n_event; // Number of events (surface crossings) the ray has undergone + + // GeometryState scalar fields + int n_coord; + int cell_instance; + int n_coord_last; + int material; + int material_last; + double sqrtkT; + double sqrtkT_last; + + // GeometryState vector fields (sized to model::n_coord_levels at runtime) + // LocalCoord is POD, so we can send it as contiguous bytes + vector coord; + + // cell_last_ array + vector cell_last; + +#ifdef OPENMC_DAGMC_ENABLED + // DAGMC fields - fixed-size array to avoid variable-length vector + Direction last_dir; + moab::EntityHandle handles[MAX_N_HANDLES]; + int n_handles; // Actual number of valid handles (may be less than + // MAX_N_HANDLES) +#endif +}; + +// Container for MPI exchange +struct RayExchangeData { + Position position; + Direction direction; + double distance_travelled; + int surface; + bool is_active; + uint64_t ray_id; + int n_event; // Number of events (surface crossings) the ray has undergone + + // GeometryState scalar fields + int n_coord; + int cell_instance; + int n_coord_last; + int material; + int material_last; + double sqrtkT; + double sqrtkT_last; + +#ifdef OPENMC_DAGMC_ENABLED + // DAGMC fields - fixed-size array to avoid variable-length vector + Direction last_dir; + moab::EntityHandle handles[MAX_N_HANDLES]; + int n_handles; // Actual number of valid handles (may be less than + // MAX_N_HANDLES) +#endif +}; + +// Forward declare +class FlatSourceDomain; + /* * The RandomRay class encompasses data and methods for transporting random rays * through the model. It is a small extension of the Particle class. @@ -38,24 +109,32 @@ class RandomRay : public Particle { SourceRegionHandle& srh, double distance, bool is_active, Position r); void initialize_ray(uint64_t ray_id, FlatSourceDomain* domain); + void restart_ray(FlatSourceDomain* domain, RayExchangeData& data, + float* angular_flux, LocalCoord* coord, int* cell_last_data); uint64_t transport_history_based_single_ray(); SourceSite sample_prng(); SourceSite sample_halton(); SourceSite sample_s2(); + bool has_left_subdomain(); + void pack_ray_for_buffer(double distance_buffer, Position position_buffer); + //---------------------------------------------------------------------------- // Static data members static double distance_inactive_; // Inactive (dead zone) ray length static double distance_active_; // Active ray length static unique_ptr ray_source_; // Starting source for ray sampling static RandomRaySourceShape source_shape_; // Flag for linear source + static RandomRayGeomDim geom_dim_; // Flag for 2D vs 3D geometry static RandomRaySampleMethod sample_method_; // Flag for sampling method //---------------------------------------------------------------------------- // Public data members vector angular_flux_; + RayBufferContainer exchange_data_; bool ray_trace_only_ {false}; // If true, only perform geometry operations + int owner_rank_ {C_NONE}; // Rank that owns this ray based on its position private: //---------------------------------------------------------------------------- @@ -72,6 +151,7 @@ class RandomRay : public Particle { double distance_travelled_ {0}; bool is_active_ {false}; bool is_alive_ {true}; + bool is_local_ {true}; }; // class RandomRay } // namespace openmc diff --git a/include/openmc/random_ray/random_ray_simulation.h b/include/openmc/random_ray/random_ray_simulation.h index e186c549f92..c338003edef 100644 --- a/include/openmc/random_ray/random_ray_simulation.h +++ b/include/openmc/random_ray/random_ray_simulation.h @@ -3,6 +3,7 @@ #include "openmc/random_ray/flat_source_domain.h" #include "openmc/random_ray/linear_source_domain.h" +#include "openmc/random_ray/ray_bank.h" namespace openmc { @@ -24,12 +25,15 @@ class RandomRaySimulation { void prepare_local_fixed_sources_adjoint(); void prepare_adjoint_simulation(bool from_forward); void simulate(); - void output_simulation_results() const; + void output_simulation_results(); void instability_check( int64_t n_hits, double k_eff, double& avg_miss_rate) const; void print_results_random_ray(uint64_t total_geometric_intersections, double avg_miss_rate, int negroups, int64_t n_source_regions, - int64_t n_external_source_regions) const; + int64_t n_external_source_regions, uint64_t avg_num_comms, + double max_load_imbalance) const; + void transport_sweep(); + void transport_sweep_decomp(RayBank& RB); //---------------------------------------------------------------------------- // Accessors @@ -52,6 +56,13 @@ class RandomRaySimulation { // Number of energy groups int negroups_; + // Average number of ray communications between rank per batch + uint64_t avg_num_communication_rounds_ {0}; + + // Tracks whether geometry-dependent one-time setup has already run for + // this simulation object across forward/adjoint solves. + bool geometry_setup_complete_ {false}; + }; // class RandomRaySimulation //============================================================================ diff --git a/include/openmc/random_ray/ray_bank.h b/include/openmc/random_ray/ray_bank.h new file mode 100644 index 00000000000..9606adde147 --- /dev/null +++ b/include/openmc/random_ray/ray_bank.h @@ -0,0 +1,67 @@ +#ifndef OPENMC_RAY_BANK_H +#define OPENMC_RAY_BANK_H + +#include "openmc/random_ray/flat_source_domain.h" +#include "openmc/random_ray/random_ray.h" +#include "openmc/random_ray/source_region.h" +#include "openmc/vector.h" +#ifdef OPENMC_MPI +#include +#endif + +namespace openmc { + +class RayBank { +public: + //---------------------------------------------------------------------------- + // Constructors + RayBank(); + + //---------------------------------------------------------------------------- + // Methods + void buffer_ray_data_to_send(RandomRay& ray, FlatSourceDomain* domain); + void update(FlatSourceDomain* domain); + int ray_bank_size(); + void communicate_rays(); + void communicate_message_metadata(); + void update_my_ray_list(FlatSourceDomain* domain); + bool is_any_ray_alive(); + + //---------------------------------------------------------------------------- + // Public data members + vector my_ray_list_; + +private: + //---------------------------------------------------------------------------- + // Private data members + int total_receiving_rays_; + int negroups_; + + // Per-rank send buffers for ray data, including geometry state and angular + // flux + struct RankSendBuffers { + vector ray_data; + vector angular_flux; + vector coord; + vector cell_last; + int count = 0; // Number of rays buffered for this rank + }; + std::unordered_map ray_send_buffer_; + int reserved_buffer_size_ = 32; // Initial reserved size for send buffers, can + // be tuned based on expected ray counts + + // Vector that contains the number of rays to be received from each rank + vector num_messages_receiving_; + vector num_messages_sending_; + + // vectors that received ray data + vector received_ray_data_; + vector received_angular_flux_data_; + vector received_coord_; + vector received_cell_last_; + +}; // class RayBank + +} // namespace openmc + +#endif // OPENMC_RAY_BANK_H diff --git a/include/openmc/random_ray/source_region.h b/include/openmc/random_ray/source_region.h index 1d2bbe1e8dc..62abe5377aa 100644 --- a/include/openmc/random_ray/source_region.h +++ b/include/openmc/random_ray/source_region.h @@ -165,6 +165,7 @@ class SourceRegionHandle { Position* centroid_t_; MomentMatrix* mom_matrix_; MomentMatrix* mom_matrix_t_; + SourceRegionKey* key_; // 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. @@ -254,6 +255,9 @@ class SourceRegionHandle { MomentMatrix& mom_matrix_t() { return *mom_matrix_t_; } const MomentMatrix mom_matrix_t() const { return *mom_matrix_t_; } + SourceRegionKey& key() { return *key_; } + const SourceRegionKey key() const { return *key_; } + std::unordered_set& volume_task() { return *volume_task_; @@ -311,24 +315,13 @@ class SourceRegionHandle { }; // class SourceRegionHandle -class SourceRegion { +class ScalarSourceRegionFields { public: - //---------------------------------------------------------------------------- - // Constructors - SourceRegion(int negroups, bool is_linear); - SourceRegion() = default; - - //---------------------------------------------------------------------------- - // Public Data members - - //--------------------------------------- - // Scalar fields int material_ {0}; //!< Index in openmc::model::materials array int temperature_idx_ { 0}; //!< Index into the MGXS array representing temperature double density_mult_ {1.0}; //!< A density multiplier queried from the cell //!< corresponding to the source region. - OpenMPMutex lock_; double volume_ { 0.0}; //!< Volume (computed from the sum of ray crossing lengths) double volume_t_ {0.0}; //!< Volume totaled over all iterations @@ -343,6 +336,7 @@ class SourceRegion { // Mesh that subdivides this source region int mesh_ {C_NONE}; //!< Index in openmc::model::meshes array that subdivides //!< this source region + SourceRegionKey key_ {0, 0}; //!< The key (base source region + mesh bin) int64_t parent_sr_ {C_NONE}; //!< Index of a parent source region Position position_ { 0.0, 0.0, 0.0}; //!< A position somewhere inside the region @@ -355,6 +349,30 @@ class SourceRegion { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; //!< The spatial moment matrix MomentMatrix mom_matrix_t_ {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; //!< The spatial moment matrix accumulated over all iterations +}; + +class SourceRegion { +public: + //---------------------------------------------------------------------------- + // Constructors + SourceRegion(int negroups, bool is_linear); + SourceRegion(const SourceRegionHandle& handle); + SourceRegion() = default; + + //---------------------------------------------------------------------------- + // Methods + void merge(SourceRegion& sr_add, bool is_linear); + + //---------------------------------------------------------------------------- + // Public Data members + + //--------------------------------------- + // Scalar fields + + OpenMPMutex lock_; + + // Container with all scalar fields of a source region + ScalarSourceRegionFields scalars_; // A set of volume tally tasks. This more complicated data structure is // convenient for ensuring that volumes are only tallied once per source @@ -474,6 +492,9 @@ class SourceRegionContainer { return mom_matrix_t_[sr]; } + SourceRegionKey& key(int64_t sr) { return key_[sr]; } + const SourceRegionKey key(int64_t sr) const { return key_[sr]; } + MomentArray& source_gradients(int64_t sr, int g) { return source_gradients_[index(sr, g)]; @@ -662,6 +683,7 @@ class SourceRegionContainer { vector centroid_t_; vector mom_matrix_; vector mom_matrix_t_; + vector key_; // 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. diff --git a/include/openmc/source.h b/include/openmc/source.h index 51b54a1d106..63e2e244f93 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -92,6 +92,8 @@ class Source { static unique_ptr create(pugi::xml_node node); + bool satisfies_spatial_constraints(Position r) const; + protected: // Strategy used for rejecting sites when constraints are applied. KILL means // that sites are always accepted but if they don't satisfy constraints, they @@ -104,7 +106,6 @@ class Source { // Methods for constraints void read_constraints(pugi::xml_node node); - bool satisfies_spatial_constraints(Position r) const; bool satisfies_energy_constraints(double E) const; bool satisfies_time_constraints(double time) const; diff --git a/include/openmc/timer.h b/include/openmc/timer.h index d928aad4560..86e7ffc0b6e 100644 --- a/include/openmc/timer.h +++ b/include/openmc/timer.h @@ -32,6 +32,13 @@ extern Timer time_event_surface_crossing; extern Timer time_event_collision; extern Timer time_event_death; extern Timer time_update_src; +extern Timer time_ray_comms; +extern Timer time_decomposition_handling; +extern Timer time_load_balance; +extern Timer time_ray_buffering; +extern Timer time_generate_voronoi_centers; +extern Timer time_source_region_exchange; +extern Timer time_mpi_imbalance; } // namespace simulation diff --git a/src/random_lcg.cpp b/src/random_lcg.cpp index f2a81fc1c7e..5bdedc59083 100644 --- a/src/random_lcg.cpp +++ b/src/random_lcg.cpp @@ -56,7 +56,7 @@ SkipAheadCoefficients future_seed_coefficients(uint64_t n) //============================================================================== // 64 bit implementation of the PCG-RXS-M-XS 64-bit state / 64-bit output -// geneator Adapted from: https://github.com/imneme/pcg-c, in particular +// generator Adapted from: https://github.com/imneme/pcg-c, in particular // https://github.com/imneme/pcg-c/blob/83252d9c23df9c82ecb42210afed61a7b42402d7/include/pcg_variants.h#L188-L192 // @techreport{oneill:pcg2014, // title = "PCG: A Family of Simple Fast Space-Efficient Statistically Good diff --git a/src/random_ray/decomposition_map.cpp b/src/random_ray/decomposition_map.cpp new file mode 100644 index 00000000000..7e208940e48 --- /dev/null +++ b/src/random_ray/decomposition_map.cpp @@ -0,0 +1,1147 @@ +#include "openmc/random_ray/decomposition_map.h" + +#include "openmc/cell.h" +#include "openmc/constants.h" +#include "openmc/message_passing.h" +#include "openmc/mgxs_interface.h" +#include "openmc/random_lcg.h" +#include "openmc/random_ray/flat_source_domain.h" +#include "openmc/random_ray/random_ray.h" +#include "openmc/simulation.h" +#include "openmc/timer.h" +#include "openmc/vector.h" +#include + +#ifdef OPENMC_MPI + +namespace openmc { + +namespace mpi { +DecompositionMap decomp_map; +} + +// Constructor +DecompositionMap::DecompositionMap() {} + +void DecompositionMap::initialize() +{ + negroups_ = data::mg.num_energy_groups_; + estimated_rank_load_fractions_.resize(mpi::n_procs, 0.0); + measured_rank_load_fractions_.resize(mpi::n_procs, 0.0); + estimated_rank_load_totals_.resize(mpi::n_procs, 0.0); + target_load_ = 1.0 / mpi::n_procs; + rank_weights_.resize(mpi::n_procs, 1.0); + + spatial_box_ = dynamic_cast( + dynamic_cast(RandomRay::ray_source_.get())->space()); + + double x_length = + spatial_box_->upper_right().x - spatial_box_->lower_left().x; + double y_length = + spatial_box_->upper_right().y - spatial_box_->lower_left().y; + double z_length = + spatial_box_->upper_right().z - spatial_box_->lower_left().z; + + max_domain_length_ = + sqrt(x_length * x_length + y_length * y_length + z_length * z_length); + + is_linear_ = RandomRay::source_shape_ != RandomRaySourceShape::FLAT; + + // Count the number of source regions in the model + n_base_sr_ = 0; + for (const auto& c : model::cells) { + if (c->type_ == Fill::MATERIAL) { + n_base_sr_ += c->n_instances(); + } + } + + num_base_source_region_RT_.resize(n_base_sr_, 0); + num_mesh_bin_RT_.resize(n_base_sr_, 0); + ray_tracing_cost_.resize(n_base_sr_); + volume_base_sr_.resize(n_base_sr_, 0.0); +} + +void DecompositionMap::generate_rank_centers() +{ + + // Calculate grid points that are used for Voronoi cells + int grid_points_total = grid_points_per_rank_ * mpi::n_procs; + if (mpi::master) + printf("Calculating %d grid points for Voronoi tessellation...\n", + grid_points_total); + calculate_grid_points(grid_points_total); + + // Initialize points with random positions + initialize_voronoi_centers(); + + double err = 1.0; + double precision = 1e-3; // corresponding to 0.001 cm position change + int it = 0; + int max_iterations = 100; + + // Lloyd's algorithm to move Voronoi centers to centroids of their cells + // https://en.wikipedia.org/wiki/Lloyd%27s_algorithm + while (err > precision && it < max_iterations) { + // Reset error to determine maximum movement below + err = 0.0; + + vector position_sum_per_rank(mpi::n_procs, Position(0, 0, 0)); + vector num_points_per_rank(mpi::n_procs, 0); + + // Compute Voronoi cells by summing up all position values of mesh grid + // points that are closest to a Voronoi center + calculate_voronoi(position_sum_per_rank, num_points_per_rank); + + for (int rank = 0; rank < mpi::n_procs; rank++) { + + // Calculate centroid of the cell + Position centroid = calculate_centroids( + position_sum_per_rank[rank], num_points_per_rank[rank], rank); + + // Calculate movement for convergence check + double movement = (centroid - rank_centers_[rank]).norm(); + + // Move rank center to centroid + rank_centers_[rank] = centroid; + + // Record maximum movement + if (movement > err) { + err = movement; + } + } + + it++; + } + + if (mpi::master) { + if (it == max_iterations) { + warning("Lloyd's algorithm did not converge within the maximum number of " + "iterations."); + } else { + printf("Lloyd's algorithm converged in %d iterations.\n", it); + } + } +} + +// Calculate grid points needed for calculating Voronoi cells +void DecompositionMap::calculate_grid_points(int grid_points_total) +{ + + // Calculate length along each dimension + vector domain_length(3); + domain_length[0] = + spatial_box_->upper_right().x - spatial_box_->lower_left().x; + domain_length[1] = + spatial_box_->upper_right().y - spatial_box_->lower_left().y; + domain_length[2] = + spatial_box_->upper_right().z - spatial_box_->lower_left().z; + + double volume = domain_length[0] * domain_length[1] * domain_length[2]; + + // For each dimension, determine grid points along that direction based on + // aspect ratio: domain_length / volume^(1/3) = grid_points_dimension / + // grid_points_total^(1/3). Check if any dimension is so distorted that it + // would only receive minimum of 1 grid point and flag that direction to + // correct total number of grid points. + int excluded_dimension = -1; + vector grid_points_per_dimension(3); + for (int i = 0; i < 3; i++) { + double grid_points_estimate = + ((domain_length[i] / cbrt(volume)) * cbrt(grid_points_total)); + + if (grid_points_estimate > 1) { + grid_points_per_dimension[i] = round(grid_points_estimate); + } else { + excluded_dimension = i; + grid_points_per_dimension[i] = 1; + } + } + + // If problem is 2D, exclude z direction + if (RandomRay::geom_dim_ == RandomRayGeomDim::TWO_DIM) { + excluded_dimension = 2; + grid_points_per_dimension[2] = 1; + } + + // If one dimension is excluded, recalculate grid points in other two + // dimensions based on area. + if (excluded_dimension != -1) { + double area = 1.0; + + for (int i = 0; i < 3; i++) { + if (i == excluded_dimension) + continue; + area *= domain_length[i]; + } + + for (int i = 0; i < 3; i++) { + if (i == excluded_dimension) + continue; + grid_points_per_dimension[i] = + round((domain_length[i] / sqrt(area)) * sqrt(grid_points_total)); + } + } + + // Adjust grid points in each dimension to match actual total number of grid + // points to the number of grid points requested + double new_total = grid_points_per_dimension[0] * + grid_points_per_dimension[1] * + grid_points_per_dimension[2]; + double adjustment; + if (excluded_dimension != -1) { + // When one dimension is excluded, use square root for the other two + // dimensions + adjustment = sqrt(grid_points_total / new_total); + } else { + // When all dimensions are used, use cubic root + adjustment = cbrt(grid_points_total / new_total); + } + + // Multiply adjustment factor + for (int i = 0; i < 3; i++) { + if (i == excluded_dimension) + continue; + grid_points_per_dimension[i] = + round(grid_points_per_dimension[i] * adjustment); + } + + // Calculate spacing between grid points in each dimension + vector delta_value(3, 0.0); + + for (int i = 0; i < 3; i++) { + if (grid_points_per_dimension[i] > 1) { + delta_value[i] = + (domain_length[i] - 2 * TINY_BIT) / (grid_points_per_dimension[i] - 1); + } + } + + // Initialize point at center of domain (in case of only 1 grid point in a + // dimension) + double x = spatial_box_->lower_left().x + domain_length[0] * 0.5; + double y = spatial_box_->lower_left().y + domain_length[1] * 0.5; + double z = spatial_box_->lower_left().z + domain_length[2] * 0.5; + + // Generate all grid points + for (int i = 0; i < grid_points_per_dimension[0]; i++) { + if (grid_points_per_dimension[0] > 1) { + x = spatial_box_->lower_left().x + TINY_BIT + i * delta_value[0]; + } + for (int j = 0; j < grid_points_per_dimension[1]; j++) { + if (grid_points_per_dimension[1] > 1) { + y = spatial_box_->lower_left().y + TINY_BIT + j * delta_value[1]; + } + for (int k = 0; k < grid_points_per_dimension[2]; k++) { + if (grid_points_per_dimension[2] > 1) { + z = spatial_box_->lower_left().z + TINY_BIT + k * delta_value[2]; + } + // Add grid point + grid_points_.push_back({x, y, z}); + } + } + } + + // Check if mesh grid points are inside spatial domain, which can be different + // from a box. If not, erase them. + for (int i = grid_points_.size() - 1; i >= 0; i--) { + Position xi = grid_points_[i]; + + bool is_inside_domain = + RandomRay::ray_source_->satisfies_spatial_constraints(xi); + + if (!is_inside_domain) { + grid_points_.erase(grid_points_.begin() + i); + } + } + + if (mpi::master && grid_points_.size() < grid_points_total) { + warning(fmt::format("Spatial constraints reduced grid points for Voronoi " + "tesselation from {} to {}.", + grid_points_total, grid_points_.size())); + } +} + +// Places random points in the spatial domain. +// Each point corresponds to the initial center of a rank. +void DecompositionMap::initialize_voronoi_centers() +{ + rank_centers_.resize(mpi::n_procs); + + uint64_t seed = openmc_get_seed(); + int rank_cnt = 0; + + // Sample random positions to start with + while (rank_cnt < mpi::n_procs) { + + double x = prn(&seed); + double y = prn(&seed); + double z = 0.0; + if (RandomRay::geom_dim_ == RandomRayGeomDim::THREE_DIM) { + z = prn(&seed); + } else { + z = 0.5; // Mid-plane in z direction + } + + Position xi {x, y, z}; + + // Make a small shift in position to avoid geometry floating point issues at + // boundaries + Position shift {FP_COINCIDENT, FP_COINCIDENT, FP_COINCIDENT}; + xi = (spatial_box_->lower_left() + shift) + + xi * ((spatial_box_->upper_right() - shift) - + (spatial_box_->lower_left() + shift)); + + bool is_inside_domain = + RandomRay::ray_source_->satisfies_spatial_constraints(xi); + + if (is_inside_domain) { + rank_centers_[rank_cnt] = xi; + rank_cnt++; + } + } +} + +// Determine the distance of each mesh grid point to all rank centers. +// Sum up the positions of all mesh grid points that are closest to a given rank +// center for computation of centroid later. Record number of grid points per +// rank center. +void DecompositionMap::calculate_voronoi( + vector& position_sum_per_rank, vector& num_points_per_rank) +{ + +// Assign each point to the closest rank center +#pragma omp parallel for schedule(static) + for (int p = 0; p < grid_points_.size(); p++) { + Position point = grid_points_[p]; + int closest_rank = C_NONE; + double min_distance = INFTY; + + // Find closest rank center + for (int rank = 0; rank < mpi::n_procs; rank++) { + double dist = (point - rank_centers_[rank]).norm(); + // Power Voronoi diagram uses squared distances + dist = dist * dist - rank_weights_[rank]; + + if (dist < min_distance) { + min_distance = dist; + closest_rank = rank; + } + } + + if (mpi::master && closest_rank == C_NONE) { + fatal_error("Could not find closest rank for Voronoi cell point " + + std::to_string(p) + "."); + } + +// Accumulate point coordinates for the closest rank +#pragma omp atomic + position_sum_per_rank[closest_rank].x += point.x; +#pragma omp atomic + position_sum_per_rank[closest_rank].y += point.y; +#pragma omp atomic + position_sum_per_rank[closest_rank].z += point.z; + +// Record number of mesh grid points for closest rank +#pragma omp atomic + num_points_per_rank[closest_rank]++; + } +} + +Position DecompositionMap::calculate_centroids( + const Position position_sum, const int num_points, int rank) +{ + + // check if any points have been recorded in rank + if (num_points == 0) { + fatal_error("Rank " + std::to_string(rank) + + " has no Voronoi cell points. This indicates that the number " + "of grid points for the Voronoi tesselation is too coarse. " + "Requires source code modificaiton to fix."); + } + + Position centroid = position_sum; + + // Divide by number of points + double n = static_cast(num_points); + centroid.x /= n; + centroid.y /= n; + centroid.z /= n; + + return centroid; +} + +bool DecompositionMap::any_discovered_source_regions( + ParallelMap& + discovered_source_regions) +{ + + simulation::time_decomposition_handling.start(); + + int flag = 0; + if (discovered_source_regions.begin() != discovered_source_regions.end()) { + flag = 1; + } + + MPI_Allreduce(MPI_IN_PLACE, &flag, 1, MPI_INT, MPI_MAX, mpi::intracomm); + + return flag > 0; + + simulation::time_decomposition_handling.start(); +} + +void DecompositionMap::exchange_sr_info( + ParallelMap& + discovered_source_regions) +{ + + // Communicate maps + for (int rank = 0; rank < mpi::n_procs; rank++) { + + // Send size + uint64_t bcast_size = 0; + if (rank == mpi::rank) { + for (const auto& pair : discovered_source_regions) { + // Only broadcast source regions that have non-zero volume, i.e. regions + // discovered during active phase of ray + if (pair.second.scalars_.volume_ > 0.0) { + bcast_size++; + } + } + } + + MPI_Bcast(&bcast_size, 1, MPI_UINT64_T, rank, mpi::intracomm); + + if (bcast_size > 0) { + + vector local_base_ids(bcast_size); + vector local_mesh_bins(bcast_size); + + if (rank == mpi::rank) { + // fill in vectors to be sent + int i = 0; + for (const auto& pair : discovered_source_regions) { + if (pair.second.scalars_.volume_ > 0.0) { + SourceRegionKey sr_key = pair.first; + local_base_ids[i] = sr_key.base_source_region_id; + local_mesh_bins[i] = sr_key.mesh_bin; + i++; + } + } + } + + // Broadcast all data + MPI_Bcast( + local_base_ids.data(), bcast_size, MPI_INT64_T, rank, mpi::intracomm); + MPI_Bcast( + local_mesh_bins.data(), bcast_size, MPI_INT64_T, rank, mpi::intracomm); + + // Update subdomain map + for (int j = 0; j < bcast_size; j++) { + + SourceRegionKey sr_key(local_base_ids[j], local_mesh_bins[j]); + + // check if already in map or not, i.e. has someone else discovered that + // region already? + if (subdomain_map_.find(sr_key) == subdomain_map_.end()) { + subdomain_map_[sr_key] = rank; + } else { + int resident_rank = subdomain_map_[sr_key]; // current owner + int challenger_rank = rank; // current broadcasting rank + + double ratio_resident = calculate_load_ratio(resident_rank); + double ratio_challenger = calculate_load_ratio(challenger_rank); + double resident_rank_load = + ratio_resident * estimated_rank_load_totals_[resident_rank]; + double challenger_rank_load = + ratio_challenger * estimated_rank_load_totals_[challenger_rank]; + + // If load of challenger rank is lower, assign source region to that + // rank, otherwise resident keeps it + int sender; + int receiver; + if (challenger_rank_load < resident_rank_load) { + subdomain_map_[sr_key] = challenger_rank; + sender = resident_rank; + receiver = challenger_rank; + } else { + // resident keeps it + sender = challenger_rank; + receiver = resident_rank; + } + + // Broadcast load of exchanged source region such that each rank can + // update load balance + double bcast_load = 0; + if (mpi::rank == sender) { + SourceRegion& contested_sr = discovered_source_regions[sr_key]; + double volume_sr = contested_sr.scalars_.volume_; + bcast_load = + C1_ * + (contested_sr.scalars_.n_hits_ / simulation::current_batch) * + negroups_ + + volume_sr * ray_tracing_cost_[sr_key.base_source_region_id]; + } + MPI_Bcast(&bcast_load, 1, MPI_DOUBLE, sender, mpi::intracomm); + + double load_change_fraction = + bcast_load / estimated_load_sum_; // load fraction update + estimated_rank_load_fractions_[sender] -= load_change_fraction; + estimated_rank_load_fractions_[receiver] += load_change_fraction; + double load_change_total = bcast_load; // load total update + estimated_rank_load_totals_[sender] -= load_change_total; + estimated_rank_load_totals_[receiver] += load_change_total; + + // Communicate source region data and merge on receiver side + if (mpi::rank == sender) { + SourceRegion& sr = discovered_source_regions[sr_key]; + send_sr_data(receiver, sr); + + // clear old source region data from discovered regions map + discovered_source_regions.erase(sr_key); + } + if (mpi::rank == receiver) { + SourceRegion& sr = discovered_source_regions[sr_key]; + SourceRegion sr_recv(negroups_, is_linear_); + receive_sr_data(sender, sr_recv); + + sr.merge(sr_recv, is_linear_); + } + } + } + } + } +} + +void DecompositionMap::send_sr_data(int receiver, SourceRegion& sr_send) +{ + + int num_scalar_messages = 1; + //! NOTE: update if new vector fields are added to SourceExchangeVectors + //! struct + int num_vector_messages = 4; + if (is_linear_) { + num_vector_messages += 4; + } + if (settings::run_mode == RunMode::FIXED_SOURCE) { + num_vector_messages += 1; + } + int num_requests = num_scalar_messages + num_vector_messages; + + vector requests(num_requests); + int req_idx = 0; + + // Send scalar data to receiver + MPI_Isend(&sr_send.scalars_, sizeof(ScalarSourceRegionFields), MPI_BYTE, + receiver, 1, mpi::intracomm, &requests[req_idx]); + req_idx++; + + // Send vector data to receiver + // Tags hardcoded to avoid confusion if new fields are not added sequentially + MPI_Isend(sr_send.scalar_flux_old_.data(), negroups_, MPI_DOUBLE, receiver, 2, + mpi::intracomm, &requests[req_idx]); + req_idx++; + + MPI_Isend(sr_send.scalar_flux_new_.data(), negroups_, MPI_DOUBLE, receiver, 3, + mpi::intracomm, &requests[req_idx]); + req_idx++; + + MPI_Isend(sr_send.source_.data(), negroups_, MPI_FLOAT, receiver, 4, + mpi::intracomm, &requests[req_idx]); + req_idx++; + + if (settings::run_mode == RunMode::FIXED_SOURCE) { + MPI_Isend(sr_send.external_source_.data(), negroups_, MPI_FLOAT, receiver, + 5, mpi::intracomm, &requests[req_idx]); + req_idx++; + } + + MPI_Isend(sr_send.scalar_flux_final_.data(), negroups_, MPI_DOUBLE, receiver, + 6, mpi::intracomm, &requests[req_idx]); + req_idx++; + + if (is_linear_) { + MPI_Isend(sr_send.source_gradients_.data(), 3 * negroups_, MPI_DOUBLE, + receiver, 7, mpi::intracomm, &requests[req_idx]); + req_idx++; + + MPI_Isend(sr_send.flux_moments_old_.data(), 3 * negroups_, MPI_DOUBLE, + receiver, 8, mpi::intracomm, &requests[req_idx]); + req_idx++; + + MPI_Isend(sr_send.flux_moments_new_.data(), 3 * negroups_, MPI_DOUBLE, + receiver, 9, mpi::intracomm, &requests[req_idx]); + req_idx++; + + MPI_Isend(sr_send.flux_moments_t_.data(), 3 * negroups_, MPI_DOUBLE, + receiver, 10, mpi::intracomm, &requests[req_idx]); + req_idx++; + } + + if (req_idx != num_requests) { + fatal_error(fmt::format( + "Number of MPI requests does not match number of messages sent." + "Check if num_vector_messages corresponds to number of transferred " + "vectors.")); + } + + // Wait for all communication to complete + MPI_Waitall(num_requests, requests.data(), MPI_STATUSES_IGNORE); +} + +void DecompositionMap::receive_sr_data(int sender, SourceRegion& sr_recv) +{ + + // Receive scalar data from sender + MPI_Recv(&sr_recv.scalars_, sizeof(ScalarSourceRegionFields), MPI_BYTE, + sender, 1, mpi::intracomm, MPI_STATUS_IGNORE); + + // Receive vector data from sender + MPI_Recv(sr_recv.scalar_flux_old_.data(), negroups_, MPI_DOUBLE, sender, 2, + mpi::intracomm, MPI_STATUS_IGNORE); + MPI_Recv(sr_recv.scalar_flux_new_.data(), negroups_, MPI_DOUBLE, sender, 3, + mpi::intracomm, MPI_STATUS_IGNORE); + MPI_Recv(sr_recv.source_.data(), negroups_, MPI_FLOAT, sender, 4, + mpi::intracomm, MPI_STATUS_IGNORE); + + if (settings::run_mode == RunMode::FIXED_SOURCE) { + MPI_Recv(sr_recv.external_source_.data(), negroups_, MPI_FLOAT, sender, 5, + mpi::intracomm, MPI_STATUS_IGNORE); + } + + MPI_Recv(sr_recv.scalar_flux_final_.data(), negroups_, MPI_DOUBLE, sender, 6, + mpi::intracomm, MPI_STATUS_IGNORE); + + if (is_linear_) { + MPI_Recv(sr_recv.source_gradients_.data(), 3 * negroups_, MPI_DOUBLE, + sender, 7, mpi::intracomm, MPI_STATUS_IGNORE); + MPI_Recv(sr_recv.flux_moments_old_.data(), 3 * negroups_, MPI_DOUBLE, + sender, 8, mpi::intracomm, MPI_STATUS_IGNORE); + MPI_Recv(sr_recv.flux_moments_new_.data(), 3 * negroups_, MPI_DOUBLE, + sender, 9, mpi::intracomm, MPI_STATUS_IGNORE); + MPI_Recv(sr_recv.flux_moments_t_.data(), 3 * negroups_, MPI_DOUBLE, sender, + 10, mpi::intracomm, MPI_STATUS_IGNORE); + } +} + +int DecompositionMap::find_owner(SourceRegionKey sr_key, Position r, + ParallelMap& + discovered_source_regions) +{ + + // Check if source region key is in subdomain map + auto it = subdomain_map_.find(sr_key); + if (it != subdomain_map_.end()) { + return it->second; + } + + // Check if already recorded in newly discovered source regions + discovered_source_regions.lock(sr_key); + bool sr_key_discovered = discovered_source_regions.contains(sr_key); + discovered_source_regions.unlock(sr_key); + if (sr_key_discovered) { + return mpi::rank; + } + + // If not found in either map, check which rank owns source + // region beased on location + int closest_rank = find_closest_rank(r, true); + return closest_rank; +} + +int DecompositionMap::find_closest_rank(Position r, bool test_all_ranks) +{ + + int closest_rank = C_NONE; + double min_distance = INFTY; + vector test_ranks; + + if (test_all_ranks) { + test_ranks.resize(mpi::n_procs); + std::iota(test_ranks.begin(), test_ranks.end(), + 0); // fill with 0, 1, ..., n_procs-1 + } else { + // convert unordered set of neighboring ranks to vector and add self rank + test_ranks = vector(mpi::decomp_map.my_neighbors_.begin(), + mpi::decomp_map.my_neighbors_.end()); + test_ranks.push_back(mpi::rank); + } + + // Find closest rank center + for (int rank : test_ranks) { + double dist = (r - rank_centers_[rank]).norm(); + // Distance function corresponding to weighted power Voronoi diagram + dist = dist * dist - rank_weights_[rank]; + if (dist < min_distance) { + min_distance = dist; + closest_rank = rank; + } + } + + if (mpi::master && closest_rank == C_NONE) { + fatal_error( + "Could not find closest rank for new source region at position (" + + std::to_string(r.x) + ", " + std::to_string(r.y) + ", " + + std::to_string(r.z) + ")."); + } + + return closest_rank; +} + +void DecompositionMap::calculate_rank_load( + FlatSourceDomain* domain, double batch_transport_time) +{ + + // Reset local volumes of base source regions, which might change when source + // regions change rank ownership + std::fill(volume_base_sr_.begin(), volume_base_sr_.end(), 0.0); + + // Add volumes of newly discovered source regions + vector mesh_bins_per_base_sr_local(n_base_sr_, 0); + for (const auto& [sr_key, sr] : domain->discovered_source_regions_) { + volume_base_sr_[sr_key.base_source_region_id] += sr.scalars_.volume_; + } + + // Add volumes of known source regions + for (int64_t sr = 0; sr < domain->n_source_regions(); sr++) { + SourceRegionKey sr_key = domain->source_regions_.key(sr); + uint64_t base_sr = sr_key.base_source_region_id; + volume_base_sr_[base_sr] += domain->source_regions_.volume_t(sr); + volume_base_sr_[base_sr] += domain->source_regions_.volume(sr); + } + +// Calculate ray tracing cost per base source region +#pragma omp parallel for + for (uint64_t bsr = 0; bsr < n_base_sr_; bsr++) { + if (volume_base_sr_[bsr] > 0.0) { + ray_tracing_cost_[bsr] = + (C2_ * static_cast(num_base_source_region_RT_[bsr]) + + C3_ * static_cast(num_mesh_bin_RT_[bsr])) / + volume_base_sr_[bsr]; + } else { + ray_tracing_cost_[bsr] = 0.0; + } + } + + // Accumulate load of known source regions + double local_estimated_load = 0.0; +#pragma omp parallel for reduction(+ : local_estimated_load) + for (int64_t sr = 0; sr < domain->n_source_regions(); sr++) { + SourceRegionKey sr_key = domain->source_regions_.key(sr); + uint64_t base_sr = sr_key.base_source_region_id; + + // Calculate volume of unique source region to weight volume-dependent ray + // tracing cost + double volume_sr = + domain->source_regions_.volume_t(sr) + domain->source_regions_.volume(sr); + + // Calculate load of source region + double load_sr = + C1_ * (domain->source_regions_.n_hits(sr) / simulation::current_batch) * + negroups_ + + volume_sr * ray_tracing_cost_[base_sr]; + + // Accumulate to local estimated load + local_estimated_load += load_sr; + } + + // Accumulate load of newly discovered source regions + for (const auto& [sr_key, sr] : domain->discovered_source_regions_) { + uint64_t base_sr = sr_key.base_source_region_id; + double volume_sr = sr.scalars_.volume_; + double load_sr = + C1_ * (sr.scalars_.n_hits_ / simulation::current_batch) * negroups_ + + volume_sr * ray_tracing_cost_[base_sr]; + local_estimated_load += load_sr; + } + + // Communicate estimated load across ranks + MPI_Allgather(&local_estimated_load, 1, MPI_DOUBLE, + estimated_rank_load_totals_.data(), 1, MPI_DOUBLE, mpi::intracomm); + estimated_load_sum_ = std::accumulate(estimated_rank_load_totals_.begin(), + estimated_rank_load_totals_.end(), 0.0); + + // Communicate measured load across ranks + MPI_Allgather(&batch_transport_time, 1, MPI_DOUBLE, + measured_rank_load_fractions_.data(), 1, MPI_DOUBLE, mpi::intracomm); + double measured_load_sum = + std::accumulate(measured_rank_load_fractions_.begin(), + measured_rank_load_fractions_.end(), 0.0); + + // Calculate fractions + for (int rank = 0; rank < mpi::n_procs; rank++) { + estimated_rank_load_fractions_[rank] = + estimated_rank_load_totals_[rank] / estimated_load_sum_; + measured_rank_load_fractions_[rank] = + measured_rank_load_fractions_[rank] / measured_load_sum; + } + + // Reset ray trace counters + fill(num_base_source_region_RT_.begin(), num_base_source_region_RT_.end(), 0); + fill(num_mesh_bin_RT_.begin(), num_mesh_bin_RT_.end(), 0); +} + +void DecompositionMap::balance_load(FlatSourceDomain* domain) +{ + + // Optimization parameters + int max_iterations = 200; + int it_outer = 0; + double adaptation_factor = 1; + double min_adaptation_factor = 0.01; + double max_adaptation_factor = 2; + double avg_rank_distance = + max_domain_length_ / + cbrt(mpi::n_procs); // rough estimate of average distance between ranks + double weight_scale = + avg_rank_distance * avg_rank_distance * optimization_history_factor_; + double beta = 0.6; // momentum damping + bool check_all_ranks = true; + + vector weight_change(mpi::n_procs, 0.0); + vector combined_rank_load(mpi::n_procs, 0.0); + vector load_ratio(mpi::n_procs, 0.0); + + // Combine estimated load with measured load ratios + for (int rank = 0; rank < mpi::n_procs; rank++) { + load_ratio[rank] = calculate_load_ratio(rank); + combined_rank_load[rank] = + load_ratio[rank] * estimated_rank_load_totals_[rank]; + } + + double combined_load_sum = + std::accumulate(combined_rank_load.begin(), combined_rank_load.end(), 0.0); + + for (int rank = 0; rank < mpi::n_procs; rank++) { + combined_rank_load[rank] = (combined_rank_load[rank] / combined_load_sum); + } + double max_load = + *std::max_element(combined_rank_load.begin(), combined_rank_load.end()); + double max_imbalance = (max_load - target_load_) / target_load_; + double prev_imbalance = max_imbalance; + + // History tracking + vector imbalance_history; + vector> weight_history; + imbalance_history.push_back(max_imbalance); + weight_history.push_back(rank_weights_); + + // Change weights to equalize load based on combined load estimates + while (max_imbalance > imbalance_tolerance_ && it_outer < max_iterations) { + + it_outer++; + + for (int rank = 0; rank < mpi::n_procs; rank++) { + double corr = ((combined_rank_load[rank] - target_load_) / target_load_) * + weight_scale; + weight_change[rank] = + beta * weight_change[rank] + + (1.0 - beta) * corr; // keep some inertia from previous changes to + // prevent oscillations + double damping = std::clamp(max_imbalance / prev_imbalance, 0.1, + 1.0); // dampening factor based on whether we are getting closer to + // convergence or not, prevents big jumps, if too big a change, + // more dampening is applied + rank_weights_[rank] -= adaptation_factor * damping * weight_change[rank]; + } + + if (simulation::current_batch > 1) { + // Check distances to all ranks only in first batch, otherwise only + // recorded neighbors + check_all_ranks = false; + } + + // Calculate new load after weight update + update_load(domain, check_all_ranks, combined_rank_load, load_ratio); + double max_load = + *std::max_element(combined_rank_load.begin(), combined_rank_load.end()); + max_imbalance = (max_load - target_load_) / target_load_; + + // Store imbalance history + imbalance_history.push_back(max_imbalance); + weight_history.push_back(rank_weights_); + + // Adaptive factor + if (max_imbalance > prev_imbalance) + adaptation_factor = + std::max(adaptation_factor * 0.5, min_adaptation_factor); + else + adaptation_factor = + std::min(adaptation_factor * 1.05, max_adaptation_factor); + + prev_imbalance = max_imbalance; + } + + // Check convergence and adjust history optimization factor depending on + // failure mode to enable better convergence in the following batch + if (it_outer == max_iterations) { + if (mpi::master) { + warning("MPI load balancing has not converged after " + + std::to_string(max_iterations) + " iterations."); + } + + // Check if oscillating or simply slow convergence + int direction_changes = 0; + + for (int i = 1; i < imbalance_history.size(); i++) { + // Calculate change + double change = imbalance_history[i] - imbalance_history[i - 1]; + + // Count direction changes + if (i > 1) { + double prev_change = + imbalance_history[i - 1] - imbalance_history[i - 2]; + if (change * prev_change < 0) { + direction_changes++; + } + } + } + + // Check for oscillations + double oscillation_ratio = + (double)direction_changes / (imbalance_history.size() - 2); + + if (oscillation_ratio > 0.4) { + // decrease weight for next batch if oscillating + optimization_history_factor_ = + std::max(optimization_history_factor_ * 0.5, min_adaptation_factor); + } else { + // increase weight for faster convergence if too slow + optimization_history_factor_ = + std::min(optimization_history_factor_ * 1.2, max_adaptation_factor); + } + + // Check which iteration had the best imbalance and revert to those weights + auto min_it = + std::min_element(imbalance_history.begin(), imbalance_history.end()); + int best_index = std::distance(imbalance_history.begin(), min_it); + double best_imbalance = *min_it; + rank_weights_ = weight_history[best_index]; + + if (mpi::master) { + printf( + "Best imbalance during optimization was %.2f%% at iteration %d. \n", + best_imbalance * 100.0, best_index); + } + + if (best_index == 0) { + // if no improvement at all, just keep current decomposition and return + // without redistributing + return; + } + } else { + optimization_history_factor_ = 1.0; // reset history factor if converged + if (mpi::master) { + printf("MPI load balancing converged after %d iterations. Max. " + "imbalance: %.2f%% \n", + it_outer, max_imbalance * 100.0); + } + } + + // Redistribute source regions according to new weights determined in + // optimization + redistribute_source_regions(domain); +} + +void DecompositionMap::update_load(FlatSourceDomain* domain, + bool check_all_ranks, vector& combined_rank_load, + vector& load_ratio) +{ + + vector load(mpi::n_procs, 0); + +// Add up source region hits to respective rank depending of position of +// centroid of each source region +#pragma omp parallel + { + // number of hits per thread + vector thread_load(mpi::n_procs, 0); + +#pragma omp for + for (int64_t sr = 0; sr < domain->n_source_regions(); sr++) { + Position centroid = domain->source_regions_.centroid(sr); + int owner = find_closest_rank(centroid, check_all_ranks); + double volume_sr = domain->source_regions_.volume_t(sr); + thread_load[owner] += + load_ratio[owner] * + (C1_ * + (domain->source_regions_.n_hits(sr) / simulation::current_batch) * + negroups_ + + volume_sr * ray_tracing_cost_[domain->source_regions_.key(sr) + .base_source_region_id]); + } + +// Combine results from different threads +#pragma omp critical(combining_loads) + { + for (int i = 0; i < mpi::n_procs; i++) { + load[i] += thread_load[i]; + } + } + } + + // Communicate new load estimates across ranks + MPI_Allreduce(MPI_IN_PLACE, load.data(), mpi::n_procs, MPI_DOUBLE, MPI_SUM, + mpi::intracomm); + double load_sum = std::accumulate(load.begin(), load.end(), 0.0); + + // Update new combined rank load fractions + for (int rank = 0; rank < mpi::n_procs; rank++) { + combined_rank_load[rank] = load[rank] / load_sum; + } +} + +void DecompositionMap::redistribute_source_regions(FlatSourceDomain* domain) +{ + + // Map of source regions to be sent to other ranks + std::unordered_map> sr_send_list; + // Number of source regions to be received from other ranks + vector num_sr_receiving(mpi::n_procs, 0); + + // Local source region container that contains updated list + SourceRegionContainer source_regions_new(negroups_, is_linear_); + + // Each rank identifies source regions that need to be transferred to new + // owner and updates subdomain map accordingly + for (int64_t sr = 0; sr < domain->n_source_regions(); sr++) { + Position centroid = domain->source_regions_.centroid(sr); + int owner = find_closest_rank(centroid, true); + + // If owner changed, write source region to list of outbound source regions + if (owner != mpi::rank) { + sr_send_list[owner].push_back(sr); + } + // If owner did not change, add source region to new local container + else { + source_regions_new.push_back( + domain->source_regions_.get_source_region_handle(sr)); + } + } + + // Each rank informs other ranks about ownership changes to update subdomain + // map + for (int rank = 0; rank < mpi::n_procs; rank++) { + + // Send size + int bcast_size = 0; + if (rank == mpi::rank) { + for (const auto& pair : sr_send_list) { + bcast_size += pair.second.size(); + } + } + MPI_Bcast(&bcast_size, 1, MPI_INT, rank, mpi::intracomm); + + if (bcast_size > 0) { + vector rank_ids(bcast_size); + vector base_ids(bcast_size); + vector mesh_bins(bcast_size); + + // Current owner prepares communication and fills vectors to be sent + if (rank == mpi::rank) { + int i = 0; + for (auto& pair : sr_send_list) { + + int receiver = pair.first; // new owner + vector& sr_indices = + pair.second; // vector of source region indices + + // Iterate through all source regions for this key + for (int sr_idx : sr_indices) { + SourceRegionKey sr_key = domain->source_regions_.key(sr_idx); + rank_ids[i] = receiver; + base_ids[i] = sr_key.base_source_region_id; + mesh_bins[i] = sr_key.mesh_bin; + i++; + } + } + } + + // Broadcast source region data + MPI_Bcast(rank_ids.data(), bcast_size, MPI_INT, rank, mpi::intracomm); + MPI_Bcast(base_ids.data(), bcast_size, MPI_INT64_T, rank, mpi::intracomm); + MPI_Bcast( + mesh_bins.data(), bcast_size, MPI_INT64_T, rank, mpi::intracomm); + + // Every rank updates subdomain map + for (int j = 0; j < bcast_size; j++) { + + // Re-assemble source region key and extract new owner + SourceRegionKey sr_key(base_ids[j], mesh_bins[j]); + int rank_new = rank_ids[j]; + + // Update subdomain_map with new owner + subdomain_map_[sr_key] = rank_new; + + // If calling rank is new owner, increase count of source regions coming + // from sending rank (current broadcaster) + if (mpi::rank == rank_new) { + num_sr_receiving[rank] += 1; + } + } + } + } + + // Clear source_region_map_ + domain->source_region_map_.clear(); + + // Send source region data to new owner + for (auto& pair : sr_send_list) { + int receiver = pair.first; // destination rank + vector& sr_indices = pair.second; // vector of source region indices + + // Iterate through all source regions for this key + for (int sr_idx : sr_indices) { + SourceRegionKey sr_key = domain->source_regions_.key(sr_idx); + SourceRegionHandle srh = + domain->source_regions_.get_source_region_handle(sr_idx); + SourceRegion sr(srh); + send_sr_data(receiver, sr); + } + } + + // Record starting source region ID for tally reinitialization later + int64_t start_sr_id = source_regions_new.n_source_regions(); + + // Receive source regions + for (int sender = 0; sender < mpi::n_procs; sender++) { + + if (sender == mpi::rank) { + if (num_sr_receiving[sender] > 0) { + fatal_error("Rank sends source regions to itself. Rank should not " + "receive source regions from itself."); + } + continue; // skip self + } + + int num_sr = num_sr_receiving[sender]; + for (int i = 0; i < num_sr; ++i) { + SourceRegion sr_recv(negroups_, is_linear_); + receive_sr_data(sender, sr_recv); + source_regions_new.push_back(sr_recv); + } + } + + // Update source regions in domain to new container + domain->source_regions_ = source_regions_new; + + // Update source region map + for (int64_t sr = 0; sr < domain->n_source_regions(); sr++) { + SourceRegionKey key = domain->source_regions_.key(sr); + domain->source_region_map_[key] = sr; + } + + // Reinitialise tallies + domain->convert_source_regions_to_tallies(start_sr_id); +} + +double DecompositionMap::calculate_load_ratio(int rank) +{ + if (estimated_rank_load_fractions_[rank] > 0.0) { + return measured_rank_load_fractions_[rank] / + estimated_rank_load_fractions_[rank]; + } else { + return 1.0; + } +} + +} // namespace openmc +#endif // OPENMC_MPI diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 83128fdaa05..f3ef5fec78d 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -9,6 +9,7 @@ #include "openmc/mgxs_interface.h" #include "openmc/output.h" #include "openmc/plot.h" +#include "openmc/random_ray/decomposition_map.h" #include "openmc/random_ray/random_ray.h" #include "openmc/simulation.h" #include "openmc/tallies/filter.h" @@ -84,6 +85,7 @@ void FlatSourceDomain::batch_reset() // 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_.centroid_iteration(sr) = {0.0, 0.0, 0.0}; source_regions_.volume(sr) = 0.0; source_regions_.volume_sq(sr) = 0.0; } @@ -183,6 +185,7 @@ void FlatSourceDomain::normalize_scalar_flux_and_volumes( // update the simulation-averaged cell-wise volume estimates #pragma omp parallel for for (int64_t sr = 0; sr < n_source_regions(); sr++) { + source_regions_.centroid_t(sr) += source_regions_.centroid_iteration(sr); source_regions_.volume_t(sr) += source_regions_.volume(sr); source_regions_.volume_sq_t(sr) += source_regions_.volume_sq(sr); source_regions_.volume_naive(sr) = @@ -191,6 +194,11 @@ void FlatSourceDomain::normalize_scalar_flux_and_volumes( source_regions_.volume_sq_t(sr) / source_regions_.volume_t(sr); 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; + } } } @@ -365,6 +373,18 @@ void FlatSourceDomain::compute_k_eff() p[sr] = sr_fission_source_new; } + // Sum up fission rates across all ranks +#ifdef OPENMC_MPI + if (mpi::n_procs > 1) { + simulation::time_decomposition_handling.start(); + MPI_Allreduce( + MPI_IN_PLACE, &fission_rate_old, 1, MPI_DOUBLE, MPI_SUM, mpi::intracomm); + MPI_Allreduce( + MPI_IN_PLACE, &fission_rate_new, 1, MPI_DOUBLE, MPI_SUM, mpi::intracomm); + simulation::time_decomposition_handling.stop(); + } +#endif + double k_eff_new = k_eff_ * (fission_rate_new / fission_rate_old); double H = 0.0; @@ -382,6 +402,18 @@ void FlatSourceDomain::compute_k_eff() } } +#ifdef OPENMC_MPI + if (mpi::n_procs > 1) { + simulation::time_decomposition_handling.start(); + if (mpi::master) { + MPI_Reduce(MPI_IN_PLACE, &H, 1, MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm); + } else { + MPI_Reduce(&H, nullptr, 1, MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm); + } + simulation::time_decomposition_handling.stop(); + } +#endif + // Adds entropy value to shared entropy vector in openmc namespace. simulation::entropy.push_back(H); @@ -585,6 +617,13 @@ double FlatSourceDomain::compute_fixed_source_normalization_factor() const } } +#ifdef OPENMC_MPI + if (mpi::n_procs > 1) { + MPI_Allreduce(MPI_IN_PLACE, &simulation_external_source_strength, 1, + MPI_DOUBLE, MPI_SUM, mpi::intracomm); + } +#endif + // Step 2 is to determine the total user-specified external source strength double user_external_source_strength = 0.0; for (auto& ext_source : model::external_sources) { @@ -722,6 +761,15 @@ void FlatSourceDomain::random_ray_tally() // see what index that score corresponds to. If that score is a flux score, // then we divide it by volume. if (volume_normalized_flux_tallies_) { +#ifdef OPENMC_MPI + if (mpi::n_procs > 1) { + for (auto& volumes : tally_volumes_) { + MPI_Allreduce(MPI_IN_PLACE, volumes.data(), + static_cast(volumes.size()), MPI_DOUBLE, MPI_SUM, + mpi::intracomm); + } + } +#endif for (int i = 0; i < model::tallies.size(); i++) { Tally& tally {*model::tallies[i]}; #pragma omp parallel for @@ -815,6 +863,7 @@ void FlatSourceDomain::output_to_vtk() const // Relate voxel spatial locations to random ray source regions vector voxel_indices(Nx * Ny * Nz); + vector voxel_indices_key(Nx * Ny * Nz); vector voxel_positions(Nx * Ny * Nz); vector weight_windows(Nx * Ny * Nz); float min_weight = 1e20; @@ -835,6 +884,7 @@ void FlatSourceDomain::output_to_vtk() const bool found = exhaustive_find_cell(p); if (!found) { + voxel_indices_key[z * Ny * Nx + y * Nx + x] = {-1, -1}; voxel_indices[z * Ny * Nx + y * Nx + x] = -1; voxel_positions[z * Ny * Nx + y * Nx + x] = sample; weight_windows[z * Ny * Nx + y * Nx + x] = 0.0; @@ -848,6 +898,7 @@ void FlatSourceDomain::output_to_vtk() const sr = it->second; } + voxel_indices_key[z * Ny * Nx + y * Nx + x] = sr_key; voxel_indices[z * Ny * Nx + y * Nx + x] = sr; voxel_positions[z * Ny * Nx + y * Nx + x] = sample; @@ -889,7 +940,6 @@ void FlatSourceDomain::output_to_vtk() const std::fprintf(plot, "LOOKUP_TABLE default\n"); for (int i = 0; i < Nx * Ny * Nz; i++) { int64_t fsr = voxel_indices[i]; - int64_t source_element = fsr * negroups_ + g; float flux = 0; if (fsr >= 0) { flux = evaluate_flux_at_point(voxel_positions[i], fsr, g); @@ -952,7 +1002,6 @@ void FlatSourceDomain::output_to_vtk() const int temp = source_regions_.temperature_idx(fsr); if (mat != MATERIAL_VOID) { 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); double sigma_f = sigma_f_[(mat * ntemperature_ + temp) * negroups_ + g] * @@ -1006,6 +1055,443 @@ void FlatSourceDomain::output_to_vtk() const } } +// Variant of output_to_vtk() for domain decomposition case. Every rank +// contributes the data of its own subdomain, reduced onto the master rank for +// file output. +#ifdef OPENMC_MPI +void FlatSourceDomain::output_to_vtk_decomp() const +{ + + if (mpi::master) { + // Rename .h5 plot filename(s) to .vtk filenames + for (int p = 0; p < model::plots.size(); p++) { + PlottableInterface* plot = model::plots[p].get(); + plot->path_plot() = + plot->path_plot().substr(0, plot->path_plot().find_last_of('.')) + + ".vtk"; + } + + // Print header information + print_plot(); + } + + // Outer loop over plots + for (int p = 0; p < model::plots.size(); p++) { + + // Get handle to OpenMC plot object and extract params + Plot* openmc_plot = dynamic_cast(model::plots[p].get()); + + // Random ray plots only support voxel plots + if (!openmc_plot) { + warning(fmt::format("Plot {} is invalid plot type -- only voxel plotting " + "is allowed in random ray mode.", + p)); + continue; + } else if (openmc_plot->type_ != Plot::PlotType::voxel) { + warning(fmt::format("Plot {} is invalid plot type -- only voxel plotting " + "is allowed in random ray mode.", + p)); + continue; + } + + int Nx = openmc_plot->pixels_[0]; + int Ny = openmc_plot->pixels_[1]; + int Nz = openmc_plot->pixels_[2]; + Position origin = openmc_plot->origin_; + Position width = openmc_plot->width_; + Position ll = origin - width / 2.0; + double x_delta = width.x / Nx; + double y_delta = width.y / Ny; + double z_delta = width.z / Nz; + std::string filename = openmc_plot->path_plot(); + + // Tag plots written during the forward solve of an adjoint run + if (solve_ == RandomRaySolve::FORWARD_FOR_ADJOINT) { + auto dot = filename.find_last_of('.'); + filename = filename.substr(0, dot) + ".forward" + filename.substr(dot); + } + + // Perform sanity checks on file size + uint64_t bytes = Nx * Ny * Nz * (negroups_ + 1 + 1 + 1) * sizeof(float); + write_message(5, "Processing plot {}: {}... (Estimated size is {} MB)", + openmc_plot->id(), filename, bytes / 1.0e6); + if (bytes / 1.0e9 > 1.0) { + if (mpi::master) { + warning( + "Voxel plot specification is very large (>1 GB). Plotting may be " + "slow."); + } + } else if (bytes / 1.0e9 > 100.0) { + if (mpi::master) { + fatal_error( + "Voxel plot specification is too large (>100 GB). Exiting."); + } + } + + // Relate voxel spatial locations to random ray source regions + vector voxel_indices(Nx * Ny * Nz); + vector voxel_positions(Nx * Ny * Nz); + vector weight_windows(Nx * Ny * Nz); + vector my_voxel_ids; + float min_weight = 1e20; +#pragma omp parallel for collapse(3) reduction(min : min_weight) + for (int z = 0; z < Nz; z++) { + for (int y = 0; y < Ny; y++) { + for (int x = 0; x < Nx; x++) { + Position sample; + sample.z = ll.z + z_delta / 2.0 + z * z_delta; + sample.y = ll.y + y_delta / 2.0 + y * y_delta; + sample.x = ll.x + x_delta / 2.0 + x * x_delta; + Particle p; + p.r() = sample; + p.r_last() = sample; + p.E() = 1.0; + p.E_last() = 1.0; + p.u() = {1.0, 0.0, 0.0}; + + bool found = exhaustive_find_cell(p); + if (!found) { + voxel_indices[z * Ny * Nx + y * Nx + x] = -1; + voxel_positions[z * Ny * Nx + y * Nx + x] = sample; + weight_windows[z * Ny * Nx + y * Nx + x] = 0.0; + continue; + } + + SourceRegionKey sr_key = lookup_source_region_key(p); + int64_t sr = -1; + auto it_sr = source_region_map_.find(sr_key); + if (it_sr != source_region_map_.end()) { + sr = it_sr->second; + } + + voxel_indices[z * Ny * Nx + y * Nx + x] = sr; + voxel_positions[z * Ny * Nx + y * Nx + x] = sample; + + // Assumed master rank = 0 + int assigned_rank = 0; + // Which rank is responsible + auto it = mpi::decomp_map.subdomain_map_.find(sr_key); + if (it != mpi::decomp_map.subdomain_map_.end()) { + assigned_rank = it->second; + } + if (assigned_rank == mpi::rank) { +#pragma omp critical(create_my_voxel_ids) + { + my_voxel_ids.push_back(z * Ny * Nx + y * Nx + x); + } + } + + if (variance_reduction::weight_windows.size() == 1) { + auto [ww_found, ww] = + variance_reduction::weight_windows[0]->get_weight_window(p); + float weight = ww.lower_weight; + weight_windows[z * Ny * Nx + y * Nx + x] = weight; + if (weight < min_weight) + min_weight = weight; + } + } + } + } + + double source_normalization_factor = + compute_fixed_source_normalization_factor(); + + // Open file for writing + std::FILE* plot = nullptr; + if (mpi::master) { + plot = std::fopen(filename.c_str(), "wb"); + + // Write vtk metadata + std::fprintf(plot, "# vtk DataFile Version 2.0\n"); + std::fprintf(plot, "Dataset File\n"); + std::fprintf(plot, "BINARY\n"); + std::fprintf(plot, "DATASET STRUCTURED_POINTS\n"); + std::fprintf(plot, "DIMENSIONS %d %d %d\n", Nx, Ny, Nz); + std::fprintf(plot, "ORIGIN %lf %lf %lf\n", ll.x, ll.y, ll.z); + std::fprintf(plot, "SPACING %lf %lf %lf\n", x_delta, y_delta, z_delta); + std::fprintf(plot, "POINT_DATA %d\n", Nx * Ny * Nz); + } + + int vector_size = Nx * Ny * Nz; + vector vector_out_float(vector_size, 0.0); + vector vector_out_int(vector_size, 0); + + int64_t num_neg = 0; + int64_t num_samples = 0; + float min_flux = 0.0; + float max_flux = -1.0e20; + + // Plot multigroup flux data + for (int g = 0; g < negroups_; g++) { + + for (int voxel_id : my_voxel_ids) { + int64_t fsr = voxel_indices[voxel_id]; + float flux = 0; + if (fsr >= 0) { + flux = evaluate_flux_at_point(voxel_positions[voxel_id], fsr, g); + if (flux < 0.0) + flux = FlatSourceDomain::evaluate_flux_at_point( + voxel_positions[voxel_id], fsr, g); + } + if (flux < 0.0) { + num_neg++; + if (flux < min_flux) { + min_flux = flux; + } + } + if (flux > max_flux) + max_flux = flux; + num_samples++; + vector_out_float[voxel_id] = flux; + } + + if (mpi::master) { + MPI_Reduce(MPI_IN_PLACE, vector_out_float.data(), vector_size, + MPI_FLOAT, MPI_SUM, 0, mpi::intracomm); + MPI_Reduce( + MPI_IN_PLACE, &num_neg, 1, MPI_INT64_T, MPI_SUM, 0, mpi::intracomm); + MPI_Reduce(MPI_IN_PLACE, &num_samples, 1, MPI_INT64_T, MPI_SUM, 0, + mpi::intracomm); + } else { + MPI_Reduce(vector_out_float.data(), nullptr, vector_size, MPI_FLOAT, + MPI_SUM, 0, mpi::intracomm); + MPI_Reduce( + &num_neg, nullptr, 1, MPI_INT64_T, MPI_SUM, 0, mpi::intracomm); + MPI_Reduce( + &num_samples, nullptr, 1, MPI_INT64_T, MPI_SUM, 0, mpi::intracomm); + } + + if (mpi::master) { + std::fprintf(plot, "SCALARS flux_group_%d float\n", g); + std::fprintf(plot, "LOOKUP_TABLE default\n"); + for (float value : vector_out_float) { + float print_value = convert_to_big_endian(value); + std::fwrite(&print_value, sizeof(float), 1, plot); + } + } + + fill(vector_out_float.begin(), vector_out_float.end(), 0.0); + } + + // Slightly negative fluxes can be normal when sampling corners of linear + // source regions. However, very common and high magnitude negative fluxes + // may indicate numerical instability. + if (mpi::master && num_neg > 0) { + warning(fmt::format("{} plot samples ({:.4f}%) contained negative fluxes " + "(minumum found = {:.2e} maximum_found = {:.2e})", + num_neg, (100.0 * num_neg) / num_samples, min_flux, max_flux)); + } + + // Plot FSRs + for (int voxel_id : my_voxel_ids) { + int fsr = voxel_indices[voxel_id]; + float value = future_prn(10, fsr); + vector_out_float[voxel_id] = value; + } + + if (mpi::master) { + MPI_Reduce(MPI_IN_PLACE, vector_out_float.data(), vector_size, MPI_FLOAT, + MPI_SUM, 0, mpi::intracomm); + } else { + MPI_Reduce(vector_out_float.data(), nullptr, vector_size, MPI_FLOAT, + MPI_SUM, 0, mpi::intracomm); + } + + if (mpi::master) { + std::fprintf(plot, "SCALARS FSRs float\n"); + std::fprintf(plot, "LOOKUP_TABLE default\n"); + for (float value : vector_out_float) { + float print_value = convert_to_big_endian(value); + std::fwrite(&print_value, sizeof(float), 1, plot); + } + } + + fill(vector_out_float.begin(), vector_out_float.end(), 0.0); + + // Plot Materials + for (int voxel_id : my_voxel_ids) { + int mat = -1; + int fsr = voxel_indices[voxel_id]; + if (fsr >= 0) { + mat = source_regions_.material(fsr); + } + vector_out_int[voxel_id] = mat + 1; // To avoid -1 for void (MPI_SUM) + } + + if (mpi::master) { + MPI_Reduce(MPI_IN_PLACE, vector_out_int.data(), vector_size, MPI_INT, + MPI_SUM, 0, mpi::intracomm); + } else { + MPI_Reduce(vector_out_int.data(), nullptr, vector_size, MPI_INT, MPI_SUM, + 0, mpi::intracomm); + } + + if (mpi::master) { + std::fprintf(plot, "SCALARS Materials int\n"); + std::fprintf(plot, "LOOKUP_TABLE default\n"); + + for (int value : vector_out_int) { + int print_value = convert_to_big_endian(value); + std::fwrite(&print_value, sizeof(int), 1, plot); + } + } + + fill(vector_out_int.begin(), vector_out_int.end(), 0); + + // Plot rank subdomains + for (int voxel_id : my_voxel_ids) { + int rank_id = mpi::rank; + float value = future_prn(10, rank_id); + vector_out_float[voxel_id] = value; + } + + if (mpi::master) { + MPI_Reduce(MPI_IN_PLACE, vector_out_float.data(), vector_size, MPI_FLOAT, + MPI_SUM, 0, mpi::intracomm); + } else { + MPI_Reduce(vector_out_float.data(), nullptr, vector_size, MPI_FLOAT, + MPI_SUM, 0, mpi::intracomm); + } + + if (mpi::master) { + std::fprintf(plot, "SCALARS rank_subdomains float\n"); + std::fprintf(plot, "LOOKUP_TABLE default\n"); + + for (float value : vector_out_float) { + float print_value = convert_to_big_endian(value); + std::fwrite(&print_value, sizeof(float), 1, plot); + } + } + + fill(vector_out_float.begin(), vector_out_float.end(), 0.0); + + // Plot measured load based on transport sweep timers + for (int voxel_id : my_voxel_ids) { + float value = mpi::decomp_map.measured_rank_load_fractions_[mpi::rank]; + vector_out_float[voxel_id] = value; + } + + if (mpi::master) { + MPI_Reduce(MPI_IN_PLACE, vector_out_float.data(), vector_size, MPI_FLOAT, + MPI_SUM, 0, mpi::intracomm); + } else { + MPI_Reduce(vector_out_float.data(), nullptr, vector_size, MPI_FLOAT, + MPI_SUM, 0, mpi::intracomm); + } + + if (mpi::master) { + std::fprintf(plot, "SCALARS measured_load float\n"); + std::fprintf(plot, "LOOKUP_TABLE default\n"); + + for (float value : vector_out_float) { + float print_value = convert_to_big_endian(value); + std::fwrite(&print_value, sizeof(float), 1, plot); + } + } + + fill(vector_out_float.begin(), vector_out_float.end(), 0.0); + + // Plot fission source + if (settings::run_mode == RunMode::EIGENVALUE) { + for (int voxel_id : my_voxel_ids) { + int64_t fsr = voxel_indices[voxel_id]; + float total_fission = 0.0; + if (fsr >= 0) { + int mat = source_regions_.material(fsr); + int temp = source_regions_.temperature_idx(fsr); + if (mat != MATERIAL_VOID) { + for (int g = 0; g < negroups_; g++) { + float flux = + evaluate_flux_at_point(voxel_positions[voxel_id], fsr, g); + double sigma_f = + sigma_f_[(mat * ntemperature_ + temp) * negroups_ + g] * + source_regions_.density_mult(fsr); + total_fission += sigma_f * flux; + } + } + } + vector_out_float[voxel_id] = total_fission; + } + } else { + for (int voxel_id : my_voxel_ids) { + int64_t fsr = voxel_indices[voxel_id]; + int mat = source_regions_.material(fsr); + int temp = source_regions_.temperature_idx(fsr); + float total_external = 0.0f; + if (fsr >= 0) { + for (int g = 0; g < negroups_; g++) { + // External sources are already divided by sigma_t, so we need to + // multiply it back to get the true external source. + double sigma_t = 1.0; + if (mat != MATERIAL_VOID) { + sigma_t = sigma_t_[(mat * ntemperature_ + temp) * negroups_ + g] * + source_regions_.density_mult(fsr); + } + total_external += source_regions_.external_source(fsr, g) * sigma_t; + } + } + vector_out_float[voxel_id] = total_external; + } + } + + if (mpi::master) { + MPI_Reduce(MPI_IN_PLACE, vector_out_float.data(), vector_size, MPI_FLOAT, + MPI_SUM, 0, mpi::intracomm); + } else { + MPI_Reduce(vector_out_float.data(), nullptr, vector_size, MPI_FLOAT, + MPI_SUM, 0, mpi::intracomm); + } + + if (mpi::master) { + if (settings::run_mode == RunMode::EIGENVALUE) { + std::fprintf(plot, "SCALARS total_fission_source float\n"); + } else { + std::fprintf(plot, "SCALARS external_source float\n"); + } + std::fprintf(plot, "LOOKUP_TABLE default\n"); + for (float value : vector_out_float) { + float print_value = convert_to_big_endian(value); + std::fwrite(&print_value, sizeof(float), 1, plot); + } + } + + fill(vector_out_float.begin(), vector_out_float.end(), 0.0); + + // Plot weight window data + if (variance_reduction::weight_windows.size() == 1) { + for (int voxel_id : my_voxel_ids) { + float weight = weight_windows[voxel_id]; + if (weight == 0.0) + weight = min_weight; + vector_out_float[voxel_id] = weight; + } + + if (mpi::master) { + MPI_Reduce(MPI_IN_PLACE, vector_out_float.data(), vector_size, + MPI_FLOAT, MPI_SUM, 0, mpi::intracomm); + } else { + MPI_Reduce(vector_out_float.data(), nullptr, vector_size, MPI_FLOAT, + MPI_SUM, 0, mpi::intracomm); + } + + if (mpi::master) { + std::fprintf(plot, "SCALARS weight_window_lower float\n"); + std::fprintf(plot, "LOOKUP_TABLE default\n"); + + for (float value : vector_out_float) { + float print_value = convert_to_big_endian(value); + std::fwrite(&print_value, sizeof(float), 1, plot); + } + } + } + + if (mpi::master) { + std::fclose(plot); + } + } +} +#endif // OPENMC_MPI + void FlatSourceDomain::apply_external_source_to_source_region( int src_idx, SourceRegionHandle& srh) { @@ -1303,7 +1789,6 @@ void FlatSourceDomain::set_fw_adjoint_sources() source_regions_.external_source_present(sr) = 0; } } - // Divide the fixed source term by sigma t (to save time when applying each // iteration) #pragma omp parallel for @@ -1630,6 +2115,7 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( SourceRegion* sr_ptr = discovered_source_regions_.emplace(sr_key, {negroups_, is_linear}); SourceRegionHandle handle {*sr_ptr}; + handle.key() = sr_key; // Determine the material int gs_i_cell = gs.lowest_coord().cell(); @@ -1655,7 +2141,6 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( handle.material() = material; handle.temperature_idx() = temp; - handle.density_mult() = cell.density_mult(gs.cell_instance()); // Store the mesh index (if any) assigned to this source region @@ -1717,7 +2202,7 @@ void FlatSourceDomain::finalize_discovered_source_regions() // Extract keys for entries with a valid volume. vector keys; for (const auto& pair : discovered_source_regions_) { - if (pair.second.volume_ > 0.0) { + if (pair.second.scalars_.volume_ > 0.0) { keys.push_back(pair.first); } } @@ -1849,4 +2334,59 @@ int64_t FlatSourceDomain::lookup_mesh_bin(int64_t sr, Position r) const return mesh_bin; } +bool FlatSourceDomain::is_geometry_3D() +{ + // Get spatial box of ray_source_ + SpatialBox* sb = dynamic_cast( + dynamic_cast(RandomRay::ray_source_.get())->space()); + + double x_length = sb->upper_right().x - sb->lower_left().x; + double y_length = sb->upper_right().y - sb->lower_left().y; + double z_length = sb->upper_right().z - sb->lower_left().z; + + int num_xy_points = 100; + int num_z_points = 100; + uint64_t seed = openmc_get_seed(); + + for (int i = 0; i < num_xy_points; i++) { + Position sample; + sample.x = sb->lower_left().x + x_length * prn(&seed); + sample.y = sb->lower_left().y + y_length * prn(&seed); + + SourceRegionKey sr_key_prev {-1, -1}; + bool check_key = false; + + for (int j = 0; j < num_z_points; j++) { + sample.z = sb->lower_left().z + z_length * prn(&seed); + + Particle p; + p.r() = sample; + p.r_last() = sample; + p.E() = 1.0; + p.E_last() = 1.0; + p.u() = {0.0, 0.0, 1.0}; + + bool found = exhaustive_find_cell(p); + if (!found) { + continue; + } + + SourceRegionKey sr_key = lookup_source_region_key(p); + + // Check if sr_key has changed in z-direction + if (check_key && + (sr_key.base_source_region_id != sr_key_prev.base_source_region_id || + sr_key.mesh_bin != sr_key_prev.mesh_bin)) { + return true; + } + + // Set check_key to true after first sr_key has been loaded + sr_key_prev = sr_key; + check_key = true; + } + } + + return false; +} + } // namespace openmc diff --git a/src/random_ray/linear_source_domain.cpp b/src/random_ray/linear_source_domain.cpp index b4701ed1fa9..60664e12ca5 100644 --- a/src/random_ray/linear_source_domain.cpp +++ b/src/random_ray/linear_source_domain.cpp @@ -25,7 +25,6 @@ void LinearSourceDomain::batch_reset() FlatSourceDomain::batch_reset(); #pragma omp parallel for 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 diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp index dde5023e44f..951d58ffd10 100644 --- a/src/random_ray/random_ray.cpp +++ b/src/random_ray/random_ray.cpp @@ -1,5 +1,6 @@ #include "openmc/random_ray/random_ray.h" +#include "openmc/cell.h" #include "openmc/constants.h" #include "openmc/geometry.h" #include "openmc/message_passing.h" @@ -9,11 +10,11 @@ #include "openmc/search.h" #include "openmc/settings.h" #include "openmc/simulation.h" - #include #include "openmc/distribution_spatial.h" #include "openmc/random_dist.h" +#include "openmc/random_ray/decomposition_map.h" #include "openmc/source.h" namespace openmc { @@ -239,6 +240,7 @@ double RandomRay::distance_inactive_; double RandomRay::distance_active_; unique_ptr RandomRay::ray_source_; RandomRaySourceShape RandomRay::source_shape_ {RandomRaySourceShape::FLAT}; +RandomRayGeomDim RandomRay::geom_dim_ {RandomRayGeomDim::THREE_DIM}; RandomRaySampleMethod RandomRay::sample_method_ {RandomRaySampleMethod::PRNG}; RandomRay::RandomRay() @@ -261,12 +263,15 @@ RandomRay::RandomRay(uint64_t ray_id, FlatSourceDomain* domain) : RandomRay() uint64_t RandomRay::transport_history_based_single_ray() { using namespace openmc; + int n_start = n_event(); + while (alive()) { event_advance_ray(); + if (!alive()) break; event_cross_surface(); - // If ray has too many events, display warning and kill it + if (n_event() >= settings::max_particle_events) { warning("Ray " + std::to_string(id()) + " underwent maximum number of events, terminating ray."); @@ -274,7 +279,8 @@ uint64_t RandomRay::transport_history_based_single_ray() } } - return n_event(); + int delta_n = n_event() - n_start; + return delta_n; } // Transports ray across a single source region @@ -288,6 +294,18 @@ void RandomRay::event_advance_ray() boundary() = distance_to_boundary(*this); double distance = boundary().distance(); +#ifdef OPENMC_MPI + if (mpi::n_procs > 1) { + // If domain decomposition is being used, update counter for + // ray trace operations in source region for load estimation + int64_t sr = domain_->lookup_base_source_region_idx(*this); + for (int i = 0; i < n_coord(); i++) { + Cell& c {*model::cells[coord(i).cell()]}; + mpi::decomp_map.num_base_source_region_RT_[sr] += c.n_surfaces(); + } + } +#endif + if (distance < 0.0) { mark_as_lost("Negative transport distance detected for particle " + std::to_string(id())); @@ -304,8 +322,8 @@ void RandomRay::event_advance_ray() wgt() = 0.0; } - distance_travelled_ += distance; attenuate_flux(distance, true); + distance_travelled_ += distance; } else { // If the ray is still in the dead zone, need to check if it // has entered the active phase. If so, split into two segments (one @@ -313,9 +331,14 @@ void RandomRay::event_advance_ray() // first part of the active length) and attenuate each. Otherwise, if the // full length of the segment is within the dead zone, attenuate as normal. if (distance_travelled_ + distance >= distance_inactive_) { - is_active_ = true; double distance_dead = distance_inactive_ - distance_travelled_; attenuate_flux(distance_dead, false); + is_active_ = true; + distance_travelled_ = 0.0; + + if (has_left_subdomain()) { + return; + } double distance_alive = distance - distance_dead; @@ -328,8 +351,8 @@ void RandomRay::event_advance_ray() attenuate_flux(distance_alive, true, distance_dead); distance_travelled_ = distance_alive; } else { - distance_travelled_ += distance; attenuate_flux(distance, false); + distance_travelled_ += distance; } } @@ -344,6 +367,10 @@ void RandomRay::attenuate_flux(double distance, bool is_active, double offset) // Lookup base source region index int64_t sr = domain_->lookup_base_source_region_idx(*this); + // Initialize values needed to buffer ray for domain decomposition + double mesh_partial_length = 0.0; + double tiny_multiplier = 0.0; + // Perform ray tracing across mesh // Determine the mesh index for the base source region, if any int mesh_idx = domain_->lookup_mesh_idx(sr); @@ -375,17 +402,84 @@ void RandomRay::attenuate_flux(double distance, bool is_active, double offset) // Loop over all mesh bins and attenuate flux for (int b = 0; b < mesh_bins_.size(); b++) { double physical_length = reduced_distance * mesh_fractional_lengths_[b]; + +#ifdef OPENMC_MPI + if (mpi::n_procs > 1) { + mpi::decomp_map.num_mesh_bin_RT_[sr] += 1; + } +#endif + + // Very flat angles can result in very small physical lengths, + // despite the TINY_BIT adjustment for Position start. If this happens at + // an MPI boundary, this can cause rays to bounce back and forth + // indefinitely. Very small lengths are therefore skipped. + if (physical_length <= TINY_BIT) { + start += physical_length * u(); + continue; + } + attenuate_flux_inner( physical_length, is_active, sr, mesh_bins_[b], start); + start += physical_length * u(); + + // If ray has left MPI subdomain, stop transport + // and calculate position + if (has_left_subdomain()) { + for (int i = 0; i <= b - 1; i++) { + mesh_partial_length += mesh_fractional_lengths_[i]; + } + + if (b > 0) { + // If ray is stopped within mesh of base source region, + // need to add TINY_BIT to account for deleted length + tiny_multiplier = 1.0; + // Reset last surface crossed to none if ray is stopped + // within mesh of base source region + surface() = 0; + } + + mesh_partial_length = + tiny_multiplier * TINY_BIT + reduced_distance * mesh_partial_length; + break; + } } } + + // If ray has left my subdomain, buffer ray state + if (has_left_subdomain()) { + Position position_buffer = r() + (offset + mesh_partial_length) * u(); + double distance_buffer = distance_travelled_ + mesh_partial_length; + +#ifdef OPENMC_DAGMC_ENABLED + history().rollback_last_intersection(); +#endif + pack_ray_for_buffer(distance_buffer, position_buffer); + wgt() = 0.0; + } } void RandomRay::attenuate_flux_inner( double distance, bool is_active, int64_t sr, int mesh_bin, Position r) { SourceRegionKey sr_key {sr, mesh_bin}; + +#ifdef OPENMC_MPI + if (mpi::n_procs > 1) { + // Check which rank owns the source region at the current position + Position midpoint = r + u() * (distance / 2.0); + int owner = mpi::decomp_map.find_owner(SourceRegionKey(sr, mesh_bin), + midpoint, domain_->discovered_source_regions_); + + // If current rank is not the owner return and mark as not local. + if (owner != mpi::rank) { + is_local_ = false; + owner_rank_ = owner; + return; + } + } +#endif + SourceRegionHandle srh; srh = domain_->get_subdivided_source_region_handle(sr_key, r, u()); if (srh.is_numerical_fp_artifact_) { @@ -450,6 +544,7 @@ void RandomRay::attenuate_flux_flat_source( // If ray is in the active phase (not in dead zone), make contributions to // source region bookkeeping + Position midpoint = r + u() * (distance / 2.0); // Aquire lock for source region srh.lock(); @@ -464,6 +559,7 @@ void RandomRay::attenuate_flux_flat_source( // Accomulate volume (ray distance) into this iteration's estimate // of the source region's volume srh.volume() += distance; + srh.centroid_iteration() += midpoint * distance; srh.n_hits() += 1; } @@ -471,7 +567,6 @@ void RandomRay::attenuate_flux_flat_source( // Tally valid position inside the source region (e.g., midpoint of // the ray) if not done already if (!srh.position_recorded()) { - Position midpoint = r + u() * (distance / 2.0); srh.position() = midpoint; srh.position_recorded() = 1; } @@ -493,6 +588,8 @@ void RandomRay::attenuate_flux_flat_source_void( // source region bookkeeping if (is_active) { + Position midpoint = r + u() * (distance / 2.0); + // Aquire lock for source region srh.lock(); @@ -505,13 +602,13 @@ void RandomRay::attenuate_flux_flat_source_void( // Accomulate volume (ray distance) into this iteration's estimate // of the source region's volume srh.volume() += distance; + srh.centroid_iteration() += midpoint * distance; srh.volume_sq() += distance * distance; srh.n_hits() += 1; // Tally valid position inside the source region (e.g., midpoint of // the ray) if not done already if (!srh.position_recorded()) { - Position midpoint = r + u() * (distance / 2.0); srh.position() = midpoint; srh.position_recorded() = 1; } @@ -625,7 +722,7 @@ void RandomRay::attenuate_flux_linear_source( // If ray is in the active phase (not in dead zone), make contributions to // source region bookkeeping - if (is_active) { + if (is_active_) { // Accumulate deltas into the new estimate of source region flux for this // iteration for (int g = 0; g < negroups_; g++) { @@ -720,7 +817,7 @@ void RandomRay::attenuate_flux_linear_source_void( // If ray is in the active phase (not in dead zone), make contributions to // source region bookkeeping - if (is_active) { + if (is_active_) { // Compute an estimate of the spatial moments matrix for the source // region based on parameters from this ray's crossing MomentMatrix moment_matrix_estimate; @@ -768,9 +865,107 @@ void RandomRay::attenuate_flux_linear_source_void( } } +void RandomRay::restart_ray(FlatSourceDomain* domain, RayExchangeData& data, + float* angular_flux, LocalCoord* coord, int* cell_last_data) +{ + + domain_ = domain; + distance_travelled_ = data.distance_travelled; + owner_rank_ = mpi::rank; + ntemperature_ = domain->ntemperature_; + + // Restore particle event counter from the transmitted ray + // This preserves the event count across MPI rank boundaries + n_event() = data.n_event; + + is_active_ = data.is_active; + + wgt() = 1.0; + + // set identifier for particle + id() = data.ray_id; + + // Restore GeometryState scalar fields + n_coord() = data.n_coord; + cell_instance() = data.cell_instance; + n_coord_last() = data.n_coord_last; + material() = data.material; + material_last() = data.material_last; + sqrtkT() = data.sqrtkT; + sqrtkT_last() = data.sqrtkT_last; + surface() = data.surface; + + // Restore LocalCoord vector data (coord_) + // The vectors were already sized to model::n_coord_levels in the + // GeometryState constructor LocalCoord is POD so we can just copy the entire + // structure + const int n_coord_max = model::n_coord_levels; + for (int i = 0; i < n_coord_max; i++) { + this->coord(i) = coord[i]; + cell_last(i) = cell_last_data[i]; + } + + // Override the position and direction at ALL coordinate levels with the + // buffered values These may have been adjusted when the ray left the previous + // subdomain We need to update all levels to maintain consistency in the + // coordinate hierarchy + Position delta_r = data.position - r(); + for (int i = 0; i < n_coord(); i++) { + this->coord(i).r() += delta_r; + } + + u() = data.direction; + +#ifdef OPENMC_DAGMC_ENABLED + // Restore DAGMC fields + // Restore last_dir and rebuild the history by adding entities in reverse + // order + last_dir() = data.last_dir; + for (int i = data.n_handles - 1; i >= 0; i--) { + history().add_entity(data.handles[i]); + } +#endif + + // Set particle type and energy (for random ray, these are not actually used) + type() = ParticleType::neutron(); + E() = 0.0; + + // No need to call exhaustive_find_cell() since we have the full geometry + // state! Just verify we have valid cell information + if (lowest_coord().cell() == C_NONE) { + this->mark_as_lost("Received particle " + std::to_string(id()) + + " with invalid cell information"); + } + + // Set birth cell attribute if not set + if (cell_born() == C_NONE) + cell_born() = lowest_coord().cell(); + + // Set ray's angular flux to value before subdomain change + if (distance_travelled_ > 0.0 || is_active_) { + for (int g = 0; g < negroups_; g++) { + angular_flux_[g] = angular_flux[g]; + } + } + // Initialize ray's starting angular flux to starting location's isotropic + // source + else { + SourceRegionKey sr_key = domain_->lookup_source_region_key(*this); + SourceRegionHandle srh = + domain_->get_subdivided_source_region_handle(sr_key, r(), u()); + + if (!srh.is_numerical_fp_artifact_) { + for (int g = 0; g < negroups_; g++) { + angular_flux_[g] = srh.source(g); + } + } + } +} + void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain) { domain_ = domain; + owner_rank_ = mpi::rank; ntemperature_ = domain->ntemperature_; // Reset particle event counter @@ -815,11 +1010,27 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain) } SourceRegionKey sr_key = domain_->lookup_source_region_key(*this); + +#ifdef OPENMC_MPI + if (mpi::n_procs > 1) { + // Check if ray sampling site belongs to subdomain + owner_rank_ = mpi::decomp_map.find_owner( + sr_key, r(), domain_->discovered_source_regions_); + + if (owner_rank_ != mpi::rank) { + for (int g = 0; g < negroups_; g++) { + angular_flux_[g] = 0.0; + } + pack_ray_for_buffer(0.0, r()); + is_local_ = false; + return; + } + } +#endif + SourceRegionHandle srh = domain_->get_subdivided_source_region_handle(sr_key, r(), u()); - // Initialize ray's starting angular flux to starting location's isotropic - // source if (!srh.is_numerical_fp_artifact_) { for (int g = 0; g < negroups_; g++) { angular_flux_[g] = srh.source(g); @@ -877,6 +1088,65 @@ SourceSite RandomRay::sample_halton() return site; } +bool RandomRay::has_left_subdomain() +{ + return !is_local_; +} + +void RandomRay::pack_ray_for_buffer( + double distance_buffer, Position position_buffer) +{ + exchange_data_.position = position_buffer; + exchange_data_.direction = u(); + exchange_data_.angular_flux = angular_flux_; + exchange_data_.distance_travelled = distance_buffer; + exchange_data_.surface = surface(); + exchange_data_.is_active = is_active_; + exchange_data_.ray_id = id(); + exchange_data_.n_event = n_event(); + + // Pack GeometryState scalar fields + exchange_data_.n_coord = n_coord(); + exchange_data_.cell_instance = cell_instance(); + exchange_data_.n_coord_last = n_coord_last(); + exchange_data_.material = material(); + exchange_data_.material_last = material_last(); + exchange_data_.sqrtkT = sqrtkT(); + exchange_data_.sqrtkT_last = sqrtkT_last(); + + // Pack GeometryState vector fields + // LocalCoord is POD, so we can just copy the entire vector + // We always pack model::n_coord_levels elements to ensure consistent sizes + const int n_coord_max = model::n_coord_levels; + + exchange_data_.coord.resize(n_coord_max); + exchange_data_.cell_last.resize(n_coord_max); + + for (int i = 0; i < n_coord_max; i++) { + exchange_data_.coord[i] = coord(i); + exchange_data_.cell_last[i] = cell_last(i); + } + +#ifdef OPENMC_DAGMC_ENABLED + // Pack DAGMC fields + // Extract up to MAX_N_HANDLES from the ray history by rolling back + exchange_data_.last_dir = last_dir(); + exchange_data_.n_handles = 0; + + for (int i = 0; i < MAX_N_HANDLES; i++) { + moab::EntityHandle handle; + if (history().get_last_intersection(handle) == moab::MB_SUCCESS) { + exchange_data_.handles[i] = handle; + exchange_data_.n_handles++; + history().rollback_last_intersection(); + } else { + // No more handles in history + break; + } + } +#endif +} + SourceSite RandomRay::sample_s2() { // set random number seed diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index 6d6f2c743de..ed8a7fe0542 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -1,14 +1,17 @@ #include "openmc/random_ray/random_ray_simulation.h" #include "openmc/capi.h" +#include "openmc/constants.h" #include "openmc/eigenvalue.h" #include "openmc/geometry.h" #include "openmc/message_passing.h" #include "openmc/mgxs_interface.h" #include "openmc/output.h" #include "openmc/plot.h" +#include "openmc/random_ray/decomposition_map.h" #include "openmc/random_ray/flat_source_domain.h" #include "openmc/random_ray/random_ray.h" +#include "openmc/random_ray/ray_bank.h" #include "openmc/simulation.h" #include "openmc/source.h" #include "openmc/tallies/filter.h" @@ -16,6 +19,9 @@ #include "openmc/tallies/tally_scoring.h" #include "openmc/timer.h" #include "openmc/weight_windows.h" +#include +// #include +// #include namespace openmc { @@ -261,17 +267,6 @@ void validate_random_ray_inputs() } } - // Warn about slow MPI domain replication, if detected - /////////////////////////////////////////////////////////////////// -#ifdef OPENMC_MPI - if (mpi::n_procs > 1) { - warning( - "MPI parallelism is not supported by the random ray solver. All work " - "will be performed by rank 0. Domain decomposition may be implemented in " - "the future to provide efficient MPI scaling."); - } -#endif - // Warn about instability resulting from linear sources in small regions // when generating weight windows with FW-CADIS and an overlaid mesh. /////////////////////////////////////////////////////////////////// @@ -296,6 +291,7 @@ void openmc_finalize_random_ray() FlatSourceDomain::mesh_domain_map_.clear(); RandomRay::ray_source_.reset(); RandomRay::source_shape_ = RandomRaySourceShape::FLAT; + RandomRay::geom_dim_ = RandomRayGeomDim::THREE_DIM; RandomRay::sample_method_ = RandomRaySampleMethod::PRNG; } @@ -392,6 +388,13 @@ void RandomRaySimulation::simulate() // Begin main simulation timer simulation::time_total.start(); +#ifdef OPENMC_MPI + if (mpi::n_procs > 1) { + // Initialize subdomains for MPI ranks + mpi::decomp_map.initialize(); + } +#endif + // Reset per-solve accumulators, as simulate() may run more than once on the // same object (e.g. forward then adjoint when generating weight windows) avg_miss_rate_ = 0.0; @@ -403,81 +406,116 @@ void RandomRaySimulation::simulate() initialize_batch(); initialize_generation(); - // MPI not supported in random ray solver, so all work is done by rank 0 - // TODO: Implement domain decomposition for MPI parallelism - if (mpi::master) { + // Reset total starting particle weight used for normalizing tallies + simulation::total_weight = 1.0; + + // Update source term (scattering + fission) + domain_->update_all_neutron_sources(); + + // Reset scalar fluxes, iteration volume tallies, and region hit flags + // to zero + domain_->batch_reset(); + + // Check if geometry is 2D or 3D and set the appropriate flag in + // RandomRay. Generate Voronoi cells for MPI rank subdomains if domain + // decomposed. This only needs to happen once for this simulation object. + if (!geometry_setup_complete_) { - // Reset total starting particle weight used for normalizing tallies - simulation::total_weight = 1.0; + // Check if problem is 3D + if (!domain_->is_geometry_3D()) { + RandomRay::geom_dim_ = RandomRayGeomDim::TWO_DIM; + } - // Update source term (scattering + fission) - domain_->update_all_neutron_sources(); + // Generate Voronoi cells, each of which corresponds to a rank subdomain +#ifdef OPENMC_MPI + if (mpi::n_procs > 1) { + simulation::time_generate_voronoi_centers.start(); + mpi::decomp_map.generate_rank_centers(); + simulation::time_generate_voronoi_centers.stop(); + } +#endif - // Reset scalar fluxes, iteration volume tallies, and region hit flags - // to zero - domain_->batch_reset(); + geometry_setup_complete_ = true; + } - // At the beginning of the simulation, if mesh subdivision is in use, we - // need to swap the main source region container into the base container, - // as the main source region container will be used to hold the true - // subdivided source regions. The base container will therefore only - // contain the external source region information, the mesh indices, - // material properties, and initial guess values for the flux/source. + // Start timer for transport + simulation::time_transport.start(); - // Start timer for transport - simulation::time_transport.start(); + // Transport sweep over all random rays for the iteration +#ifdef OPENMC_MPI + if (mpi::n_procs > 1) { -// Transport sweep over all random rays for the iteration -#pragma omp parallel for schedule(dynamic) \ - reduction(+ : total_geometric_intersections_) - for (int i = 0; i < settings::n_particles; i++) { - RandomRay ray(i, domain_.get()); - total_geometric_intersections_ += - ray.transport_history_based_single_ray(); + // Create ray bank to store rays + RayBank RB; + + transport_sweep_decomp(RB); + + // Check if any new source regions discovered and if so, exchange + // discovered cell data between ranks + if (mpi::decomp_map.any_discovered_source_regions( + domain_->discovered_source_regions_)) { + simulation::time_source_region_exchange.start(); + mpi::decomp_map.exchange_sr_info(domain_->discovered_source_regions_); + MPI_Barrier(mpi::intracomm); + simulation::time_source_region_exchange.stop(); } - simulation::time_transport.stop(); + } else { + transport_sweep(); + } +#else + transport_sweep(); +#endif - // Add any newly discovered source regions to the main source region - // container. - domain_->finalize_discovered_source_regions(); + // Add any newly discovered source regions to the main source region + // container + domain_->finalize_discovered_source_regions(); - // Normalize scalar flux and update volumes - domain_->normalize_scalar_flux_and_volumes( - settings::n_particles * RandomRay::distance_active_); + // Normalize scalar flux and update volumes + domain_->normalize_scalar_flux_and_volumes( + settings::n_particles * RandomRay::distance_active_); - // Add source to scalar flux, compute number of FSR hits - int64_t n_hits = domain_->add_source_to_scalar_flux(); +#ifdef OPENMC_MPI + if (mpi::n_procs > 1 && simulation::current_batch <= ITER_LOAD_BALANCE) { + // Balance load between MPI ranks by exchanging source regions + simulation::time_load_balance.start(); + mpi::decomp_map.balance_load(domain_.get()); + MPI_Barrier(mpi::intracomm); + simulation::time_load_balance.stop(); + } +#endif - // Apply transport stabilization factors - domain_->apply_transport_stabilization(); + // Add source to scalar flux, compute number of FSR hits + int64_t n_hits = domain_->add_source_to_scalar_flux(); - if (settings::run_mode == RunMode::EIGENVALUE) { - // Compute random ray k-eff - domain_->compute_k_eff(); + // Apply transport stabilization factors + domain_->apply_transport_stabilization(); - // Store random ray k-eff into OpenMC's native k-eff variable - global_tally_tracklength = domain_->k_eff_; - } + if (settings::run_mode == RunMode::EIGENVALUE) { + // Compute random ray k-eff + domain_->compute_k_eff(); + + // Store random ray k-eff into OpenMC's native k-eff variable + global_tally_tracklength = domain_->k_eff_; + } - // Execute all tallying tasks, if this is an active batch - if (simulation::current_batch > settings::n_inactive) { + // Execute all tallying tasks, if this is an active batch + if (simulation::current_batch > settings::n_inactive) { - // Add this iteration's scalar flux estimate to final accumulated - // estimate - domain_->accumulate_iteration_flux(); + // Add this iteration's scalar flux estimate to final accumulated + // estimate + domain_->accumulate_iteration_flux(); - // Use above mapping to contribute FSR flux data to appropriate - // tallies - domain_->random_ray_tally(); - } + // Use above mapping to contribute FSR flux data to appropriate + // tallies + domain_->random_ray_tally(); + } - // Set phi_old = phi_new - domain_->flux_swap(); + // Set phi_old = phi_new + domain_->flux_swap(); - // Check for any obvious insabilities/nans/infs - instability_check(n_hits, domain_->k_eff_, avg_miss_rate_); - } // End MPI master work + // Check for any obvious insabilities/nans/infs + instability_check(n_hits, domain_->k_eff_, avg_miss_rate_); // Finalize the current batch finalize_generation(); @@ -507,16 +545,85 @@ void RandomRaySimulation::simulate() output_simulation_results(); } -void RandomRaySimulation::output_simulation_results() const +void RandomRaySimulation::output_simulation_results() { + + // Compute values for diagnostics + int64_t total_n_source_regions = domain_->n_source_regions(); + int64_t total_n_external_source_regions = domain_->n_external_source_regions_; + double max_load_imbalance = 0.0; + +#ifdef OPENMC_MPI + if (mpi::n_procs > 1) { + + // Average number of ray communications between ranks per batch + avg_num_communication_rounds_ = static_cast( + std::round(avg_num_communication_rounds_ / settings::n_batches)); + + // Exchange intersection data + if (mpi::master) { + MPI_Reduce(MPI_IN_PLACE, &total_geometric_intersections_, 1, MPI_UINT64_T, + MPI_SUM, 0, mpi::intracomm); + } else { + MPI_Reduce(&total_geometric_intersections_, nullptr, 1, MPI_UINT64_T, + MPI_SUM, 0, mpi::intracomm); + } + + // Exchange source region data + if (mpi::master) { + MPI_Reduce(MPI_IN_PLACE, &total_n_source_regions, 1, MPI_INT64_T, MPI_SUM, + 0, mpi::intracomm); + MPI_Reduce(MPI_IN_PLACE, &total_n_external_source_regions, 1, MPI_INT64_T, + MPI_SUM, 0, mpi::intracomm); + } else { + MPI_Reduce(&total_n_source_regions, nullptr, 1, MPI_INT64_T, MPI_SUM, 0, + mpi::intracomm); + MPI_Reduce(&total_n_external_source_regions, nullptr, 1, MPI_INT64_T, + MPI_SUM, 0, mpi::intracomm); + } + + // Determine load imbalance based on measured transport times, every rank + // needs to know about this for plotting purposes + double time_transport_total = simulation::time_transport.elapsed() - + simulation::time_ray_buffering.elapsed(); + MPI_Allgather(&time_transport_total, 1, MPI_DOUBLE, + mpi::decomp_map.measured_rank_load_fractions_.data(), 1, MPI_DOUBLE, + mpi::intracomm); + double measured_load_sum = + std::accumulate(mpi::decomp_map.measured_rank_load_fractions_.begin(), + mpi::decomp_map.measured_rank_load_fractions_.end(), 0.0); + for (int rank = 0; rank < mpi::n_procs; rank++) { + mpi::decomp_map.measured_rank_load_fractions_[rank] = + mpi::decomp_map.measured_rank_load_fractions_[rank] / measured_load_sum; + } + if (mpi::master) { + double max_load_measured = + *std::max_element(mpi::decomp_map.measured_rank_load_fractions_.begin(), + mpi::decomp_map.measured_rank_load_fractions_.end()); + max_load_imbalance = (max_load_measured - mpi::decomp_map.target_load_) / + mpi::decomp_map.target_load_; + } + } +#endif + // Print random ray results if (mpi::master) { print_results_random_ray(total_geometric_intersections_, - avg_miss_rate_ / settings::n_batches, negroups_, - domain_->n_source_regions(), domain_->n_external_source_regions_); - if (model::plots.size() > 0) { + avg_miss_rate_ / settings::n_batches, negroups_, total_n_source_regions, + total_n_external_source_regions, avg_num_communication_rounds_, + max_load_imbalance); + } + + if (model::plots.size() > 0) { +#ifdef OPENMC_MPI + if (mpi::n_procs > 1) { + domain_->output_to_vtk_decomp(); + } else { domain_->output_to_vtk(); } +#else + domain_->output_to_vtk(); +#endif } } @@ -525,12 +632,34 @@ void RandomRaySimulation::output_simulation_results() const void RandomRaySimulation::instability_check( int64_t n_hits, double k_eff, double& avg_miss_rate) const { - double percent_missed = ((domain_->n_source_regions() - n_hits) / - static_cast(domain_->n_source_regions())) * - 100.0; - avg_miss_rate += percent_missed; + + uint64_t n_source_regions = domain_->n_source_regions(); + +#ifdef OPENMC_MPI + if (mpi::n_procs > 1) { + // Reduce n_hits and n_source_regions on master rank to compute + // global miss rate + simulation::time_decomposition_handling.start(); + if (mpi::master) { + MPI_Reduce( + MPI_IN_PLACE, &n_hits, 1, MPI_INT64_T, MPI_SUM, 0, mpi::intracomm); + MPI_Reduce(MPI_IN_PLACE, &n_source_regions, 1, MPI_INT64_T, MPI_SUM, 0, + mpi::intracomm); + } else { + MPI_Reduce(&n_hits, nullptr, 1, MPI_INT64_T, MPI_SUM, 0, mpi::intracomm); + MPI_Reduce( + &n_source_regions, nullptr, 1, MPI_INT64_T, MPI_SUM, 0, mpi::intracomm); + } + simulation::time_decomposition_handling.stop(); + } +#endif if (mpi::master) { + double percent_missed = + ((n_source_regions - n_hits) / static_cast(n_source_regions)) * + 100.0; + avg_miss_rate += percent_missed; + if (percent_missed > 10.0) { warning(fmt::format( "Very high FSR miss rate detected ({:.3f}%). Instability may occur. " @@ -553,17 +682,30 @@ void RandomRaySimulation::instability_check( // Print random ray simulation results void RandomRaySimulation::print_results_random_ray( uint64_t total_geometric_intersections, double avg_miss_rate, int negroups, - int64_t n_source_regions, int64_t n_external_source_regions) const + int64_t n_source_regions, int64_t n_external_source_regions, + uint64_t avg_num_communication_rounds, double max_load_imbalance) const { using namespace simulation; if (settings::verbosity >= 6) { double total_integrations = total_geometric_intersections * negroups; - double time_per_integration = - simulation::time_transport.elapsed() / total_integrations; + + // Transport time varies between ranks, so we need to compute the total + // transport time as the sum of the max transport time across all ranks, + // i.e. the transport time of the master rank plus the time the master rank + // spends waiting for slowest other rank to finish. + double time_transport_total = time_transport.elapsed() - + time_ray_buffering.elapsed() + + time_mpi_imbalance.elapsed(); + + double time_per_integration = time_transport_total / total_integrations; + double time_domain_decomposition = + time_decomposition_handling.elapsed() + + time_generate_voronoi_centers.elapsed() + time_load_balance.elapsed() + + time_source_region_exchange.elapsed() - time_transport_total; double misc_time = time_total.elapsed() - time_update_src.elapsed() - - time_transport.elapsed() - time_tallies.elapsed() - - time_bank_sendrecv.elapsed(); + time_transport_total - time_tallies.elapsed() - + time_bank_sendrecv.elapsed() - time_domain_decomposition; header("Simulation Statistics", 4); fmt::print( @@ -591,6 +733,14 @@ void RandomRaySimulation::print_results_random_ray( fmt::print(" Avg per Iteration = {:.4e}\n", total_integrations / settings::n_batches); + if (mpi::n_procs > 1) { + fmt::print(" MPI Ranks = {}\n", mpi::n_procs); + fmt::print(" Avg Ray Subdomain Crossings = {}\n", + avg_num_communication_rounds); + fmt::print(" Maximum Load Imbalance = {:.2f}%\n", + max_load_imbalance * 100.0); + } + std::string estimator; switch (domain_->volume_estimator_) { case RandomRayVolumeEstimator::SIMULATION_AVERAGED: @@ -651,11 +801,30 @@ void RandomRaySimulation::print_results_random_ray( show_time("Total time for initialization", time_initialize.elapsed()); show_time("Reading cross sections", time_read_xs.elapsed(), 1); show_time("Total simulation time", time_total.elapsed()); - show_time("Transport sweep only", time_transport.elapsed(), 1); + if (mpi::n_procs > 1) { + show_time( + "Transport sweep only (incl. rank wait time)", time_transport_total, 1); + } else { + show_time("Transport sweep only", time_transport_total, 1); + } show_time("Source update only", time_update_src.elapsed(), 1); show_time("Tally conversion only", time_tallies.elapsed(), 1); - show_time("MPI source reductions only", time_bank_sendrecv.elapsed(), 1); + if (mpi::n_procs > 1) { + double time_decomp_misc = + time_domain_decomposition - time_source_region_exchange.elapsed() - + time_generate_voronoi_centers.elapsed() - time_ray_comms.elapsed() - + time_load_balance.elapsed(); + show_time("Decomposition handling", time_domain_decomposition, 1); + show_time("Ray communication", time_ray_comms.elapsed(), 2); + show_time( + "Exchanging contested SRs", time_source_region_exchange.elapsed(), 2); + show_time("Load balancing", time_load_balance.elapsed(), 2); + show_time("Generating Voronoi centers", + time_generate_voronoi_centers.elapsed(), 2); + show_time("Other decomposition routines", time_decomp_misc, 2); + } show_time("Other iteration routines", misc_time, 1); + if (settings::run_mode == RunMode::EIGENVALUE) { show_time("Time in inactive batches", time_inactive.elapsed()); } @@ -672,6 +841,125 @@ void RandomRaySimulation::print_results_random_ray( } } +void RandomRaySimulation::transport_sweep() +{ + + // Start timer for transport + simulation::time_transport.start(); + +// Transport sweep over all random rays for the iteration +#pragma omp parallel for schedule(dynamic) \ + reduction(+ : total_geometric_intersections_) + for (int i = 0; i < settings::n_particles; i++) { + RandomRay ray(i, domain_.get()); + total_geometric_intersections_ += ray.transport_history_based_single_ray(); + } + + simulation::time_transport.stop(); +} + +// Transport sweep for the domain decompososition case +#ifdef OPENMC_MPI +void RandomRaySimulation::transport_sweep_decomp(RayBank& RB) +{ + + double start_time_transport = simulation::time_transport.elapsed(); + double start_time_ray_buffering = simulation::time_ray_buffering.elapsed(); + + simulation::time_decomposition_handling.start(); + +// Create rays and add them to ray bank +#pragma omp parallel for schedule(static) + for (int i = 0; i < simulation::work_per_rank; i++) { + uint64_t id = simulation::work_index[mpi::rank] + i; + RandomRay ray(id, domain_.get()); + + // Add ray to ray bank if it starts in my subdomain + if (!ray.has_left_subdomain()) { +#pragma omp critical(raybank) + { + RB.my_ray_list_.push_back(ray); + } + } + // Put ray straight into buffer if it starts outside my subdomain + else { +#pragma omp critical(raybuffer) + { + RB.buffer_ray_data_to_send(ray, domain_.get()); + } + } + } + + // If no ray is alive at this stage, it means that all of them have been + // buffered because they were sampled in foregin subdomain. This requires a + // ray bank update here. + if (!RB.is_any_ray_alive()) { + RB.update(domain_.get()); + } + + int num_communication_rounds = 0; + + // Move rays across ranks until they are terminated + while (RB.is_any_ray_alive()) { + + // Start timer for transport + simulation::time_transport.start(); + +#pragma omp parallel for schedule(dynamic) \ + reduction(+ : total_geometric_intersections_) + for (int i = 0; i < RB.ray_bank_size(); i++) { + RandomRay& ray = RB.my_ray_list_[i]; + total_geometric_intersections_ += + ray.transport_history_based_single_ray(); + + // If ray has left my subdomain, buffer ray state + if (ray.has_left_subdomain()) { +#pragma omp critical(raybuffer) + { + simulation::time_ray_buffering.start(); + RB.buffer_ray_data_to_send(ray, domain_.get()); + simulation::time_ray_buffering.stop(); + } + } + } + simulation::time_transport.stop(); + + // Capture wait time resulting from other transport sweeps + simulation::time_mpi_imbalance.start(); + MPI_Barrier(mpi::intracomm); + simulation::time_mpi_imbalance.stop(); + + // Update ray bank by communicating rays in buffer to new owner ranks and + // removing terminated ranks + simulation::time_ray_comms.start(); + RB.update(domain_.get()); + MPI_Barrier(mpi::intracomm); + simulation::time_ray_comms.stop(); + + num_communication_rounds++; + } + + // Calculate load per rank based on number of hits in each source region that + // a rank owns + double batch_ray_buffering_time = + simulation::time_ray_buffering.elapsed() - start_time_ray_buffering; + double batch_transport_time = simulation::time_transport.elapsed() - + start_time_transport - batch_ray_buffering_time; + + // Calculate rank load fractions for load balancing + if (simulation::current_batch <= ITER_LOAD_BALANCE) { + mpi::decomp_map.calculate_rank_load(domain_.get(), batch_transport_time); + } + + avg_num_communication_rounds_ += num_communication_rounds; + + // Reset ray bank list for next batch + RB.my_ray_list_.resize(0); + + simulation::time_decomposition_handling.stop(); +} +#endif // OPENMC_MPI + } // namespace openmc //============================================================================== diff --git a/src/random_ray/ray_bank.cpp b/src/random_ray/ray_bank.cpp new file mode 100644 index 00000000000..07f8f4f92b5 --- /dev/null +++ b/src/random_ray/ray_bank.cpp @@ -0,0 +1,278 @@ +#include "openmc/random_ray/ray_bank.h" + +#include + +#include "openmc/geometry.h" +#include "openmc/message_passing.h" +#include "openmc/mgxs_interface.h" +#include "openmc/random_ray/decomposition_map.h" +#include "openmc/random_ray/random_ray.h" +#include "openmc/random_ray/source_region.h" +#include "openmc/timer.h" +#include + +#ifdef OPENMC_MPI + +namespace openmc { + +// Constructor +RayBank::RayBank() +{ + negroups_ = data::mg.num_energy_groups_; + num_messages_receiving_.resize(mpi::n_procs, 0); + num_messages_sending_.resize(mpi::n_procs, 0); +} + +// Buffer ray that has left my subdomain +void RayBank::buffer_ray_data_to_send(RandomRay& ray, FlatSourceDomain* domain) +{ + + // Get rank to send ray to + int rank = ray.owner_rank_; + + if (rank == mpi::rank) { + warning(fmt::format("Ray {} at position ({:.5e}, {:.5e}, {:.5e})" + "is being sent to the same rank {}. This may indicate " + "an error in the decomposition map.", + ray.exchange_data_.ray_id, ray.exchange_data_.position.x, + ray.exchange_data_.position.y, ray.exchange_data_.position.z, rank)); + } + + // Get or create the send buffer for this rank + auto& buffers = ray_send_buffer_[rank]; + + // Get the maximum coordinate levels + const int n_coord_max = model::n_coord_levels; + + // Reserve space if this is the first ray for this rank + // This is just a heuristic to reduce reallocations + if (buffers.count == 0) { + buffers.ray_data.reserve(reserved_buffer_size_); + buffers.angular_flux.reserve(reserved_buffer_size_ * negroups_); + buffers.coord.reserve(reserved_buffer_size_ * n_coord_max); + buffers.cell_last.reserve(reserved_buffer_size_ * n_coord_max); + + buffers.ray_data.reserve(32); + buffers.angular_flux.reserve(32 * negroups_); + buffers.coord.reserve(32 * n_coord_max); + buffers.cell_last.reserve(32 * n_coord_max); + } + + // Pack RayExchangeData directly into send buffer + RayBufferContainer& rbc = ray.exchange_data_; + RayExchangeData rd; + + rd.position = rbc.position; + rd.direction = rbc.direction; + rd.distance_travelled = rbc.distance_travelled; + rd.surface = rbc.surface; + rd.is_active = rbc.is_active; + rd.ray_id = rbc.ray_id; + rd.n_event = rbc.n_event; + rd.n_coord = rbc.n_coord; + rd.cell_instance = rbc.cell_instance; + rd.n_coord_last = rbc.n_coord_last; + rd.material = rbc.material; + rd.material_last = rbc.material_last; + rd.sqrtkT = rbc.sqrtkT; + rd.sqrtkT_last = rbc.sqrtkT_last; + +#ifdef OPENMC_DAGMC_ENABLED + rd.last_dir = rbc.last_dir; + rd.n_handles = rbc.n_handles; + std::memcpy( + rd.handles, rbc.handles, MAX_N_HANDLES * sizeof(moab::EntityHandle)); +#endif + + buffers.ray_data.push_back(rd); + + // Pack angular flux array directly + buffers.angular_flux.insert(buffers.angular_flux.end(), + rbc.angular_flux.begin(), rbc.angular_flux.end()); + + // Pack coord and cell_last arrays directly + buffers.coord.insert(buffers.coord.end(), rbc.coord.begin(), rbc.coord.end()); + buffers.cell_last.insert( + buffers.cell_last.end(), rbc.cell_last.begin(), rbc.cell_last.end()); + + buffers.count++; +} + +// Update ray bank +void RayBank::update(FlatSourceDomain* domain) +{ + + // Empty ray list because rays have either died or are in buffer to be sent to + // other ranks + my_ray_list_.resize(0); + + // Communicate number of rays to be sent/received between ranks + communicate_message_metadata(); + + // Send and receive ray data between MPI ranks + communicate_rays(); + + // Add received rays to ray list of that rank + update_my_ray_list(domain); +} + +int RayBank::ray_bank_size() +{ + int ray_bank_size = my_ray_list_.size(); + return ray_bank_size; +} + +// Tells each rank how many rays to receive from who and how many rays to send +// to who +void RayBank::communicate_message_metadata() +{ + + // Ensure all values are zero in vector for receiving counts + fill(num_messages_receiving_.begin(), num_messages_receiving_.end(), 0); + fill(num_messages_sending_.begin(), num_messages_sending_.end(), 0); + + // Fill the sending counts + for (auto& [rank, buffers] : ray_send_buffer_) { + num_messages_sending_[rank] = buffers.count; + } + + // Exchange message counts with all ranks + MPI_Alltoall(num_messages_sending_.data(), 1, MPI_INT, + num_messages_receiving_.data(), 1, MPI_INT, mpi::intracomm); + + total_receiving_rays_ = accumulate( + num_messages_receiving_.begin(), num_messages_receiving_.end(), 0); +} + +void RayBank::communicate_rays() +{ + + // Get the maximum coordinate levels + const int n_coord_max = model::n_coord_levels; + + // Allocate receiving buffers + received_ray_data_.resize(total_receiving_rays_); + received_angular_flux_data_.resize(total_receiving_rays_ * negroups_); + received_coord_.resize(total_receiving_rays_ * n_coord_max); + received_cell_last_.resize(total_receiving_rays_ * n_coord_max); + + // Calculate total number of MPI requests needed + // 4 messages per sending rank + 4 messages per receiving rank + int num_send_ranks = ray_send_buffer_.size(); + int num_recv_ranks = 0; + for (int i = 0; i < mpi::n_procs; i++) { + if (num_messages_receiving_[i] > 0) + num_recv_ranks++; + } + + int total_requests = num_send_ranks * 4 + num_recv_ranks * 4; + vector requests(total_requests); + int req_idx = 0; + + // Post all non-blocking receives first to allow for potential overlap + // of communication and packing of send buffers + int recv_offset = 0; + for (int sending_rank = 0; sending_rank < mpi::n_procs; sending_rank++) { + int num_rays_receiving = num_messages_receiving_[sending_rank]; + if (num_rays_receiving == 0) + continue; + + MPI_Irecv(&received_ray_data_[recv_offset], + num_rays_receiving * sizeof(RayExchangeData), MPI_BYTE, sending_rank, 1, + mpi::intracomm, &requests[req_idx++]); + + MPI_Irecv(&received_angular_flux_data_[recv_offset * negroups_], + num_rays_receiving * negroups_, MPI_FLOAT, sending_rank, 2, + mpi::intracomm, &requests[req_idx++]); + + MPI_Irecv(&received_coord_[recv_offset * n_coord_max], + num_rays_receiving * n_coord_max * sizeof(LocalCoord), MPI_BYTE, + sending_rank, 3, mpi::intracomm, &requests[req_idx++]); + + MPI_Irecv(&received_cell_last_[recv_offset * n_coord_max], + num_rays_receiving * n_coord_max, MPI_INT, sending_rank, 4, + mpi::intracomm, &requests[req_idx++]); + + recv_offset += num_rays_receiving; + } + + // Post all non-blocking sends! + for (auto& [receiving_rank, buffers] : ray_send_buffer_) { + int num_rays = buffers.count; + + // Check neighbor list and add if not already known + // Filter out rays that are sampled elsewhere + bool has_transported_rays = false; + for (int i = 0; i < num_rays; i++) { + if (buffers.ray_data[i].distance_travelled > 0.0 || + buffers.ray_data[i].is_active) { + has_transported_rays = true; + break; + } + } + if (has_transported_rays) { + mpi::decomp_map.my_neighbors_.insert(receiving_rank); + } + + // Send all 4 data arrays - already packed, no intermediate copying! + MPI_Isend(buffers.ray_data.data(), num_rays * sizeof(RayExchangeData), + MPI_BYTE, receiving_rank, 1, mpi::intracomm, &requests[req_idx++]); + + MPI_Isend(buffers.angular_flux.data(), num_rays * negroups_, MPI_FLOAT, + receiving_rank, 2, mpi::intracomm, &requests[req_idx++]); + + MPI_Isend(buffers.coord.data(), num_rays * n_coord_max * sizeof(LocalCoord), + MPI_BYTE, receiving_rank, 3, mpi::intracomm, &requests[req_idx++]); + + MPI_Isend(buffers.cell_last.data(), num_rays * n_coord_max, MPI_INT, + receiving_rank, 4, mpi::intracomm, &requests[req_idx++]); + } + + // Wait for all communication to complete + MPI_Waitall(req_idx, requests.data(), MPI_STATUSES_IGNORE); + + // Clear send buffers + ray_send_buffer_.clear(); +} + +void RayBank::update_my_ray_list(FlatSourceDomain* domain) +{ + + my_ray_list_.resize(received_ray_data_.size()); + + const int n_coord_max = model::n_coord_levels; + +// Add re-initialized random ray objects to my_ray_list +#pragma omp parallel for + for (int i = 0; i < received_ray_data_.size(); i++) { + + // Re-initialize rays with received data, including full geometry state + my_ray_list_[i].restart_ray(domain, received_ray_data_[i], + &received_angular_flux_data_[i * negroups_], + &received_coord_[i * n_coord_max], &received_cell_last_[i * n_coord_max]); + } + + // Clear received data vectors + received_ray_data_.resize(0); + received_angular_flux_data_.resize(0); + received_coord_.resize(0); + received_cell_last_.resize(0); +} + +bool RayBank::is_any_ray_alive() +{ + + int local_rays_alive = ray_bank_size(); + int flag = 0; + if (local_rays_alive > 0) { + flag = 1; + } + + MPI_Allreduce(MPI_IN_PLACE, &flag, 1, MPI_INT, MPI_MAX, mpi::intracomm); + + return flag > 0; +} + +} // namespace openmc + +#endif // OPENMC_MPI diff --git a/src/random_ray/source_region.cpp b/src/random_ray/source_region.cpp index 78543c5ab53..a893ce401b1 100644 --- a/src/random_ray/source_region.cpp +++ b/src/random_ray/source_region.cpp @@ -10,19 +10,24 @@ namespace openmc { // SourceRegionHandle implementation //============================================================================== SourceRegionHandle::SourceRegionHandle(SourceRegion& sr) - : negroups_(sr.scalar_flux_old_.size()), material_(&sr.material_), - temperature_idx_(&sr.temperature_idx_), density_mult_(&sr.density_mult_), - is_small_(&sr.is_small_), n_hits_(&sr.n_hits_), + : negroups_(sr.scalar_flux_old_.size()), material_(&sr.scalars_.material_), + temperature_idx_(&sr.scalars_.temperature_idx_), + density_mult_(&sr.scalars_.density_mult_), + is_small_(&sr.scalars_.is_small_), n_hits_(&sr.scalars_.n_hits_), is_linear_(sr.source_gradients_.size() > 0), lock_(&sr.lock_), - volume_(&sr.volume_), volume_t_(&sr.volume_t_), volume_sq_(&sr.volume_sq_), - volume_sq_t_(&sr.volume_sq_t_), volume_naive_(&sr.volume_naive_), - position_recorded_(&sr.position_recorded_), - external_source_present_(&sr.external_source_present_), - position_(&sr.position_), centroid_(&sr.centroid_), - centroid_iteration_(&sr.centroid_iteration_), centroid_t_(&sr.centroid_t_), - mom_matrix_(&sr.mom_matrix_), mom_matrix_t_(&sr.mom_matrix_t_), - volume_task_(&sr.volume_task_), mesh_(&sr.mesh_), - parent_sr_(&sr.parent_sr_), scalar_flux_old_(sr.scalar_flux_old_.data()), + volume_(&sr.scalars_.volume_), volume_t_(&sr.scalars_.volume_t_), + volume_sq_(&sr.scalars_.volume_sq_), + volume_sq_t_(&sr.scalars_.volume_sq_t_), + volume_naive_(&sr.scalars_.volume_naive_), + position_recorded_(&sr.scalars_.position_recorded_), + external_source_present_(&sr.scalars_.external_source_present_), + position_(&sr.scalars_.position_), centroid_(&sr.scalars_.centroid_), + centroid_iteration_(&sr.scalars_.centroid_iteration_), + centroid_t_(&sr.scalars_.centroid_t_), + mom_matrix_(&sr.scalars_.mom_matrix_), + mom_matrix_t_(&sr.scalars_.mom_matrix_t_), volume_task_(&sr.volume_task_), + mesh_(&sr.scalars_.mesh_), parent_sr_(&sr.scalars_.parent_sr_), + scalar_flux_old_(sr.scalar_flux_old_.data()), scalar_flux_new_(sr.scalar_flux_new_.data()), source_(sr.source_.data()), external_source_(sr.external_source_.data()), scalar_flux_final_(sr.scalar_flux_final_.data()), @@ -30,7 +35,7 @@ SourceRegionHandle::SourceRegionHandle(SourceRegion& sr) flux_moments_old_(sr.flux_moments_old_.data()), flux_moments_new_(sr.flux_moments_new_.data()), flux_moments_t_(sr.flux_moments_t_.data()), - tally_task_(sr.tally_task_.data()) + tally_task_(sr.tally_task_.data()), key_(&sr.scalars_.key_) {} //============================================================================== @@ -61,6 +66,88 @@ SourceRegion::SourceRegion(int negroups, bool is_linear) } } +SourceRegion::SourceRegion(const SourceRegionHandle& handle) + : SourceRegion(handle.negroups_, handle.is_linear_) +{ + scalars_.material_ = handle.material(); + scalars_.is_small_ = handle.is_small(); + scalars_.temperature_idx_ = handle.temperature_idx(); + scalars_.density_mult_ = handle.density_mult(); + scalars_.n_hits_ = handle.n_hits(); + scalars_.volume_ = handle.volume(); + scalars_.volume_t_ = handle.volume_t(); + scalars_.volume_sq_ = handle.volume_sq(); + scalars_.volume_sq_t_ = handle.volume_sq_t(); + scalars_.volume_naive_ = handle.volume_naive(); + scalars_.position_recorded_ = handle.position_recorded(); + scalars_.position_ = handle.position(); + scalars_.centroid_ = handle.centroid(); + scalars_.centroid_iteration_ = handle.centroid_iteration(); + scalars_.centroid_t_ = handle.centroid_t(); + scalars_.key_ = handle.key(); + scalars_.mesh_ = handle.mesh(); + scalars_.parent_sr_ = handle.parent_sr(); + + if (handle.is_linear_) { + scalars_.mom_matrix_ = handle.mom_matrix(); + scalars_.mom_matrix_t_ = handle.mom_matrix_t(); + } + + for (int g = 0; g < scalar_flux_new_.size(); g++) { + scalar_flux_old_[g] = handle.scalar_flux_old(g); + scalar_flux_new_[g] = handle.scalar_flux_new(g); + source_[g] = handle.source(g); + scalar_flux_final_[g] = handle.scalar_flux_final(g); + if (handle.is_linear_) { + source_gradients_[g] = handle.source_gradients(g); + flux_moments_old_[g] = handle.flux_moments_old(g); + flux_moments_new_[g] = handle.flux_moments_new(g); + flux_moments_t_[g] = handle.flux_moments_t(g); + } + tally_task_[g] = handle.tally_task(g); + } + if (settings::run_mode == RunMode::FIXED_SOURCE) { + scalars_.external_source_present_ = handle.external_source_present(); + for (int g = 0; g < scalar_flux_new_.size(); g++) { + external_source_[g] = handle.external_source(g); + } + } + + volume_task_ = handle.volume_task(); +} + +// combine two source regions from different ranks together +void SourceRegion::merge(SourceRegion& sr_add, bool is_linear) +{ + + // scalar fields + scalars_.volume_ += sr_add.scalars_.volume_; + scalars_.volume_sq_ += sr_add.scalars_.volume_sq_; + scalars_.volume_naive_ += sr_add.scalars_.volume_naive_; + scalars_.n_hits_ += sr_add.scalars_.n_hits_; + scalars_.external_source_present_ = + std::max(scalars_.external_source_present_, + sr_add.scalars_.external_source_present_); + scalars_.centroid_iteration_ += sr_add.scalars_.centroid_iteration_; + if (is_linear) { + scalars_.mom_matrix_ += sr_add.scalars_.mom_matrix_; + } + +// vector fields +#pragma omp simd + for (int g = 0; g < scalar_flux_new_.size(); g++) { + scalar_flux_new_[g] += sr_add.scalar_flux_new_[g]; + scalar_flux_final_[g] += sr_add.scalar_flux_final_[g]; + + if (settings::run_mode == RunMode::FIXED_SOURCE) { + external_source_[g] += sr_add.external_source_[g]; + } + if (is_linear) { + flux_moments_new_[g] += sr_add.flux_moments_new_[g]; + } + } +} + //============================================================================== // SourceRegionContainer implementation //============================================================================== @@ -70,31 +157,32 @@ void SourceRegionContainer::push_back(const SourceRegion& sr) n_source_regions_++; // Scalar fields - material_.push_back(sr.material_); - temperature_idx_.push_back(sr.temperature_idx_); - density_mult_.push_back(sr.density_mult_); - is_small_.push_back(sr.is_small_); - n_hits_.push_back(sr.n_hits_); + material_.push_back(sr.scalars_.material_); + temperature_idx_.push_back(sr.scalars_.temperature_idx_); + density_mult_.push_back(sr.scalars_.density_mult_); + is_small_.push_back(sr.scalars_.is_small_); + n_hits_.push_back(sr.scalars_.n_hits_); lock_.push_back(sr.lock_); - volume_.push_back(sr.volume_); - volume_t_.push_back(sr.volume_t_); - volume_sq_.push_back(sr.volume_sq_); - volume_sq_t_.push_back(sr.volume_sq_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_.push_back(sr.scalars_.volume_); + volume_t_.push_back(sr.scalars_.volume_t_); + volume_sq_.push_back(sr.scalars_.volume_sq_); + volume_sq_t_.push_back(sr.scalars_.volume_sq_t_); + volume_naive_.push_back(sr.scalars_.volume_naive_); + position_recorded_.push_back(sr.scalars_.position_recorded_); + external_source_present_.push_back(sr.scalars_.external_source_present_); + position_.push_back(sr.scalars_.position_); volume_task_.push_back(sr.volume_task_); - mesh_.push_back(sr.mesh_); - parent_sr_.push_back(sr.parent_sr_); + mesh_.push_back(sr.scalars_.mesh_); + parent_sr_.push_back(sr.scalars_.parent_sr_); + key_.push_back(sr.scalars_.key_); + centroid_.push_back(sr.scalars_.centroid_); + centroid_iteration_.push_back(sr.scalars_.centroid_iteration_); + centroid_t_.push_back(sr.scalars_.centroid_t_); // 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_); + mom_matrix_.push_back(sr.scalars_.mom_matrix_); + mom_matrix_t_.push_back(sr.scalars_.mom_matrix_t_); } // Energy-dependent fields @@ -141,11 +229,12 @@ void SourceRegionContainer::assign( position_.clear(); mesh_.clear(); parent_sr_.clear(); + centroid_.clear(); + centroid_iteration_.clear(); + centroid_t_.clear(); + key_.clear(); if (is_linear_) { - centroid_.clear(); - centroid_iteration_.clear(); - centroid_t_.clear(); mom_matrix_.clear(); mom_matrix_t_.clear(); } @@ -212,11 +301,12 @@ SourceRegionHandle SourceRegionContainer::get_source_region_handle(int64_t sr) } handle.scalar_flux_final_ = &scalar_flux_final(sr, 0); handle.tally_task_ = &tally_task(sr, 0); + handle.centroid_ = ¢roid(sr); + handle.centroid_iteration_ = ¢roid_iteration(sr); + handle.centroid_t_ = ¢roid_t(sr); + handle.key_ = &key(sr); if (handle.is_linear_) { - handle.centroid_ = ¢roid(sr); - handle.centroid_iteration_ = ¢roid_iteration(sr); - handle.centroid_t_ = ¢roid_t(sr); handle.mom_matrix_ = &mom_matrix(sr); handle.mom_matrix_t_ = &mom_matrix_t(sr); handle.source_gradients_ = &source_gradients(sr, 0); @@ -230,6 +320,10 @@ SourceRegionHandle SourceRegionContainer::get_source_region_handle(int64_t sr) void SourceRegionContainer::adjoint_reset() { + // Note: key_ must NOT be reset here. Source region keys are permanent + // identifiers and required under domain decomposition to rebuild + // source_region_map_ during load balancing and to exchange source region data + // between ranks. std::fill(n_hits_.begin(), n_hits_.end(), 0); std::fill(volume_.begin(), volume_.end(), 0.0); std::fill(volume_t_.begin(), volume_t_.end(), 0.0); diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 08563324761..04ded9f5ace 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -1098,7 +1098,7 @@ void accumulate_tallies() { #ifdef OPENMC_MPI // Combine tally results onto master process - if (mpi::n_procs > 1 && settings::solver_type == SolverType::MONTE_CARLO) { + if (mpi::n_procs > 1) { reduce_tally_results(); } #endif diff --git a/src/timer.cpp b/src/timer.cpp index 6d692d4fbf6..c959da29af2 100644 --- a/src/timer.cpp +++ b/src/timer.cpp @@ -27,6 +27,13 @@ Timer time_event_surface_crossing; Timer time_event_collision; Timer time_event_death; Timer time_update_src; +Timer time_ray_comms; +Timer time_ray_buffering; +Timer time_decomposition_handling; +Timer time_load_balance; +Timer time_generate_voronoi_centers; +Timer time_source_region_exchange; +Timer time_mpi_imbalance; } // namespace simulation