Line data Source code
1 0 : // Distributed under the MIT License. 2 : // See LICENSE.txt for details. 3 : 4 : #pragma once 5 : 6 : #include <algorithm> 7 : #include <charm++.h> 8 : #include <cstddef> 9 : #include <optional> 10 : #include <unordered_map> 11 : #include <unordered_set> 12 : #include <vector> 13 : 14 : #include "Domain/Creators/Tags/Domain.hpp" 15 : #include "Domain/Creators/Tags/InitialExtents.hpp" 16 : #include "Domain/Creators/Tags/InitialRefinementLevels.hpp" 17 : #include "Domain/ElementDistribution.hpp" 18 : #include "Domain/Structure/ElementId.hpp" 19 : #include "Domain/Tags/ElementDistribution.hpp" 20 : #include "Elliptic/DiscontinuousGalerkin/Tags.hpp" 21 : #include "NumericalAlgorithms/Convergence/Tags.hpp" 22 : #include "Parallel/GlobalCache.hpp" 23 : #include "Parallel/Local.hpp" 24 : #include "Parallel/Printf.hpp" 25 : #include "Parallel/Protocols/ArrayElementsAllocator.hpp" 26 : #include "Parallel/Section.hpp" 27 : #include "Parallel/Tags/Section.hpp" 28 : #include "ParallelAlgorithms/LinearSolver/Multigrid/Hierarchy.hpp" 29 : #include "ParallelAlgorithms/LinearSolver/Multigrid/Tags.hpp" 30 : #include "Utilities/ErrorHandling/Error.hpp" 31 : #include "Utilities/PrettyType.hpp" 32 : #include "Utilities/ProtocolHelpers.hpp" 33 : #include "Utilities/Serialization/Serialize.hpp" 34 : #include "Utilities/System/ParallelInfo.hpp" 35 : #include "Utilities/TaggedTuple.hpp" 36 : 37 : namespace LinearSolver::multigrid { 38 : 39 : /*! 40 : * \brief A `Parallel::protocols::ArrayElementsAllocator` that creates array 41 : * elements to cover the initial computational domain multiple times at 42 : * different refinement levels, suitable for the 43 : * `LinearSolver::multigrid::Multigrid` algorithm 44 : * 45 : * The `initial_element_ids` function and the option-created 46 : * `domain::Tags::Domain` and `domain::Tags::InitialRefinementLevels` determine 47 : * the initial set of element IDs in the domain. This is taken as the finest 48 : * grid in the multigrid hierarchy. Coarser grids are determined by successively 49 : * applying `LinearSolver::multigrid::coarsen`, up to 50 : * `LinearSolver::multigrid::Tags::MaxLevels` grids. 51 : * 52 : * Array elements are created for all element IDs on all grids, meaning they all 53 : * share the same parallel component, action list etc. Elements are connected to 54 : * their neighbors _on the same grid_ by the 55 : * `domain::Initialization::create_initial_element` function (this is 56 : * independent of the multigrid code, the function is typically called in an 57 : * initialization action). Elements are connected to their parent and children 58 : * _across grids_ by the multigrid-tags set here and in the multigrid 59 : * initialization actions (see `LinearSolver::multigrid::parent_id` and 60 : * `LinearSolver::multigrid::child_ids`). 61 : * 62 : * This allocator also creates two sets of sections (see `Parallel::Section`): 63 : * - `Parallel::Tags::Section<ElementArray, 64 : * LinearSolver::multigrid::Tags::MultigridLevel>`: One section per grid. All 65 : * elements are part of a section with this tag. 66 : * - `Parallel::Tags::Section<ElementArray, 67 : * LinearSolver::multigrid::Tags::IsFinestGrid>`: A single section that holds 68 : * only the elements on the finest grid. Holds `std::nullopt` on all other 69 : * elements. 70 : * 71 : * The elements are distributed on processors using the 72 : * `domain::BlockZCurveProcDistribution` for every grid independently. An 73 : * unordered set of `size_t`s can be passed to the `apply` function which 74 : * represents physical processors to avoid placing elements on. 75 : */ 76 : template <size_t Dim, typename OptionsGroup> 77 1 : struct ElementsAllocator 78 : : tt::ConformsTo<Parallel::protocols::ArrayElementsAllocator> { 79 : template <typename ElementArray> 80 0 : using array_allocation_tags = 81 : tmpl::list<Parallel::Tags::Section<ElementArray, Tags::MultigridLevel>, 82 : Parallel::Tags::Section<ElementArray, Tags::IsFinestGrid>>; 83 : 84 : template <typename ElementArray, typename Metavariables, 85 : typename... InitializationTags> 86 0 : static void apply(Parallel::CProxy_GlobalCache<Metavariables>& global_cache, 87 : const tuples::TaggedTuple<InitializationTags...>& 88 : original_initialization_items, 89 : const tuples::tagged_tuple_from_typelist< 90 : typename ElementArray:: 91 : array_allocation_tags>& /*array_allocation_items*/ 92 : = {}, 93 : const std::unordered_set<size_t>& procs_to_ignore = {}) { 94 : // Copy the initialization items so we can adjust them on each refinement 95 : // level 96 : auto initialization_items = 97 : deserialize<tuples::TaggedTuple<InitializationTags...>>( 98 : serialize<tuples::TaggedTuple<InitializationTags...>>( 99 : original_initialization_items) 100 : .data()); 101 : auto& local_cache = *Parallel::local_branch(global_cache); 102 : auto& element_array = 103 : Parallel::get_parallel_component<ElementArray>(local_cache); 104 : const auto& domain = get<domain::Tags::Domain<Dim>>(local_cache); 105 : const auto& initial_extents = 106 : get<domain::Tags::InitialExtents<Dim>>(initialization_items); 107 : auto& initial_refinement_levels = 108 : get<domain::Tags::InitialRefinementLevels<Dim>>(initialization_items); 109 : auto& children_refinement_levels = 110 : get<Tags::ChildrenRefinementLevels<Dim>>(initialization_items); 111 : auto& parent_refinement_levels = 112 : get<Tags::ParentRefinementLevels<Dim>>(initialization_items); 113 : const auto& quadrature = 114 : Parallel::get<elliptic::dg::Tags::Quadrature>(local_cache); 115 : const std::optional<domain::ElementWeight>& element_weight = 116 : get<domain::Tags::ElementDistribution>(local_cache); 117 : std::optional<size_t> max_levels = 118 : get<Tags::MaxLevels<OptionsGroup>>(local_cache); 119 : const size_t number_of_procs = 120 : Parallel::number_of_procs<size_t>(local_cache); 121 : if (max_levels == 0) { 122 : ERROR_NO_TRACE( 123 : "The 'MaxLevels' option includes the finest grid, so '0' is not a " 124 : "valid value. Set the option to '1' to effectively disable " 125 : "multigrid."); 126 : } 127 : const size_t num_iterations = 128 : get<Convergence::Tags::Iterations<OptionsGroup>>(local_cache); 129 : if (UNLIKELY(num_iterations == 0)) { 130 : Parallel::printf( 131 : "%s is disabled (zero iterations). Only creating a single grid.\n", 132 : pretty_type::name<OptionsGroup>()); 133 : max_levels = 1; 134 : } 135 : const auto& blocks = domain.blocks(); 136 : size_t multigrid_level = 0; 137 : do { 138 : // Store the current grid as child grid before coarsening it 139 : children_refinement_levels = initial_refinement_levels; 140 : initial_refinement_levels = parent_refinement_levels; 141 : // Construct coarsened (parent) grid 142 : if (not max_levels.has_value() or multigrid_level < *max_levels - 1) { 143 : parent_refinement_levels = 144 : LinearSolver::multigrid::coarsen(initial_refinement_levels); 145 : } 146 : // Create element IDs for all elements on this level 147 : std::vector<ElementId<Dim>> element_ids{}; 148 : for (const auto& block : blocks) { 149 : const std::vector<ElementId<Dim>> block_element_ids = 150 : initial_element_ids(block.id(), 151 : initial_refinement_levels[block.id()], 152 : multigrid_level); 153 : element_ids.insert(element_ids.begin(), block_element_ids.begin(), 154 : block_element_ids.end()); 155 : } 156 : // Create an array section for this refinement level 157 : std::vector<CkArrayIndex> array_indices(element_ids.size()); 158 : std::transform( 159 : element_ids.begin(), element_ids.end(), array_indices.begin(), 160 : [](const ElementId<Dim>& local_element_id) { 161 : return Parallel::ArrayIndex<ElementId<Dim>>(local_element_id); 162 : }); 163 : using MultigridLevelSection = 164 : Parallel::Section<ElementArray, Tags::MultigridLevel>; 165 : get<Parallel::Tags::Section<ElementArray, Tags::MultigridLevel>>( 166 : initialization_items) = MultigridLevelSection{ 167 : multigrid_level, MultigridLevelSection::cproxy_section::ckNew( 168 : element_array.ckGetArrayID(), 169 : array_indices.data(), array_indices.size())}; 170 : using FinestGridSection = 171 : Parallel::Section<ElementArray, Tags::IsFinestGrid>; 172 : get<Parallel::Tags::Section<ElementArray, Tags::IsFinestGrid>>( 173 : initialization_items) = 174 : multigrid_level == 0 175 : ? std::make_optional(FinestGridSection{ 176 : true, FinestGridSection::cproxy_section::ckNew( 177 : element_array.ckGetArrayID(), 178 : array_indices.data(), array_indices.size())}) 179 : : std::nullopt; 180 : // Create the elements for this refinement level and distribute them among 181 : // processors 182 : const size_t num_of_procs_to_use = 183 : static_cast<size_t>(sys::number_of_procs()) - procs_to_ignore.size(); 184 : // Distributed with weighted space filling curve 185 : if (element_weight.has_value()) { 186 : const std::unordered_map<ElementId<Dim>, double> element_costs = 187 : domain::get_element_costs(blocks, initial_refinement_levels, 188 : initial_extents, element_weight.value(), 189 : quadrature); 190 : const domain::BlockZCurveProcDistribution<Dim> element_distribution{ 191 : element_costs, num_of_procs_to_use, 192 : blocks, initial_refinement_levels, 193 : initial_extents, procs_to_ignore}; 194 : 195 : for (const auto& element_id : element_ids) { 196 : const size_t target_proc = 197 : element_distribution.get_proc_for_element(element_id); 198 : element_array(element_id) 199 : .insert(global_cache, initialization_items, target_proc); 200 : } 201 : } else { 202 : // Distributed with round-robin 203 : size_t which_proc = 0; 204 : for (const auto& element_id : element_ids) { 205 : while (procs_to_ignore.find(which_proc) != procs_to_ignore.end()) { 206 : which_proc = which_proc + 1 == number_of_procs ? 0 : which_proc + 1; 207 : } 208 : 209 : element_array(element_id) 210 : .insert(global_cache, initialization_items, which_proc); 211 : 212 : which_proc = which_proc + 1 == number_of_procs ? 0 : which_proc + 1; 213 : } 214 : } 215 : Parallel::printf( 216 : "%s level %zu has %zu elements in %zu blocks distributed on %d " 217 : "procs.\n", 218 : pretty_type::name<OptionsGroup>(), multigrid_level, 219 : element_ids.size(), blocks.size(), num_of_procs_to_use); 220 : ++multigrid_level; 221 : } while (initial_refinement_levels != parent_refinement_levels); 222 : element_array.doneInserting(); 223 : } 224 : }; 225 : 226 : } // namespace LinearSolver::multigrid