Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <map>
8 : #include <optional>
9 : #include <unordered_map>
10 : #include <unordered_set>
11 :
12 : #include "DataStructures/ApplyMatrices.hpp"
13 : #include "DataStructures/DataBox/DataBox.hpp"
14 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
15 : #include "DataStructures/FixedHashMap.hpp"
16 : #include "DataStructures/Matrix.hpp"
17 : #include "Domain/Creators/Tags/InitialRefinementLevels.hpp"
18 : #include "Domain/Structure/ChildSize.hpp"
19 : #include "Domain/Structure/ElementId.hpp"
20 : #include "IO/Logging/Tags.hpp"
21 : #include "IO/Observer/Tags.hpp"
22 : #include "NumericalAlgorithms/Convergence/Tags.hpp"
23 : #include "NumericalAlgorithms/Spectral/Projection.hpp"
24 : #include "Parallel/AlgorithmExecution.hpp"
25 : #include "Parallel/GlobalCache.hpp"
26 : #include "Parallel/InboxInserters.hpp"
27 : #include "Parallel/Invoke.hpp"
28 : #include "Parallel/Printf/Printf.hpp"
29 : #include "ParallelAlgorithms/Actions/Goto.hpp"
30 : #include "ParallelAlgorithms/Amr/Protocols/Projector.hpp"
31 : #include "ParallelAlgorithms/Amr/Tags.hpp"
32 : #include "ParallelAlgorithms/LinearSolver/Multigrid/Actions/RestrictFields.hpp"
33 : #include "ParallelAlgorithms/LinearSolver/Multigrid/Hierarchy.hpp"
34 : #include "ParallelAlgorithms/LinearSolver/Multigrid/Tags.hpp"
35 : #include "ParallelAlgorithms/LinearSolver/Tags.hpp"
36 : #include "Utilities/ConstantExpressions.hpp"
37 : #include "Utilities/ErrorHandling/Assert.hpp"
38 : #include "Utilities/GetOutput.hpp"
39 : #include "Utilities/PrettyType.hpp"
40 : #include "Utilities/ProtocolHelpers.hpp"
41 : #include "Utilities/TaggedTuple.hpp"
42 : #include "Utilities/TypeTraits/IsA.hpp"
43 :
44 : namespace LinearSolver::multigrid::detail {
45 :
46 : /// \cond
47 : template <typename FieldsTag, typename OptionsGroup, typename SourceTag>
48 : struct SendCorrectionToFinerGrid;
49 : template <typename FieldsTag, typename OptionsGroup, typename SourceTag>
50 : struct SkipBottomSolver;
51 : template <typename FieldsTag, typename OptionsGroup, typename SourceTag>
52 : struct SkipPostSmoothingAtBottom;
53 : /// \endcond
54 :
55 : struct PostSmoothingBeginLabel {};
56 :
57 : template <size_t Dim, typename FieldsTag, typename OptionsGroup,
58 : typename SourceTag>
59 : struct InitializeElement : tt::ConformsTo<amr::protocols::Projector> {
60 : private:
61 : using VolumeDataVars =
62 : typename Tags::VolumeDataForOutput<OptionsGroup, FieldsTag>::type;
63 :
64 : public: // Iterable action
65 : using simple_tags_from_options =
66 : tmpl::list<Tags::ChildrenRefinementLevels<Dim>,
67 : Tags::ParentRefinementLevels<Dim>>;
68 : using simple_tags =
69 : tmpl::list<amr::Tags::ParentId<Dim>, amr::Tags::ChildIds<Dim>,
70 : amr::Tags::ParentMesh<Dim>,
71 : LinearSolver::Tags::ObservationId<OptionsGroup>,
72 : Tags::VolumeDataForOutput<OptionsGroup, FieldsTag>>;
73 : using compute_tags = tmpl::list<>;
74 : using const_global_cache_tags =
75 : tmpl::list<Tags::InitialCoarseLevels<OptionsGroup>,
76 : LinearSolver::Tags::OutputVolumeData<OptionsGroup>>;
77 :
78 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
79 : typename ActionList, typename ParallelComponent>
80 : static Parallel::iterable_action_return_t apply(
81 : db::DataBox<DbTagsList>& box,
82 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
83 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
84 : const ElementId<Dim>& /*element_id*/, const ActionList /*meta*/,
85 : const ParallelComponent* const /*meta*/) {
86 : db::mutate_apply<InitializeElement>(make_not_null(&box));
87 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
88 : }
89 :
90 : public: // amr::protocols::Projector
91 : using argument_tags =
92 : tmpl::list<domain::Tags::Mesh<Dim>, domain::Tags::Element<Dim>,
93 : domain::Tags::InitialRefinementLevels<Dim>,
94 : LinearSolver::Tags::OutputVolumeData<OptionsGroup>>;
95 : using return_tags = tmpl::append<simple_tags, simple_tags_from_options>;
96 :
97 : template <typename... AmrData>
98 : static void apply(
99 : const gsl::not_null<std::optional<ElementId<Dim>>*> parent_id,
100 : const gsl::not_null<std::unordered_set<ElementId<Dim>>*> child_ids,
101 : const gsl::not_null<std::optional<Mesh<Dim>>*> parent_mesh,
102 : const gsl::not_null<size_t*> observation_id,
103 : const gsl::not_null<VolumeDataVars*> volume_data_for_output,
104 : const gsl::not_null<std::vector<std::array<size_t, Dim>>*>
105 : children_refinement_levels,
106 : const gsl::not_null<std::vector<std::array<size_t, Dim>>*>
107 : parent_refinement_levels,
108 : const Mesh<Dim>& mesh, const Element<Dim>& element,
109 : const std::vector<std::array<size_t, Dim>> initial_refinement_levels,
110 : const bool output_volume_data, const AmrData&... amr_data) {
111 : // Note: This initialization code runs on elements of the initial multigrid
112 : // hierarchy (created by LinearSolver::multigrid::ElementsAllocator) and on
113 : // elements created by AMR. Elements in a block of the initial domain are
114 : // assumed to have the same p-refinement.
115 :
116 : if constexpr (sizeof...(AmrData) == 0) {
117 : // Initialization: use initial domain to set up multigrid hierarchy
118 : const auto& element_id = element.id();
119 : const bool is_coarsest_grid =
120 : initial_refinement_levels == *parent_refinement_levels;
121 : const bool is_finest_grid =
122 : initial_refinement_levels == *children_refinement_levels;
123 : *parent_id = is_coarsest_grid
124 : ? std::nullopt
125 : : std::make_optional(multigrid::parent_id(element_id));
126 : *child_ids =
127 : is_finest_grid
128 : ? std::unordered_set<ElementId<Dim>>{}
129 : : multigrid::child_ids(
130 : element_id,
131 : (*children_refinement_levels)[element_id.block_id()]);
132 : *parent_mesh = is_coarsest_grid ? std::nullopt : std::make_optional(mesh);
133 : *observation_id = 0;
134 : } else {
135 : // These items are updated by AMR
136 : (void)parent_id;
137 : (void)child_ids;
138 : (void)parent_mesh;
139 : // These items are only needed during initialization
140 : (void)parent_refinement_levels;
141 : (void)children_refinement_levels;
142 : // Preserve state of observation ID
143 : if constexpr (tt::is_a_v<tuples::TaggedTuple, AmrData...>) {
144 : // h-refinement: copy from the parent
145 : *observation_id =
146 : get<LinearSolver::Tags::ObservationId<OptionsGroup>>(amr_data...);
147 : } else if constexpr (tt::is_a_v<std::unordered_map, AmrData...>) {
148 : // h-coarsening: copy from one of the children (doesn't matter which)
149 : *observation_id = get<LinearSolver::Tags::ObservationId<OptionsGroup>>(
150 : amr_data.begin()->second...);
151 : } else {
152 : (void)observation_id;
153 : }
154 : }
155 : // Initialize volume data output
156 : if (output_volume_data) {
157 : volume_data_for_output->initialize(mesh.number_of_grid_points());
158 : }
159 : }
160 : };
161 :
162 : // These two actions communicate and project the residual from the finer grid to
163 : // the coarser grid, storing it in the `SourceTag` on the coarser grid.
164 : template <typename FieldsTag, typename OptionsGroup,
165 : typename ResidualIsMassiveTag, typename SourceTag>
166 : using SendResidualToCoarserGrid = Actions::SendFieldsToCoarserGrid<
167 : tmpl::list<db::add_tag_prefix<LinearSolver::Tags::Residual, FieldsTag>>,
168 : OptionsGroup, ResidualIsMassiveTag, tmpl::list<SourceTag>>;
169 :
170 : template <size_t Dim, typename FieldsTag, typename OptionsGroup,
171 : typename SourceTag>
172 : using ReceiveResidualFromFinerGrid = Actions::ReceiveFieldsFromFinerGrid<
173 : Dim,
174 : tmpl::list<db::add_tag_prefix<LinearSolver::Tags::Residual, FieldsTag>>,
175 : OptionsGroup, tmpl::list<SourceTag>>;
176 :
177 : // Once the residual from the finer grid has been received and stored in the
178 : // `SourceTag`, this action prepares the pre-smoothing that will determine
179 : // an approximate solution on this grid. The pre-smoother is a separate
180 : // linear solver that runs independently after this action.
181 : template <typename FieldsTag, typename OptionsGroup, typename SourceTag,
182 : bool EnableBottomSolver>
183 : struct PreparePreSmoothing {
184 : private:
185 : using fields_tag = FieldsTag;
186 : using operator_applied_to_fields_tag =
187 : db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, fields_tag>;
188 : using source_tag = SourceTag;
189 :
190 : public:
191 : using const_global_cache_tags = tmpl::append<
192 : tmpl::list<
193 : LinearSolver::multigrid::Tags::EnablePreSmoothing<OptionsGroup>>,
194 : tmpl::conditional_t<EnableBottomSolver,
195 : tmpl::list<LinearSolver::multigrid::Tags::
196 : UseBottomSolver<OptionsGroup>>,
197 : tmpl::list<>>>;
198 :
199 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
200 : size_t Dim, typename ActionList, typename ParallelComponent>
201 : static Parallel::iterable_action_return_t apply(
202 : db::DataBox<DbTagsList>& box,
203 : tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
204 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
205 : const ElementId<Dim>& element_id, const ActionList /*meta*/,
206 : const ParallelComponent* const /*meta*/) {
207 : const size_t iteration_id =
208 : db::get<Convergence::Tags::IterationId<OptionsGroup>>(box);
209 : const bool is_coarsest_grid =
210 : not db::get<::amr::Tags::ParentId<Dim>>(box).has_value();
211 : if (UNLIKELY(db::get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
212 : ::Verbosity::Debug)) {
213 : Parallel::printf("%s %s(%zu): Prepare %s\n", element_id,
214 : pretty_type::name<OptionsGroup>(), iteration_id,
215 : (EnableBottomSolver and is_coarsest_grid)
216 : ? "bottom-solver"
217 : : "pre-smoothing");
218 : }
219 :
220 : // On coarser grids the smoother solves for a correction to the finer-grid
221 : // fields, so we set its initial guess to zero. On the finest grid we smooth
222 : // the fields directly, so there's nothing to prepare.
223 : const bool is_finest_grid = db::get<amr::Tags::ChildIds<Dim>>(box).empty();
224 : if (not is_finest_grid) {
225 : db::mutate<fields_tag, operator_applied_to_fields_tag>(
226 : [](const auto fields, const auto operator_applied_to_fields,
227 : const auto& source) {
228 : *fields = make_with_value<typename fields_tag::type>(source, 0.);
229 : // We can set the linear operator applied to the initial fields to
230 : // zero as well, since it's linear
231 : *operator_applied_to_fields =
232 : make_with_value<typename operator_applied_to_fields_tag::type>(
233 : source, 0.);
234 : },
235 : make_not_null(&box), db::get<source_tag>(box));
236 : }
237 :
238 : // Record pre-smoothing initial fields and source
239 : if (db::get<LinearSolver::Tags::OutputVolumeData<OptionsGroup>>(box)) {
240 : db::mutate<Tags::VolumeDataForOutput<OptionsGroup, FieldsTag>>(
241 : [](const auto volume_data, const auto& initial_fields,
242 : const auto& source) {
243 : volume_data->assign_subset(
244 : Variables<db::wrap_tags_in<Tags::PreSmoothingInitial,
245 : typename fields_tag::tags_list>>(
246 : initial_fields));
247 : volume_data->assign_subset(
248 : Variables<db::wrap_tags_in<Tags::PreSmoothingSource,
249 : typename fields_tag::tags_list>>(
250 : source));
251 : },
252 : make_not_null(&box), db::get<fields_tag>(box),
253 : db::get<source_tag>(box));
254 : }
255 :
256 : // Skip pre-smoothing, if requested, or use bottom solver
257 : if constexpr (EnableBottomSolver) {
258 : if (db::get<LinearSolver::multigrid::Tags::UseBottomSolver<OptionsGroup>>(
259 : box) and
260 : is_coarsest_grid) {
261 : const size_t bottom_solver_index =
262 : tmpl::index_of<ActionList, SkipBottomSolver<FieldsTag, OptionsGroup,
263 : SourceTag>>::value +
264 : 1;
265 : return {Parallel::AlgorithmExecution::Continue, bottom_solver_index};
266 : }
267 : }
268 : const size_t first_action_after_pre_smoothing_index = tmpl::index_of<
269 : ActionList,
270 : SkipPostSmoothingAtBottom<FieldsTag, OptionsGroup, SourceTag>>::value;
271 : const size_t this_action_index =
272 : tmpl::index_of<ActionList, PreparePreSmoothing>::value;
273 : return {
274 : Parallel::AlgorithmExecution::Continue,
275 : db::get<
276 : LinearSolver::multigrid::Tags::EnablePreSmoothing<OptionsGroup>>(
277 : box)
278 : ? (this_action_index + 1)
279 : : first_action_after_pre_smoothing_index};
280 : }
281 : };
282 :
283 : // Once pre-smoothing is done, skip the bottom solver that comes next in the
284 : // action list. On the coarsest grid we directly jump to the bottom solver.
285 : template <typename FieldsTag, typename OptionsGroup, typename SourceTag>
286 : struct SkipBottomSolver {
287 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
288 : size_t Dim, typename ActionList, typename ParallelComponent>
289 : static Parallel::iterable_action_return_t apply(
290 : db::DataBox<DbTagsList>& /*box*/,
291 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
292 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
293 : const ElementId<Dim>& /*element_id*/, const ActionList /*meta*/,
294 : const ParallelComponent* const /*meta*/) {
295 : const size_t first_action_after_bottom_solver_index = tmpl::index_of<
296 : ActionList,
297 : SkipPostSmoothingAtBottom<FieldsTag, OptionsGroup, SourceTag>>::value;
298 : return {Parallel::AlgorithmExecution::Continue,
299 : first_action_after_bottom_solver_index};
300 : }
301 : };
302 :
303 : // Once the pre-smoothing is done, we skip the second smoothing step on the
304 : // coarsest grid, i.e. at the "tip" of the V-cycle.
305 : template <typename FieldsTag, typename OptionsGroup, typename SourceTag>
306 : struct SkipPostSmoothingAtBottom {
307 : private:
308 : using fields_tag = FieldsTag;
309 : using residual_tag =
310 : db::add_tag_prefix<LinearSolver::Tags::Residual, fields_tag>;
311 :
312 : public:
313 : using const_global_cache_tags = tmpl::list<
314 : LinearSolver::multigrid::Tags::EnablePostSmoothingAtBottom<OptionsGroup>>;
315 :
316 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
317 : size_t Dim, typename ActionList, typename ParallelComponent>
318 : static Parallel::iterable_action_return_t apply(
319 : db::DataBox<DbTagsList>& box,
320 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
321 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
322 : const ElementId<Dim>& /*element_id*/, const ActionList /*meta*/,
323 : const ParallelComponent* const /*meta*/) {
324 : const bool is_coarsest_grid =
325 : not db::get<amr::Tags::ParentId<Dim>>(box).has_value();
326 :
327 : // Record pre-smoothing result fields and residual
328 : if (db::get<LinearSolver::Tags::OutputVolumeData<OptionsGroup>>(box)) {
329 : db::mutate<Tags::VolumeDataForOutput<OptionsGroup, FieldsTag>>(
330 : [](const auto volume_data, const auto& result_fields,
331 : const auto& residuals) {
332 : volume_data->assign_subset(
333 : Variables<db::wrap_tags_in<Tags::PreSmoothingResult,
334 : typename fields_tag::tags_list>>(
335 : result_fields));
336 : volume_data->assign_subset(
337 : Variables<db::wrap_tags_in<Tags::PreSmoothingResidual,
338 : typename fields_tag::tags_list>>(
339 : residuals));
340 : },
341 : make_not_null(&box), db::get<fields_tag>(box),
342 : db::get<residual_tag>(box));
343 : }
344 :
345 : // Skip post-smoothing on the coarsest grid, if requested
346 : const size_t first_action_after_post_smoothing_index = tmpl::index_of<
347 : ActionList,
348 : SendCorrectionToFinerGrid<FieldsTag, OptionsGroup, SourceTag>>::value;
349 : const size_t post_smoothing_begin_index =
350 : tmpl::index_of<ActionList,
351 : ::Actions::Label<PostSmoothingBeginLabel>>::value +
352 : 1;
353 : const size_t this_action_index =
354 : tmpl::index_of<ActionList, SkipPostSmoothingAtBottom>::value;
355 : return {Parallel::AlgorithmExecution::Continue,
356 : is_coarsest_grid
357 : ? (db::get<LinearSolver::multigrid::Tags::
358 : EnablePostSmoothingAtBottom<OptionsGroup>>(box)
359 : ? post_smoothing_begin_index
360 : : first_action_after_post_smoothing_index)
361 : : (this_action_index + 1)};
362 : }
363 : };
364 :
365 : template <typename FieldsTag>
366 : struct CorrectionInboxTag
367 : : public Parallel::InboxInserters::Value<CorrectionInboxTag<FieldsTag>> {
368 : using temporal_id = size_t;
369 : using type = std::map<temporal_id, typename FieldsTag::type>;
370 : };
371 :
372 : // The next two actions communicate and project the coarse-grid correction, i.e.
373 : // the solution of the post-smoother, to the finer grid. The post-smoother on
374 : // finer grids runs after receiving this coarse-grid correction. Since the
375 : // post-smoother is skipped on the coarsest level, it directly sends the
376 : // solution of the pre-smoother to the finer grid, thus kicking off the
377 : // "ascending" branch of the V-cycle.
378 : template <typename FieldsTag, typename OptionsGroup, typename SourceTag>
379 : struct SendCorrectionToFinerGrid {
380 : private:
381 : using fields_tag = FieldsTag;
382 : using residual_tag =
383 : db::add_tag_prefix<LinearSolver::Tags::Residual, fields_tag>;
384 :
385 : public:
386 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
387 : size_t Dim, typename ActionList, typename ParallelComponent>
388 : static Parallel::iterable_action_return_t apply(
389 : db::DataBox<DbTagsList>& box,
390 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
391 : Parallel::GlobalCache<Metavariables>& cache,
392 : const ElementId<Dim>& element_id, const ActionList /*meta*/,
393 : const ParallelComponent* const /*meta*/) {
394 : const auto& child_ids = db::get<amr::Tags::ChildIds<Dim>>(box);
395 :
396 : // Record post-smoothing result fields and residual
397 : if (db::get<LinearSolver::Tags::OutputVolumeData<OptionsGroup>>(box)) {
398 : db::mutate<Tags::VolumeDataForOutput<OptionsGroup, FieldsTag>>(
399 : [](const auto volume_data, const auto& result_fields,
400 : const auto& residuals) {
401 : volume_data->assign_subset(
402 : Variables<db::wrap_tags_in<Tags::PostSmoothingResult,
403 : typename fields_tag::tags_list>>(
404 : result_fields));
405 : volume_data->assign_subset(
406 : Variables<db::wrap_tags_in<Tags::PostSmoothingResidual,
407 : typename fields_tag::tags_list>>(
408 : residuals));
409 : },
410 : make_not_null(&box), db::get<fields_tag>(box),
411 : db::get<residual_tag>(box));
412 : }
413 :
414 : if (child_ids.empty()) {
415 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
416 : }
417 :
418 : const size_t iteration_id =
419 : db::get<Convergence::Tags::IterationId<OptionsGroup>>(box);
420 : if (UNLIKELY(db::get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
421 : ::Verbosity::Debug)) {
422 : Parallel::printf("%s %s(%zu): Send correction to children\n", element_id,
423 : pretty_type::name<OptionsGroup>(), iteration_id);
424 : }
425 :
426 : // Send a copy of the correction to all children
427 : auto& receiver_proxy =
428 : Parallel::get_parallel_component<ParallelComponent>(cache);
429 : for (const auto& child_id : child_ids) {
430 : auto coarse_grid_correction = db::get<fields_tag>(box);
431 : Parallel::receive_data<CorrectionInboxTag<FieldsTag>>(
432 : receiver_proxy[child_id], iteration_id,
433 : std::move(coarse_grid_correction));
434 : }
435 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
436 : }
437 : };
438 :
439 : template <size_t Dim, typename FieldsTag, typename OptionsGroup,
440 : typename SourceTag>
441 : struct ReceiveCorrectionFromCoarserGrid {
442 : private:
443 : using fields_tag = FieldsTag;
444 : using source_tag = SourceTag;
445 :
446 : public:
447 : using inbox_tags = tmpl::list<CorrectionInboxTag<FieldsTag>>;
448 :
449 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
450 : typename ActionList, typename ParallelComponent>
451 : static Parallel::iterable_action_return_t apply(
452 : db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
453 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
454 : const ElementId<Dim>& element_id, const ActionList /*meta*/,
455 : const ParallelComponent* const /*meta*/) {
456 : const auto& parent_id = db::get<amr::Tags::ParentId<Dim>>(box);
457 : // We should always have a `parent_id` at this point because we skip this
458 : // part of the algorithm on the coarsest grid with the
459 : // `SkipPostSmoothingAtBottom` action
460 : ASSERT(parent_id.has_value(),
461 : "Trying to receive data from parent but no parent is set on element "
462 : << element_id << ".");
463 : const size_t iteration_id =
464 : db::get<Convergence::Tags::IterationId<OptionsGroup>>(box);
465 :
466 : // Wait for data from coarser grid
467 : auto& inbox = tuples::get<CorrectionInboxTag<FieldsTag>>(inboxes);
468 : if (inbox.find(iteration_id) == inbox.end()) {
469 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
470 : }
471 : auto parent_correction = std::move(inbox.extract(iteration_id).mapped());
472 :
473 : if (UNLIKELY(db::get<logging::Tags::Verbosity<OptionsGroup>>(box) >=
474 : ::Verbosity::Debug)) {
475 : Parallel::printf("%s %s(%zu): Prolongate correction from parent\n",
476 : element_id, pretty_type::name<OptionsGroup>(),
477 : iteration_id);
478 : }
479 :
480 : // Apply prolongation operator
481 : const auto& mesh = db::get<domain::Tags::Mesh<Dim>>(box);
482 : const auto& parent_mesh = db::get<amr::Tags::ParentMesh<Dim>>(box);
483 : ASSERT(
484 : parent_mesh.has_value(),
485 : "Should have a parent mesh, because a parent ID is set. This element: "
486 : << element_id << ", parent element: " << *parent_id);
487 : const auto child_size =
488 : domain::child_size(element_id.segment_ids(), parent_id->segment_ids());
489 : const auto prolongated_parent_correction =
490 : [&parent_correction, &parent_mesh, &mesh, &child_size]() {
491 : if (Spectral::needs_projection(*parent_mesh, mesh, child_size)) {
492 : const auto prolongation_operator =
493 : Spectral::projection_matrix_parent_to_child(*parent_mesh, mesh,
494 : child_size);
495 : return apply_matrices(prolongation_operator, parent_correction,
496 : parent_mesh->extents());
497 : } else {
498 : return std::move(parent_correction);
499 : }
500 : }();
501 :
502 : // Add correction to the solution on this grid
503 : db::mutate<fields_tag>(
504 : [&prolongated_parent_correction](const auto fields) {
505 : *fields += prolongated_parent_correction;
506 : },
507 : make_not_null(&box));
508 :
509 : // Record post-smoothing initial fields and source
510 : if (db::get<LinearSolver::Tags::OutputVolumeData<OptionsGroup>>(box)) {
511 : db::mutate<Tags::VolumeDataForOutput<OptionsGroup, FieldsTag>>(
512 : [](const auto volume_data, const auto& initial_fields,
513 : const auto& source) {
514 : volume_data->assign_subset(
515 : Variables<db::wrap_tags_in<Tags::PostSmoothingInitial,
516 : typename fields_tag::tags_list>>(
517 : initial_fields));
518 : volume_data->assign_subset(
519 : Variables<db::wrap_tags_in<Tags::PostSmoothingSource,
520 : typename fields_tag::tags_list>>(
521 : source));
522 : },
523 : make_not_null(&box), db::get<fields_tag>(box),
524 : db::get<source_tag>(box));
525 : }
526 :
527 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
528 : }
529 : };
530 :
531 : } // namespace LinearSolver::multigrid::detail
|
