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 "DataStructures/TaggedTuple.hpp"
15 : #include "Domain/Creators/Tags/Domain.hpp"
16 : #include "Domain/Creators/Tags/InitialExtents.hpp"
17 : #include "Domain/Creators/Tags/InitialRefinementLevels.hpp"
18 : #include "Domain/ElementDistribution.hpp"
19 : #include "Domain/Structure/ElementId.hpp"
20 : #include "Domain/Tags/ElementDistribution.hpp"
21 : #include "Elliptic/DiscontinuousGalerkin/Tags.hpp"
22 : #include "NumericalAlgorithms/Convergence/Tags.hpp"
23 : #include "NumericalAlgorithms/Spectral/Basis.hpp"
24 : #include "Parallel/GlobalCache.hpp"
25 : #include "Parallel/Local.hpp"
26 : #include "Parallel/Printf/Printf.hpp"
27 : #include "Parallel/Protocols/ArrayElementsAllocator.hpp"
28 : #include "Parallel/Section.hpp"
29 : #include "Parallel/Tags/Section.hpp"
30 : #include "ParallelAlgorithms/LinearSolver/Multigrid/Hierarchy.hpp"
31 : #include "ParallelAlgorithms/LinearSolver/Multigrid/Tags.hpp"
32 : #include "Utilities/ErrorHandling/Error.hpp"
33 : #include "Utilities/PrettyType.hpp"
34 : #include "Utilities/ProtocolHelpers.hpp"
35 : #include "Utilities/Serialization/Serialize.hpp"
36 : #include "Utilities/System/ParallelInfo.hpp"
37 :
38 : namespace LinearSolver::multigrid {
39 :
40 : /*!
41 : * \brief A `Parallel::protocols::ArrayElementsAllocator` that creates array
42 : * elements to cover the initial computational domain multiple times at
43 : * different refinement levels, suitable for the
44 : * `LinearSolver::multigrid::Multigrid` algorithm
45 : *
46 : * The `initial_element_ids` function and the option-created
47 : * `domain::Tags::Domain` and `domain::Tags::InitialRefinementLevels` determine
48 : * the initial set of element IDs in the domain. This is taken as the finest
49 : * grid in the multigrid hierarchy. Coarser grids are determined by successively
50 : * applying `LinearSolver::multigrid::coarsen`, up to
51 : * `LinearSolver::multigrid::Tags::InitialCoarseLevels` coarser grids.
52 : *
53 : * Array elements are created for all element IDs on all grids, meaning they all
54 : * share the same parallel component, action list etc. Elements are connected to
55 : * their neighbors _on the same grid_ by the
56 : * `domain::create_initial_element` function (this is independent of the
57 : * multigrid code, the function is typically called in an initialization
58 : * action). Elements are connected to their parent and children _across grids_
59 : * by the multigrid-tags set here and in the multigrid initialization actions
60 : * (see `LinearSolver::multigrid::parent_id` and
61 : * `LinearSolver::multigrid::child_ids`).
62 : *
63 : * Once the multigrid hierarchy is created, the AMR infrastructure registers
64 : * the grids and creates the following sections (see `Parallel::Section`).
65 : * - `Parallel::Tags::Section<ElementArray, amr::Tags::GridIndex>`: One section
66 : * per grid. All elements are part of a section with this tag.
67 : * - `Parallel::Tags::Section<ElementArray, amr::Tags::IsFinestGrid>`: A single
68 : * section that holds only the elements on the finest grid. Holds
69 : * `std::nullopt` on all other 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 = tmpl::list<>;
81 :
82 : template <typename ElementArray, typename Metavariables,
83 : typename... InitializationTags>
84 0 : static void apply(Parallel::CProxy_GlobalCache<Metavariables>& global_cache,
85 : const tuples::TaggedTuple<InitializationTags...>&
86 : original_initialization_items,
87 : const tuples::tagged_tuple_from_typelist<
88 : typename ElementArray::
89 : array_allocation_tags>& /*array_allocation_items*/
90 : = {},
91 : const std::unordered_set<size_t>& procs_to_ignore = {}) {
92 : // Copy the initialization items so we can adjust them on each refinement
93 : // level
94 : auto initialization_items =
95 : serialize_and_deserialize<tuples::TaggedTuple<InitializationTags...>>(
96 : original_initialization_items);
97 : auto& local_cache = *Parallel::local_branch(global_cache);
98 : auto& element_array =
99 : Parallel::get_parallel_component<ElementArray>(local_cache);
100 : const auto& domain = get<domain::Tags::Domain<Dim>>(local_cache);
101 : const auto& initial_extents =
102 : get<domain::Tags::InitialExtents<Dim>>(initialization_items);
103 : auto& initial_refinement_levels =
104 : get<domain::Tags::InitialRefinementLevels<Dim>>(initialization_items);
105 : auto& children_refinement_levels =
106 : get<Tags::ChildrenRefinementLevels<Dim>>(initialization_items);
107 : auto& parent_refinement_levels =
108 : get<Tags::ParentRefinementLevels<Dim>>(initialization_items);
109 : const auto basis = Spectral::Basis::Legendre;
110 : const auto& quadrature =
111 : Parallel::get<elliptic::dg::Tags::Quadrature>(local_cache);
112 : const std::optional<domain::ElementWeight>& element_weight =
113 : get<domain::Tags::ElementDistribution>(local_cache);
114 : std::optional<size_t> max_coarse_levels =
115 : get<Tags::InitialCoarseLevels<OptionsGroup>>(local_cache);
116 : const size_t number_of_procs =
117 : Parallel::number_of_procs<size_t>(local_cache);
118 : const size_t num_iterations =
119 : get<Convergence::Tags::Iterations<OptionsGroup>>(local_cache);
120 : if (UNLIKELY(num_iterations == 0)) {
121 : Parallel::printf(
122 : "%s is disabled (zero iterations). Only creating a single grid.\n",
123 : pretty_type::name<OptionsGroup>());
124 : max_coarse_levels = 0;
125 : }
126 : const auto& blocks = domain.blocks();
127 :
128 : // Determine number of initial multigrid levels
129 : size_t num_levels = 1;
130 : for (const auto& ref_levs : initial_refinement_levels) {
131 : num_levels = std::max(
132 : num_levels, *std::max_element(ref_levs.begin(), ref_levs.end()) + 1);
133 : }
134 : if (max_coarse_levels.has_value()) {
135 : num_levels = std::min(num_levels, max_coarse_levels.value() + 1);
136 : }
137 :
138 : // Create the initial grids
139 : for (int multigrid_level = static_cast<int>(num_levels) - 1;
140 : multigrid_level >= 0; --multigrid_level) {
141 : // Store the current grid as child grid before coarsening it
142 : children_refinement_levels = initial_refinement_levels;
143 : initial_refinement_levels = parent_refinement_levels;
144 : // Construct coarsened (parent) grid
145 : if (multigrid_level > 0) {
146 : parent_refinement_levels =
147 : LinearSolver::multigrid::coarsen(initial_refinement_levels);
148 : }
149 : // Create element IDs for all elements on this level
150 : std::vector<ElementId<Dim>> element_ids{};
151 : for (const auto& block : blocks) {
152 : const std::vector<ElementId<Dim>> block_element_ids =
153 : initial_element_ids(block.id(),
154 : initial_refinement_levels[block.id()],
155 : static_cast<size_t>(multigrid_level));
156 : element_ids.insert(element_ids.begin(), block_element_ids.begin(),
157 : block_element_ids.end());
158 : }
159 : // Create the elements for this refinement level and distribute them among
160 : // processors
161 : const size_t num_of_procs_to_use =
162 : static_cast<size_t>(sys::number_of_procs()) - procs_to_ignore.size();
163 : // Distributed with weighted space filling curve
164 : if (element_weight.has_value()) {
165 : const std::unordered_map<ElementId<Dim>, double> element_costs =
166 : domain::get_element_costs(blocks, initial_refinement_levels,
167 : initial_extents, element_weight.value(),
168 : basis, quadrature);
169 : const domain::BlockZCurveProcDistribution<Dim> element_distribution{
170 : element_costs, num_of_procs_to_use, blocks,
171 : initial_refinement_levels, procs_to_ignore};
172 :
173 : for (const auto& element_id : element_ids) {
174 : const size_t target_proc =
175 : element_distribution.get_proc_for_element(element_id);
176 : element_array(element_id)
177 : .insert(global_cache, initialization_items, target_proc);
178 : }
179 : } else {
180 : // Distributed with round-robin
181 : size_t which_proc = 0;
182 : for (const auto& element_id : element_ids) {
183 : while (procs_to_ignore.find(which_proc) != procs_to_ignore.end()) {
184 : which_proc = which_proc + 1 == number_of_procs ? 0 : which_proc + 1;
185 : }
186 :
187 : element_array(element_id)
188 : .insert(global_cache, initialization_items, which_proc);
189 :
190 : which_proc = which_proc + 1 == number_of_procs ? 0 : which_proc + 1;
191 : }
192 : }
193 : }
194 : element_array.doneInserting();
195 : }
196 : };
197 :
198 : } // namespace LinearSolver::multigrid
|
