From 1bbe253678641fc92a7a8d17c70a8dc6840b5547 Mon Sep 17 00:00:00 2001 From: "W. Daryl Hawkins" Date: Tue, 2 Jun 2026 13:03:08 -0500 Subject: [PATCH] Moving sweep runtime construction out of DiscreteOrdinatesProblem --- .../discrete_ordinates_problem.cc | 441 +-------------- .../discrete_ordinates_problem.cu | 28 - .../discrete_ordinates_problem.h | 16 - .../sweep/spds/aah.h | 7 +- .../sweep/sweep_runtime_builder.cc | 511 ++++++++++++++++++ .../sweep/sweep_runtime_builder.cu | 34 ++ .../sweep/sweep_runtime_builder.h | 55 ++ 7 files changed, 611 insertions(+), 481 deletions(-) create mode 100644 modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/sweep_runtime_builder.cc create mode 100644 modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/sweep_runtime_builder.cu create mode 100644 modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/sweep_runtime_builder.h diff --git a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.cc b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.cc index 6ad697b0bd..b5a0cfb273 100644 --- a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.cc +++ b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.cc @@ -13,6 +13,7 @@ #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/spds/aah.h" #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/fluds/aah_fluds.h" #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/angle_set/aah_angle_set.h" +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/sweep_runtime_builder.h" #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep_chunks/aah_sweep_chunk.h" #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep_chunks/aah_sweep_chunk_td.h" #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep_chunks/cbc_sweep_chunk.h" @@ -129,11 +130,6 @@ DiscreteOrdinatesProblem::GetInputParameters() "boundary_conditions", {}, "An array containing tables for each boundary specification."); params.LinkParameterToBlock("boundary_conditions", "BoundaryOptionsBlock"); - params.AddOptionalParameterArray( - "directions_sweep_order_to_print", - std::vector(), - "List of direction id's for which sweep ordering info is to be printed."); - params.AddOptionalParameter( "sweep_type", "AAH", "The sweep type to use for sweep operatorations."); params.ConstrainParameterRange("sweep_type", AllowableRangeList::New({"AAH", "CBC"})); @@ -208,9 +204,7 @@ DiscreteOrdinatesProblem::Create(const ParameterBlock& params) } DiscreteOrdinatesProblem::DiscreteOrdinatesProblem(const InputParameters& params) - : LBSProblem(params), - verbose_sweep_angles_(params.GetParamVectorValue("directions_sweep_order_to_print")), - sweep_type_(params.GetParamValue("sweep_type")) + : LBSProblem(params), sweep_type_(params.GetParamValue("sweep_type")) { if (params.GetParamValue("time_dependent")) { @@ -557,9 +551,6 @@ DiscreteOrdinatesProblem::BuildRuntime() LBSProblem::BuildRuntime(); InitializeFCS(); - // Make face histogram - grid_face_histogram_ = grid_->MakeGridFaceHistogram(); - UpdateAngularFluxStorage(); const auto grid_dim = grid_->GetDimension(); @@ -1328,235 +1319,12 @@ DiscreteOrdinatesProblem::InitializeSweepDataStructures() log.Log() << program_timer.GetTimeString() << " Initializing sweep datastructures.\n"; - // Define sweep ordering groups - quadrature_unq_so_grouping_map_.clear(); - std::map, bool> quadrature_allow_cycles_map_; - for (auto& groupset : groupsets_) - { - if (quadrature_unq_so_grouping_map_.count(groupset.quadrature) == 0) - { - quadrature_unq_so_grouping_map_[groupset.quadrature] = AssociateSOsAndDirections( - grid_, *groupset.quadrature, groupset.angleagg_method, geometry_type_); - } - - if (quadrature_allow_cycles_map_.count(groupset.quadrature) == 0) - quadrature_allow_cycles_map_[groupset.quadrature] = groupset.allow_cycles; - } - - // Build sweep orderings - quadrature_spds_map_.clear(); - if (sweep_type_ == "AAH") - { - // Creating an AAH SPDS can be an expensive operation. We break it up into multiple phases so - // so that we can distribute the work across MPI ranks: - // 1) Initialize the SPDS for each angleset. This is done by all ranks. - // 2) For each SPDS, generate the feedback arc set (FAS) for the global sweep graph. Angelsets - // are distributed as evenly as possible across MPI ranks. - // 3) Gather the FAS for each SPDS on all ranks. - // 4) Build the global sweep task dependency graph (TDG) for each SPDS. - - // Initalize SPDS. All ranks initialize a SPDS for each angleset. - log.Log0Verbose1() << program_timer.GetTimeString() << " Initializing AAH SPDS."; - for (const auto& [quadrature, info] : quadrature_unq_so_grouping_map_) - { - int id = 0; - const auto& unique_so_groupings = info.first; - for (const auto& so_grouping : unique_so_groupings) - { - if (so_grouping.empty()) - continue; - - const size_t master_dir_id = so_grouping.front(); - const auto& omega = quadrature->GetOmega(master_dir_id); - const auto new_swp_order = std::make_shared( - id, omega, this->grid_, quadrature_allow_cycles_map_[quadrature], use_gpus_); - quadrature_spds_map_[quadrature].push_back(new_swp_order); - ++id; - } - } - - // Accumulate global edge weights for each SPDS on the owning rank only. - const int comm_size = opensn::mpi_comm.size(); - const int matrix_size = comm_size * comm_size; - std::vector recv_counts(opensn::mpi_comm.size(), comm_size); - std::vector recv_displacements(opensn::mpi_comm.size(), 0); - for (int loc = 0; loc < opensn::mpi_comm.size(); ++loc) - recv_displacements[loc] = loc * comm_size; - for (const auto& [quadrature, spds_list] : quadrature_spds_map_) - { - for (const auto& spds : spds_list) - { - auto aah_spds = std::static_pointer_cast(spds); - const int owner = aah_spds->GetId() % opensn::mpi_comm.size(); - - // Local contributions - weights from this rank to all others for this SPDS - const auto local_row = aah_spds->ComputeLocalLocationEdgeWeights(); - std::vector recv; - if (opensn::mpi_comm.rank() == owner) - recv.assign(matrix_size, 0.0); - opensn::mpi_comm.gather(local_row, recv, recv_counts, recv_displacements, owner); - - if (opensn::mpi_comm.rank() == owner) - aah_spds->SetGlobalEdgeWeights(recv); - } - } - - // Generate the global sweep FAS for each SPDS. This is an expensive operation. It is - // distributed via MPI so that multiple MPI ranks can compute the FAS for one or more SPDS - // independently. - log.Log0Verbose1() << program_timer.GetTimeString() << " Build global sweep FAS for each SPDS."; - for (const auto& quadrature : quadrature_spds_map_) - { - for (const auto& spds : quadrature.second) - { - auto aah_spds = std::static_pointer_cast(spds); - auto id = aah_spds->GetId(); - if (opensn::mpi_comm.rank() == (id % opensn::mpi_comm.size())) - aah_spds->BuildGlobalSweepFAS(); - } - } - - // Communicate the FAS for each SPDS to all ranks. - log.Log0Verbose1() << program_timer.GetTimeString() << " Gather FAS for each SPDS."; - std::vector local_edges_to_remove; - for (const auto& quadrature : quadrature_spds_map_) - { - for (const auto& spds : quadrature.second) - { - auto aah_spds = std::static_pointer_cast(spds); - auto id = aah_spds->GetId(); - if ((id % opensn::mpi_comm.size()) == opensn::mpi_comm.rank()) - { - auto edges_to_remove = aah_spds->GetGlobalSweepFAS(); - local_edges_to_remove.push_back(id); - local_edges_to_remove.push_back(static_cast(edges_to_remove.size())); - local_edges_to_remove.insert( - local_edges_to_remove.end(), edges_to_remove.begin(), edges_to_remove.end()); - } - } - } - - int local_size = static_cast(local_edges_to_remove.size()); - std::vector receive_counts(opensn::mpi_comm.size(), 0); - std::vector displacements(opensn::mpi_comm.size(), 0); - mpi_comm.all_gather(local_size, receive_counts); - - int total_size = 0; - for (auto i = 0; i < receive_counts.size(); ++i) - { - displacements[i] = total_size; - total_size += receive_counts[i]; - } - - std::vector global_edges_to_remove(total_size, 0); - mpi_comm.all_gather( - local_edges_to_remove, global_edges_to_remove, receive_counts, displacements); - - // Unpack the gathered data and update SPDS on all ranks. - int offset = 0; - while (offset < global_edges_to_remove.size()) - { - auto spds_id = global_edges_to_remove[offset++]; - auto num_edges = global_edges_to_remove[offset++]; - std::vector edges; - edges.reserve(num_edges); - for (auto i = 0; i < num_edges; ++i) - edges.emplace_back(global_edges_to_remove[offset++]); - - for (const auto& quadrature : quadrature_spds_map_) - { - for (const auto& spds : quadrature.second) - { - auto aah_spds = std::static_pointer_cast(spds); - if (aah_spds->GetId() == spds_id) - { - aah_spds->SetGlobalSweepFAS(edges); - break; - } - } - } - } + auto sweep_runtime = BuildSweepRuntime( + GetName(), groupsets_, grid_, sweep_type_, use_gpus_, *discretization_, grid_nodal_mappings_); - // Build TDG for each SPDS on all ranks. - log.Log0Verbose1() << program_timer.GetTimeString() << " Build global sweep TDGs."; - for (const auto& quadrature : quadrature_spds_map_) - for (const auto& spds : quadrature.second) - std::static_pointer_cast(spds)->BuildGlobalSweepTDG(); - - // Print ghosted sweep graph if requested - if (not verbose_sweep_angles_.empty()) - { - for (const auto& quadrature : quadrature_spds_map_) - { - for (const auto& spds : quadrature.second) - { - for (const int dir_id : verbose_sweep_angles_) - { - auto aah_spds = std::static_pointer_cast(spds); - if (aah_spds->GetId() == dir_id) - aah_spds->PrintGhostedGraph(); - } - } - } - } - } - else if (sweep_type_ == "CBC") - { - // Build SPDS - for (const auto& [quadrature, info] : quadrature_unq_so_grouping_map_) - { - const auto& unique_so_groupings = info.first; - for (const auto& so_grouping : unique_so_groupings) - { - if (so_grouping.empty()) - continue; - - const size_t master_dir_id = so_grouping.front(); - const auto& omega = quadrature->GetOmega(master_dir_id); - const auto new_swp_order = - std::make_shared(omega, this->grid_, quadrature_allow_cycles_map_[quadrature]); - quadrature_spds_map_[quadrature].push_back(new_swp_order); - } - } - } - else - OpenSnInvalidArgument("Unsupported sweep type \"" + sweep_type_ + "\""); - - opensn::mpi_comm.barrier(); - - // Build FLUDS templates - quadrature_fluds_commondata_map_.clear(); - if (sweep_type_ == "AAH" && use_gpus_) - { - CreateAAHD_FLUDSCommonData(); - } - else if (sweep_type_ == "AAH") - { - for (const auto& [quadrature, spds_list] : quadrature_spds_map_) - { - for (const auto& spds : spds_list) - { - quadrature_fluds_commondata_map_[quadrature].push_back( - std::make_unique( - grid_nodal_mappings_, *spds, *grid_face_histogram_)); - } - } - } - else if (sweep_type_ == "CBC" and use_gpus_) - { - CreateCBCD_FLUDSCommonData(); - } - else if (sweep_type_ == "CBC") - { - for (const auto& [quadrature, spds_list] : quadrature_spds_map_) - { - for (const auto& spds : spds_list) - { - quadrature_fluds_commondata_map_[quadrature].push_back( - std::make_unique(*spds, grid_nodal_mappings_)); - } - } - } + quadrature_unq_so_grouping_map_ = std::move(sweep_runtime.quadrature_unq_so_grouping_map); + quadrature_spds_map_ = std::move(sweep_runtime.quadrature_spds_map); + quadrature_fluds_commondata_map_ = std::move(sweep_runtime.quadrature_fluds_commondata_map); log.Log() << program_timer.GetTimeString() << " Done initializing sweep datastructures.\n"; } @@ -1579,13 +1347,6 @@ DiscreteOrdinatesProblem::ResetBoundaryCarrier() { } -void -DiscreteOrdinatesProblem::CreateAAHD_FLUDSCommonData() -{ - throw std::runtime_error( - "DiscreteOrdinatesProblem::CreateAAHD_FLUDSCommonData : OPENSN_WITH_CUDA not enabled."); -} - void DiscreteOrdinatesProblem::UpdateAAHD_FLUDSCommonDataWithBoundary() { @@ -1637,13 +1398,6 @@ DiscreteOrdinatesProblem::CopyPhiAndOutflowBackToHost() { } -void -DiscreteOrdinatesProblem::CreateCBCD_FLUDSCommonData() -{ - throw std::runtime_error( - "DiscreteOrdinatesProblem::CreateCBCD_FLUDSCommonData : OPENSN_WITH_CUDA not enabled."); -} - std::shared_ptr DiscreteOrdinatesProblem::CreateCBCD_FLUDS(std::size_t num_groups, std::size_t num_angles, @@ -1682,187 +1436,6 @@ DiscreteOrdinatesProblem::CreateCBCDSweepChunk(LBSGroupset& groupset) } #endif -std::pair -DiscreteOrdinatesProblem::AssociateSOsAndDirections(const std::shared_ptr grid, - const AngularQuadrature& quadrature, - const AngleAggregationType agg_type, - const GeometryType lbs_geo_type) -{ - - // Checks - if (quadrature.GetOmegas().empty()) - throw std::logic_error(GetName() + ": Quadrature with no omegas cannot be used"); - if (quadrature.GetWeights().empty()) - throw std::logic_error(GetName() + ": Quadrature with no weights cannot be used"); - - // Build groupings - UniqueSOGroupings unq_so_grps; - switch (agg_type) - { - // SINGLE AGGREGATION - // The simplest aggregation type. Every direction is assumed to have a unique sweep ordering. - // There are as many direction sets as there are directions. - case AngleAggregationType::SINGLE: - { - if (lbs_geo_type == GeometryType::TWOD_CYLINDRICAL) - { - // Preserve azimuthal ordering per polar level - const auto* product_quad = dynamic_cast(&quadrature); - if (product_quad) - { - for (const auto& dir_set : product_quad->GetDirectionMap()) - for (const auto dir_id : dir_set.second) - unq_so_grps.push_back({dir_id}); - } - else - { - const size_t num_dirs = quadrature.GetNumAngles(); - for (size_t n = 0; n < num_dirs; ++n) - unq_so_grps.push_back({n}); - } - } - else - { - const size_t num_dirs = quadrature.GetNumAngles(); - for (size_t n = 0; n < num_dirs; ++n) - unq_so_grps.push_back({n}); - } - break; - } // case agg_type SINGLE - - // POLAR AGGREGATION - // Aggregate all polar directions for a given azimuthal direction into a direction set. - case AngleAggregationType::POLAR: - { - // Check geometry types - if (grid->GetType() != ORTHOGONAL and grid->GetDimension() != 2 and not grid->Extruded()) - throw std::logic_error( - GetName() + ": The simulation is using polar angle aggregation for which only certain " - "geometry types are supported, i.e., ORTHOGONAL, 2D or 3D EXTRUDED"); - - // Check quadrature type - const auto quad_type = quadrature.GetType(); - if (quad_type != AngularQuadratureType::PRODUCT_QUADRATURE) - throw std::logic_error(GetName() + ": The simulation is using polar angle aggregation for " - "which only Product-type quadratures are supported"); - - // Process Product Quadrature - try - { - const auto& product_quad = dynamic_cast(quadrature); - - const auto& azimuthal_angles = product_quad.GetAzimuthalAngles(); - const auto& polar_angles = product_quad.GetPolarAngles(); - const auto num_azi = azimuthal_angles.size(); - const auto num_pol = polar_angles.size(); - - // Make two separate list of polar angles - // One upward-pointing and one downward - std::vector upward_polar_ids; - std::vector dnward_polar_ids; - for (size_t p = 0; p < num_pol; ++p) - if (polar_angles[p] > M_PI_2) - upward_polar_ids.push_back(p); - else - dnward_polar_ids.push_back(p); - - // Define lambda working for both upward and dnward polar-ids - /**Lambda to convert indices and push it onto unq_so_grps.*/ - auto MapPolarAndAzimuthalIDs = - [&product_quad, &unq_so_grps](const DirIDs& polar_ids, const size_t azimuthal_id) - { - DirIDs dir_ids; - dir_ids.reserve(polar_ids.size()); - for (const size_t p : polar_ids) - dir_ids.push_back(product_quad.GetAngleNum(p, azimuthal_id)); - unq_so_grps.push_back(std::move(dir_ids)); - }; - - // Stack id's for all azimuthal angles - for (size_t a = 0; a < num_azi; ++a) - { - if (not upward_polar_ids.empty()) - MapPolarAndAzimuthalIDs(upward_polar_ids, a); - if (not dnward_polar_ids.empty()) - MapPolarAndAzimuthalIDs(dnward_polar_ids, a); - } // for azi-id a - - } // try product quadrature - catch (const std::bad_cast& bc) - { - throw std::runtime_error( - GetName() + ": Casting the angular quadrature to the product quadrature base failed"); - } - - break; - } // case agg_type POLAR - - // AZIMUTHAL AGGREGATION - // All azimuthal direction in a quadrant/octant are assigned to a direction set - case AngleAggregationType::AZIMUTHAL: - { - // Check geometry types - if (lbs_geo_type != GeometryType::ONED_SPHERICAL and - lbs_geo_type != GeometryType::TWOD_CYLINDRICAL) - throw std::logic_error( - GetName() + ": AZIMUTHAL aggregation is only valid for TWOD_CYLINDRICAL geometry"); - - // Check quadrature type - const auto quad_type = quadrature.GetType(); - if (quad_type != AngularQuadratureType::PRODUCT_QUADRATURE) - throw std::logic_error( - GetName() + ": AZIMUTHAL aggregation is only valid for TWOD_CYLINDRICAL geometry."); - - // Process Product Quadrature - try - { - const auto& product_quad = dynamic_cast(quadrature); - - for (const auto& dir_set : product_quad.GetDirectionMap()) - { - std::vector group1; - std::vector group2; - for (const auto& dir_id : dir_set.second) - if (quadrature.GetAbscissa(dir_id).phi > M_PI_2) - group1.push_back(dir_id); - else - group2.push_back(dir_id); - - DirIDs group1_ids(group1.begin(), group1.end()); - DirIDs group2_ids(group2.begin(), group2.end()); - - unq_so_grps.push_back(std::move(group1_ids)); - unq_so_grps.push_back(std::move(group2_ids)); - } - } // try product quadrature - catch (const std::bad_cast& bc) - { - throw std::runtime_error( - GetName() + ": Casting the angular quadrature to the product quadrature base failed"); - } - - break; - } - default: - throw std::invalid_argument(GetName() + ": Called with UNDEFINED angle aggregation type"); - } // switch angle aggregation type - - // Map directions to sweep orderings - DirIDToSOMap dir_id_to_so_map; - { - size_t so_grouping_id = 0; - for (const auto& so_grouping : unq_so_grps) - { - for (const size_t dir_id : so_grouping) - dir_id_to_so_map[dir_id] = so_grouping_id; - - ++so_grouping_id; - } // for so_grouping - } // map scope - - return {unq_so_grps, dir_id_to_so_map}; -} - void DiscreteOrdinatesProblem::InitFluxDataStructures(LBSGroupset& groupset) { diff --git a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.cu b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.cu index faa2421b0b..c351b22992 100644 --- a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.cu +++ b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.cu @@ -5,11 +5,9 @@ #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/boundary/boundary_carrier.h" #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/boundary/sweep_boundary.h" #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/angle_set/aahd_angle_set.h" -#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/fluds/aahd_fluds_common_data.h" #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/fluds/aahd_fluds.h" #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep_chunks/aahd_sweep_chunk.h" #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/angle_set/cbcd_angle_set.h" -#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/fluds/cbcd_fluds_common_data.h" #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/fluds/cbcd_fluds.h" #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep_chunks/cbcd_sweep_chunk.h" #include "modules/linear_boltzmann_solvers/lbs_problem/device/device_vector_mirror.h" @@ -47,19 +45,6 @@ DiscreteOrdinatesProblem::ResetBoundaryCarrier() boundary_carrier_.reset(); } -void -DiscreteOrdinatesProblem::CreateAAHD_FLUDSCommonData() -{ - for (const auto& [quadrature, spds_list] : quadrature_spds_map_) - { - for (const auto& spds : spds_list) - { - quadrature_fluds_commondata_map_[quadrature].push_back( - std::make_unique(*spds, grid_nodal_mappings_, *discretization_)); - } - } -} - void DiscreteOrdinatesProblem::UpdateAAHD_FLUDSCommonDataWithBoundary() { @@ -103,19 +88,6 @@ DiscreteOrdinatesProblem::CreateAAHD_SweepChunk(LBSGroupset& groupset) return std::make_shared(*this, groupset); } -void -DiscreteOrdinatesProblem::CreateCBCD_FLUDSCommonData() -{ - for (const auto& [quadrature, spds_list] : quadrature_spds_map_) - { - for (const auto& spds : spds_list) - { - quadrature_fluds_commondata_map_[quadrature].push_back( - std::make_unique(*spds, grid_nodal_mappings_, *discretization_)); - } - } -} - std::shared_ptr DiscreteOrdinatesProblem::CreateCBCD_FLUDS(std::size_t num_groups, std::size_t num_angles, diff --git a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.h b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.h index 861dffb4a7..f129d9638c 100644 --- a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.h +++ b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/discrete_ordinates_problem.h @@ -267,7 +267,6 @@ class DiscreteOrdinatesProblem : public LBSProblem std::map, std::vector>> quadrature_fluds_commondata_map_; - std::vector verbose_sweep_angles_; const std::string sweep_type_; /** @} */ @@ -295,9 +294,6 @@ class DiscreteOrdinatesProblem : public LBSProblem std::size_t max_angleset_size_ = 0; /// Max group-set size. unsigned int max_groupset_size_ = 0; - - std::shared_ptr grid_face_histogram_ = nullptr; - std::vector> psi_new_local_; std::vector> psi_old_local_; std::optional sweep_chunk_mode_; @@ -325,7 +321,6 @@ class DiscreteOrdinatesProblem : public LBSProblem std::vector ComputeAngularFieldFunctionData(size_t groupset_id, unsigned int group, size_t angle) const; - void CreateAAHD_FLUDSCommonData(); void UpdateAAHD_FLUDSCommonDataWithBoundary(); std::shared_ptr CreateAAHD_FLUDS(unsigned int num_groups, std::size_t num_angles, @@ -341,7 +336,6 @@ class DiscreteOrdinatesProblem : public LBSProblem const MPICommunicatorSet& in_comm_set); std::shared_ptr CreateAAHD_SweepChunk(LBSGroupset& groupset); - void CreateCBCD_FLUDSCommonData(); std::shared_ptr CreateCBCD_FLUDS(std::size_t num_groups, std::size_t num_angles, std::size_t num_local_cells, @@ -360,16 +354,6 @@ class DiscreteOrdinatesProblem : public LBSProblem const MPICommunicatorSet& in_comm_set); std::shared_ptr CreateCBCDSweepChunk(LBSGroupset& groupset); - /** - * This routine groups angle-indices to groups sharing the same sweep ordering. It also takes - * geometry into account. - */ - std::pair - AssociateSOsAndDirections(std::shared_ptr grid, - const AngularQuadrature& quadrature, - AngleAggregationType agg_type, - GeometryType lbs_geo_type); - void UpdateAngularFluxStorage(); void UpdateAngularFluxFieldFunction(FieldFunctionGridBased& ff, size_t groupset_id, diff --git a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/spds/aah.h b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/spds/aah.h index d92d679572..11777c6eef 100644 --- a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/spds/aah.h +++ b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/spds/aah.h @@ -4,6 +4,7 @@ #pragma once #include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/spds/spds.h" +#include namespace opensn { @@ -51,19 +52,19 @@ class AAH_SPDS : public SPDS } /// Returns the global sweep FAS as a vector of edges. - std::vector GetGlobalSweepFAS() { return global_sweep_fas_; } + const std::vector& GetGlobalSweepFAS() const { return global_sweep_fas_; } /** * Sets the global sweep FAS. * \param edges A vector of edges representing the FAS. */ - void SetGlobalSweepFAS(std::vector& edges) { global_sweep_fas_ = edges; } + void SetGlobalSweepFAS(std::vector edges) { global_sweep_fas_ = std::move(edges); } /// Returns the locally accumulated location-to-location edge weights for this SPDS. std::vector ComputeLocalLocationEdgeWeights() const; /// Sets the global location-to-location edge weights for this SPDS. - void SetGlobalEdgeWeights(std::vector& weights) + void SetGlobalEdgeWeights(std::vector weights) { global_edge_weights_ = std::move(weights); } diff --git a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/sweep_runtime_builder.cc b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/sweep_runtime_builder.cc new file mode 100644 index 0000000000..a42d1757c7 --- /dev/null +++ b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/sweep_runtime_builder.cc @@ -0,0 +1,511 @@ +// SPDX-FileCopyrightText: 2026 The OpenSn Authors +// SPDX-License-Identifier: MIT + +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/sweep_runtime_builder.h" +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/fluds/aah_fluds_common_data.h" +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/fluds/cbc_fluds_common_data.h" +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/spds/aah.h" +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/spds/cbc.h" +#include "modules/linear_boltzmann_solvers/lbs_problem/groupset/lbs_groupset.h" +#include "framework/logging/log.h" +#include "framework/math/quadratures/angular/product_quadrature.h" +#include "framework/mesh/mesh_continuum/grid_face_histogram.h" +#include "framework/mesh/mesh_continuum/mesh_continuum.h" +#include "framework/runtime.h" +#include "framework/utils/error.h" +#include "framework/utils/timer.h" +#include +#include + +namespace opensn +{ + +#ifndef __OPENSN_WITH_GPU__ +namespace detail +{ + +void +BuildAAHGPUFludsCommonData(SweepRuntime& runtime, + const SpatialDiscretization& discretization, + const std::vector& grid_nodal_mappings) +{ + static_cast(runtime); + static_cast(discretization); + static_cast(grid_nodal_mappings); + throw std::runtime_error("BuildAAHGPUFludsCommonData: OPENSN_WITH_CUDA not enabled."); +} + +void +BuildCBCGPUFludsCommonData(SweepRuntime& runtime, + const SpatialDiscretization& discretization, + const std::vector& grid_nodal_mappings) +{ + static_cast(runtime); + static_cast(discretization); + static_cast(grid_nodal_mappings); + throw std::runtime_error("BuildCBCGPUFludsCommonData: OPENSN_WITH_CUDA not enabled."); +} + +} // namespace detail +#endif + +namespace +{ + +// Developer-only diagnostic for tiny sweep graphs. Populate with AAH_SPDS ids when debugging. +const std::vector SWEEP_ORDER_DIRECTIONS_TO_PRINT = {}; + +void +AppendNonEmptyGrouping(UniqueSOGroupings& unique_so_groupings, DirIDs dir_ids) +{ + if (not dir_ids.empty()) + unique_so_groupings.push_back(std::move(dir_ids)); +} + +DirIDToSOMap +BuildDirectionToSweepOrderingMap(const UniqueSOGroupings& unique_so_groupings) +{ + DirIDToSOMap dir_id_to_so_map; + for (size_t so_grouping_id = 0; so_grouping_id < unique_so_groupings.size(); ++so_grouping_id) + for (const size_t dir_id : unique_so_groupings[so_grouping_id]) + dir_id_to_so_map[dir_id] = so_grouping_id; + + return dir_id_to_so_map; +} + +std::pair +AssociateSOsAndDirections(const std::string& problem_name, + const std::shared_ptr& grid, + const AngularQuadrature& quadrature, + AngleAggregationType agg_type, + GeometryType geometry_type) +{ + if (quadrature.GetOmegas().empty()) + throw std::logic_error(problem_name + ": Quadrature with no omegas cannot be used"); + if (quadrature.GetWeights().empty()) + throw std::logic_error(problem_name + ": Quadrature with no weights cannot be used"); + + UniqueSOGroupings unique_so_groupings; + switch (agg_type) + { + case AngleAggregationType::SINGLE: + { + if (geometry_type == GeometryType::TWOD_CYLINDRICAL) + { + const auto* product_quad = dynamic_cast(&quadrature); + if (product_quad) + { + for (const auto& dir_set : product_quad->GetDirectionMap()) + for (const auto dir_id : dir_set.second) + AppendNonEmptyGrouping(unique_so_groupings, {dir_id}); + } + else + { + const size_t num_dirs = quadrature.GetNumAngles(); + for (size_t n = 0; n < num_dirs; ++n) + AppendNonEmptyGrouping(unique_so_groupings, {n}); + } + } + else + { + const size_t num_dirs = quadrature.GetNumAngles(); + for (size_t n = 0; n < num_dirs; ++n) + AppendNonEmptyGrouping(unique_so_groupings, {n}); + } + break; + } + case AngleAggregationType::POLAR: + { + if (grid->GetType() != ORTHOGONAL and grid->GetDimension() != 2 and not grid->Extruded()) + throw std::logic_error( + problem_name + + ": The simulation is using polar angle aggregation for which only certain geometry " + "types are supported, i.e., ORTHOGONAL, 2D or 3D EXTRUDED"); + + const auto quad_type = quadrature.GetType(); + if (quad_type != AngularQuadratureType::PRODUCT_QUADRATURE) + throw std::logic_error(problem_name + + ": The simulation is using polar angle aggregation for which only " + "Product-type quadratures are supported"); + + try + { + const auto& product_quad = dynamic_cast(quadrature); + + const auto& azimuthal_angles = product_quad.GetAzimuthalAngles(); + const auto& polar_angles = product_quad.GetPolarAngles(); + const auto num_azimuthal = azimuthal_angles.size(); + const auto num_polar = polar_angles.size(); + + std::vector upward_polar_ids; + std::vector downward_polar_ids; + for (size_t p = 0; p < num_polar; ++p) + if (polar_angles[p] > M_PI_2) + upward_polar_ids.push_back(p); + else + downward_polar_ids.push_back(p); + + auto MapPolarAndAzimuthalIDs = + [&product_quad, &unique_so_groupings](const DirIDs& polar_ids, size_t azimuthal_id) + { + DirIDs dir_ids; + dir_ids.reserve(polar_ids.size()); + for (const size_t p : polar_ids) + dir_ids.push_back(product_quad.GetAngleNum(p, azimuthal_id)); + AppendNonEmptyGrouping(unique_so_groupings, std::move(dir_ids)); + }; + + for (size_t a = 0; a < num_azimuthal; ++a) + { + if (not upward_polar_ids.empty()) + MapPolarAndAzimuthalIDs(upward_polar_ids, a); + if (not downward_polar_ids.empty()) + MapPolarAndAzimuthalIDs(downward_polar_ids, a); + } + } + catch (const std::bad_cast&) + { + throw std::runtime_error(problem_name + + ": Casting the angular quadrature to the product quadrature base " + "failed"); + } + + break; + } + case AngleAggregationType::AZIMUTHAL: + { + if (geometry_type != GeometryType::ONED_SPHERICAL and + geometry_type != GeometryType::TWOD_CYLINDRICAL) + throw std::logic_error(problem_name + + ": AZIMUTHAL aggregation is only valid for TWOD_CYLINDRICAL " + "geometry"); + + const auto quad_type = quadrature.GetType(); + if (quad_type != AngularQuadratureType::PRODUCT_QUADRATURE) + throw std::logic_error(problem_name + + ": AZIMUTHAL aggregation is only valid for TWOD_CYLINDRICAL " + "geometry."); + + try + { + const auto& product_quad = dynamic_cast(quadrature); + + for (const auto& dir_set : product_quad.GetDirectionMap()) + { + std::vector group1; + std::vector group2; + for (const auto& dir_id : dir_set.second) + if (quadrature.GetAbscissa(dir_id).phi > M_PI_2) + group1.push_back(dir_id); + else + group2.push_back(dir_id); + + AppendNonEmptyGrouping(unique_so_groupings, {group1.begin(), group1.end()}); + AppendNonEmptyGrouping(unique_so_groupings, {group2.begin(), group2.end()}); + } + } + catch (const std::bad_cast&) + { + throw std::runtime_error(problem_name + + ": Casting the angular quadrature to the product quadrature base " + "failed"); + } + + break; + } + default: + throw std::invalid_argument(problem_name + ": Called with UNDEFINED angle aggregation type"); + } + + return {unique_so_groupings, BuildDirectionToSweepOrderingMap(unique_so_groupings)}; +} + +void +BuildSweepOrderingGroups(SweepRuntime& runtime, + const std::string& problem_name, + const std::vector& groupsets, + const std::shared_ptr& grid, + GeometryType geometry_type, + std::map, bool>& allow_cycles_map) +{ + for (const auto& groupset : groupsets) + { + if (runtime.quadrature_unq_so_grouping_map.count(groupset.quadrature) == 0) + { + runtime.quadrature_unq_so_grouping_map[groupset.quadrature] = AssociateSOsAndDirections( + problem_name, grid, *groupset.quadrature, groupset.angleagg_method, geometry_type); + } + + if (allow_cycles_map.count(groupset.quadrature) == 0) + allow_cycles_map[groupset.quadrature] = groupset.allow_cycles; + } +} + +template +void +BuildSPDSForSweepOrderings(SweepRuntime& runtime, CreateSPDS create_spds) +{ + for (const auto& [quadrature, info] : runtime.quadrature_unq_so_grouping_map) + { + auto& spds_list = runtime.quadrature_spds_map[quadrature]; + const auto& unique_so_groupings = info.first; + for (const auto& so_grouping : unique_so_groupings) + { + const size_t master_dir_id = so_grouping.front(); + const auto& omega = quadrature->GetOmega(master_dir_id); + spds_list.push_back(create_spds(quadrature, spds_list.size(), omega)); + } + } +} + +void +BuildAAHSPDS(SweepRuntime& runtime, + const std::shared_ptr& grid, + const std::map, bool>& allow_cycles_map, + bool use_gpus) +{ + log.Log0Verbose1() << program_timer.GetTimeString() << " Initializing AAH SPDS."; + BuildSPDSForSweepOrderings( + runtime, + [&grid, &allow_cycles_map, use_gpus](const std::shared_ptr& quadrature, + size_t sweep_ordering_id, + const Vector3& omega) -> std::shared_ptr + { + return std::make_shared(static_cast(sweep_ordering_id), + omega, + grid, + allow_cycles_map.at(quadrature), + use_gpus); + }); +} + +void +BuildCBCSPDS(SweepRuntime& runtime, + const std::shared_ptr& grid, + const std::map, bool>& allow_cycles_map) +{ + BuildSPDSForSweepOrderings( + runtime, + [&grid, &allow_cycles_map](const std::shared_ptr& quadrature, + size_t, + const Vector3& omega) -> std::shared_ptr + { return std::make_shared(omega, grid, allow_cycles_map.at(quadrature)); }); +} + +std::vector> +GetAAHSPDSList(const SweepRuntime& runtime) +{ + std::vector> aah_spds_list; + // The flattened order is the MPI serialization key for gathered AAH sweep FAS records. Keep this + // traversal deterministic across ranks and in sync with ApplyAAHSweepFAS. + for (const auto& [quadrature, spds_list] : runtime.quadrature_spds_map) + for (const auto& spds : spds_list) + aah_spds_list.push_back(std::static_pointer_cast(spds)); + + return aah_spds_list; +} + +int +GetSweepGraphOwner(size_t spds_ordinal) +{ + return static_cast(spds_ordinal % static_cast(opensn::mpi_comm.size())); +} + +void +AccumulateAAHGlobalEdgeWeights(const std::vector>& spds_list) +{ + const int comm_size = opensn::mpi_comm.size(); + const int matrix_size = comm_size * comm_size; + std::vector recv_counts(comm_size, comm_size); + std::vector recv_displacements(comm_size, 0); + for (int loc = 0; loc < comm_size; ++loc) + recv_displacements[loc] = loc * comm_size; + + for (size_t spds_ordinal = 0; spds_ordinal < spds_list.size(); ++spds_ordinal) + { + const auto& spds = spds_list[spds_ordinal]; + const int owner = GetSweepGraphOwner(spds_ordinal); + + const auto local_row = spds->ComputeLocalLocationEdgeWeights(); + std::vector recv; + if (opensn::mpi_comm.rank() == owner) + recv.assign(matrix_size, 0.0); + opensn::mpi_comm.gather(local_row, recv, recv_counts, recv_displacements, owner); + + if (opensn::mpi_comm.rank() == owner) + spds->SetGlobalEdgeWeights(std::move(recv)); + } +} + +void +BuildOwnedAAHSweepFAS(const std::vector>& spds_list) +{ + log.Log0Verbose1() << program_timer.GetTimeString() << " Build global sweep FAS for each SPDS."; + for (size_t spds_ordinal = 0; spds_ordinal < spds_list.size(); ++spds_ordinal) + if (opensn::mpi_comm.rank() == GetSweepGraphOwner(spds_ordinal)) + spds_list[spds_ordinal]->BuildGlobalSweepFAS(); +} + +std::vector +GatherAAHSweepFAS(const std::vector>& spds_list) +{ + log.Log0Verbose1() << program_timer.GetTimeString() << " Gather FAS for each SPDS."; + std::vector local_edges_to_remove; + for (size_t spds_ordinal = 0; spds_ordinal < spds_list.size(); ++spds_ordinal) + { + if (opensn::mpi_comm.rank() == GetSweepGraphOwner(spds_ordinal)) + { + const auto& edges_to_remove = spds_list[spds_ordinal]->GetGlobalSweepFAS(); + local_edges_to_remove.push_back(static_cast(spds_ordinal)); + local_edges_to_remove.push_back(static_cast(edges_to_remove.size())); + local_edges_to_remove.insert( + local_edges_to_remove.end(), edges_to_remove.begin(), edges_to_remove.end()); + } + } + + int local_size = static_cast(local_edges_to_remove.size()); + std::vector receive_counts(opensn::mpi_comm.size(), 0); + std::vector displacements(opensn::mpi_comm.size(), 0); + mpi_comm.all_gather(local_size, receive_counts); + + int total_size = 0; + for (size_t i = 0; i < receive_counts.size(); ++i) + { + displacements[i] = total_size; + total_size += receive_counts[i]; + } + + std::vector global_edges_to_remove(total_size, 0); + mpi_comm.all_gather(local_edges_to_remove, global_edges_to_remove, receive_counts, displacements); + + return global_edges_to_remove; +} + +void +ApplyAAHSweepFAS(const std::vector>& spds_list, + const std::vector& global_edges_to_remove) +{ + size_t offset = 0; + while (offset < global_edges_to_remove.size()) + { + if (offset + 2 > global_edges_to_remove.size()) + OpenSnLogicalError("Malformed AAH sweep FAS payload."); + + const auto spds_ordinal = static_cast(global_edges_to_remove[offset++]); + const auto num_edges = global_edges_to_remove[offset++]; + if (num_edges < 0 or offset + static_cast(num_edges) > global_edges_to_remove.size()) + OpenSnLogicalError("Malformed AAH sweep FAS payload."); + + std::vector edges; + edges.reserve(num_edges); + for (int i = 0; i < num_edges; ++i) + edges.emplace_back(global_edges_to_remove[offset++]); + + if (spds_ordinal >= spds_list.size()) + OpenSnLogicalError("Invalid SPDS ordinal gathered while building AAH sweep graph."); + spds_list[spds_ordinal]->SetGlobalSweepFAS(std::move(edges)); + } +} + +void +BuildAAHSweepTDGs(const std::vector>& spds_list) +{ + log.Log0Verbose1() << program_timer.GetTimeString() << " Build global sweep TDGs."; + for (const auto& spds : spds_list) + spds->BuildGlobalSweepTDG(); +} + +void +PrintRequestedSweepGraphs(const std::vector>& spds_list) +{ + for (const auto& spds : spds_list) + for (const int dir_id : SWEEP_ORDER_DIRECTIONS_TO_PRINT) + if (spds->GetId() == dir_id) + spds->PrintGhostedGraph(); +} + +void +BuildAAHGlobalSweepGraph(SweepRuntime& runtime) +{ + // Creating an AAH SPDS can be expensive, so the global graph work is split across MPI ranks: + // accumulate graph edge weights, build the feedback arc set on owning ranks, gather the FAS on + // all ranks, then build each global sweep task dependency graph locally. + const auto spds_list = GetAAHSPDSList(runtime); + AccumulateAAHGlobalEdgeWeights(spds_list); + BuildOwnedAAHSweepFAS(spds_list); + ApplyAAHSweepFAS(spds_list, GatherAAHSweepFAS(spds_list)); + BuildAAHSweepTDGs(spds_list); + PrintRequestedSweepGraphs(spds_list); +} + +void +BuildAAHCPUFludsCommonData(SweepRuntime& runtime, + const std::shared_ptr& grid, + const std::vector& grid_nodal_mappings) +{ + const auto grid_face_histogram = grid->MakeGridFaceHistogram(); + for (const auto& [quadrature, spds_list] : runtime.quadrature_spds_map) + for (const auto& spds : spds_list) + runtime.quadrature_fluds_commondata_map[quadrature].push_back( + std::make_unique(grid_nodal_mappings, *spds, *grid_face_histogram)); +} + +void +BuildCBCCPUFludsCommonData(SweepRuntime& runtime, + const std::vector& grid_nodal_mappings) +{ + for (const auto& [quadrature, spds_list] : runtime.quadrature_spds_map) + for (const auto& spds : spds_list) + runtime.quadrature_fluds_commondata_map[quadrature].push_back( + std::make_unique(*spds, grid_nodal_mappings)); +} + +} // namespace + +SweepRuntime +BuildSweepRuntime(const std::string& problem_name, + const std::vector& groupsets, + const std::shared_ptr& grid, + const std::string& sweep_type, + bool use_gpus, + const SpatialDiscretization& discretization, + const std::vector& grid_nodal_mappings) +{ + SweepRuntime runtime; + std::map, bool> quadrature_allow_cycles_map; + + const auto geometry_type = grid->GetGeometryType(); + BuildSweepOrderingGroups( + runtime, problem_name, groupsets, grid, geometry_type, quadrature_allow_cycles_map); + + if (sweep_type == "AAH") + { + BuildAAHSPDS(runtime, grid, quadrature_allow_cycles_map, use_gpus); + BuildAAHGlobalSweepGraph(runtime); + } + else if (sweep_type == "CBC") + BuildCBCSPDS(runtime, grid, quadrature_allow_cycles_map); + else + OpenSnInvalidArgument("Unsupported sweep type \"" + sweep_type + "\""); + + opensn::mpi_comm.barrier(); + + if (not use_gpus) + { + if (sweep_type == "AAH") + BuildAAHCPUFludsCommonData(runtime, grid, grid_nodal_mappings); + else if (sweep_type == "CBC") + BuildCBCCPUFludsCommonData(runtime, grid_nodal_mappings); + } + else + { + if (sweep_type == "AAH") + detail::BuildAAHGPUFludsCommonData(runtime, discretization, grid_nodal_mappings); + else if (sweep_type == "CBC") + detail::BuildCBCGPUFludsCommonData(runtime, discretization, grid_nodal_mappings); + } + + return runtime; +} + +} // namespace opensn diff --git a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/sweep_runtime_builder.cu b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/sweep_runtime_builder.cu new file mode 100644 index 0000000000..09158fb967 --- /dev/null +++ b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/sweep_runtime_builder.cu @@ -0,0 +1,34 @@ +// SPDX-FileCopyrightText: 2026 The OpenSn Authors +// SPDX-License-Identifier: MIT + +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/sweep_runtime_builder.h" + +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/fluds/aahd_fluds_common_data.h" +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/fluds/cbcd_fluds_common_data.h" + +namespace opensn::detail +{ + +void +BuildAAHGPUFludsCommonData(SweepRuntime& runtime, + const SpatialDiscretization& discretization, + const std::vector& grid_nodal_mappings) +{ + for (const auto& [quadrature, spds_list] : runtime.quadrature_spds_map) + for (const auto& spds : spds_list) + runtime.quadrature_fluds_commondata_map[quadrature].push_back( + std::make_unique(*spds, grid_nodal_mappings, discretization)); +} + +void +BuildCBCGPUFludsCommonData(SweepRuntime& runtime, + const SpatialDiscretization& discretization, + const std::vector& grid_nodal_mappings) +{ + for (const auto& [quadrature, spds_list] : runtime.quadrature_spds_map) + for (const auto& spds : spds_list) + runtime.quadrature_fluds_commondata_map[quadrature].push_back( + std::make_unique(*spds, grid_nodal_mappings, discretization)); +} + +} // namespace opensn::detail diff --git a/modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/sweep_runtime_builder.h b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/sweep_runtime_builder.h new file mode 100644 index 0000000000..ce17780a5b --- /dev/null +++ b/modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/sweep_runtime_builder.h @@ -0,0 +1,55 @@ +// SPDX-FileCopyrightText: 2026 The OpenSn Authors +// SPDX-License-Identifier: MIT + +#pragma once + +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/angle_aggregation/angle_aggregation.h" +#include "modules/linear_boltzmann_solvers/discrete_ordinates_problem/sweep/fluds/fluds_common_data.h" +#include +#include +#include +#include + +namespace opensn +{ + +class GridFaceHistogram; +class LBSGroupset; +class MeshContinuum; +class SPDS; +class SpatialDiscretization; + +using SweepOrderGroupingInfo = std::pair; + +struct SweepRuntime +{ + std::map, SweepOrderGroupingInfo> + quadrature_unq_so_grouping_map; + std::map, std::vector>> + quadrature_spds_map; + std::map, std::vector>> + quadrature_fluds_commondata_map; +}; + +SweepRuntime BuildSweepRuntime(const std::string& problem_name, + const std::vector& groupsets, + const std::shared_ptr& grid, + const std::string& sweep_type, + bool use_gpus, + const SpatialDiscretization& discretization, + const std::vector& grid_nodal_mappings); + +namespace detail +{ + +void BuildAAHGPUFludsCommonData(SweepRuntime& runtime, + const SpatialDiscretization& discretization, + const std::vector& grid_nodal_mappings); + +void BuildCBCGPUFludsCommonData(SweepRuntime& runtime, + const SpatialDiscretization& discretization, + const std::vector& grid_nodal_mappings); + +} // namespace detail + +} // namespace opensn