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 :
8 : #include "DataStructures/DataBox/Prefixes.hpp"
9 : #include "IO/Observer/Actions/RegisterWithObservers.hpp"
10 : #include "IO/Observer/Helpers.hpp"
11 : #include "ParallelAlgorithms/Actions/Goto.hpp"
12 : #include "ParallelAlgorithms/LinearSolver/AsynchronousSolvers/ElementActions.hpp"
13 : #include "ParallelAlgorithms/LinearSolver/Multigrid/ElementActions.hpp"
14 : #include "ParallelAlgorithms/LinearSolver/Multigrid/ObserveVolumeData.hpp"
15 : #include "Utilities/TMPL.hpp"
16 :
17 : /// Items related to the multigrid linear solver
18 : ///
19 : /// \see `LinearSolver::multigrid::Multigrid`
20 : namespace LinearSolver::multigrid {
21 :
22 : /// A label indicating the pre-smoothing step in a V-cycle multigrid algorithm,
23 : /// i.e. the smoothing step before sending the residual to the coarser (parent)
24 : /// grid
25 1 : struct VcycleDownLabel {};
26 :
27 : /// A label indicating the post-smoothing step in a V-cycle multigrid algorithm,
28 : /// i.e. the smoothing step before sending the correction to the finer (child)
29 : /// grid
30 1 : struct VcycleUpLabel {};
31 :
32 : /*!
33 : * \brief A V-cycle geometric multgrid solver for linear equations \f$Ax = b\f$
34 : *
35 : * This linear solver iteratively corrects an initial guess \f$x_0\f$ by
36 : * restricting the residual \f$b - Ax\f$ to a series of coarser grids, solving
37 : * for a correction on the coarser grids, and then prolongating (interpolating)
38 : * the correction back to the finer grids. The solves on grids with different
39 : * scales can very effectively solve large-scale modes in the solution, which
40 : * Krylov-type linear solvers such as GMRES or Conjugate Gradients typically
41 : * struggle with. Therefore, a multigrid solver can be an effective
42 : * preconditioner for Krylov-type linear solvers (see
43 : * `LinearSolver::gmres::Gmres` and `LinearSolver::cg::ConjugateGradient`). See
44 : * \cite Briggs2000jp for an introduction to multigrid methods.
45 : *
46 : * \par Grid hierarchy
47 : * This geometric multigrid solver relies on a strategy to coarsen the
48 : * computational grid in a way that removes small-scale modes. We currently
49 : * h-coarsen the initial domain, meaning that we create multigrid levels by
50 : * successively combining two elements into one along every dimension of the
51 : * grid. We only p-coarsen the grid in the sense that we choose the smaller of
52 : * the two polynomial degrees when combining elements. This strategy follows
53 : * \cite Vincent2019qpd. See `LinearSolver::multigrid::ElementsAllocator` and
54 : * `LinearSolver::multigrid::coarsen` for the code that creates the initial
55 : * multigrid hierarchy.
56 : * Once the initial grids are created, AMR steps can create incrementally finer
57 : * grids. For this to work, AMR must be configured with the compile-time option
58 : * `keep_coarse_grids = true` (see `amr::protocols::AmrMetavariables`).
59 : *
60 : * \par Inter-mesh operators
61 : * The algorithm relies on operations that project data between grids. Residuals
62 : * are projected from finer to coarser grids ("restriction") and solutions are
63 : * projected from coarser to finer grids ("prolongation"). We use the standard
64 : * \f$L_2\f$-projections implemented in
65 : * `Spectral::projection_matrix_child_to_parent` and
66 : * `Spectral::projection_matrix_parent_to_child`, where "child" means the finer
67 : * grid and "parent" means the coarser grid. Note that the residuals \f$b -
68 : * Ax\f$ may or may not already include a mass matrix and Jacobian factors,
69 : * depending on the implementation of the linear operator \f$A\f$. To account
70 : * for this, provide a tag as the `ResidualIsMassiveTag` that holds a `bool`.
71 : * See section 3 in \cite Fortunato2019jl for details on the inter-mesh
72 : * operators.
73 : *
74 : * \par Smoother and bottom solver
75 : * On every level of the grid hierarchy the multigrid solver relies on a
76 : * "smoother" that performs an approximate linear solve on that level. The
77 : * smoother can be any linear solver that solves the `smooth_fields_tag` for the
78 : * source in the `smooth_source_tag`. Useful smoothers are asynchronous linear
79 : * solvers that parallelize well, such as `LinearSolver::Schwarz::Schwarz`. The
80 : * multigrid algorithm doesn't assume anything about the smoother, but requires
81 : * only that the `PreSmootherActions` and the `PostSmootherActions` passed to
82 : * `solve` leave the `smooth_fields_tag` in a state that represents an
83 : * approximate solution to the `smooth_source_tag`, and that
84 : * `db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo,
85 : * smooth_fields_tag>` is left up-to-date as well.
86 : * The smoother can assume that, on entry, these two tags represent the initial
87 : * fields and the linear operator applied to the initial fields, respectively.
88 : * Here's an example of setting up a smoother for the multigrid solver:
89 : *
90 : * \snippet Test_MultigridAlgorithm.cpp setup_smoother
91 : *
92 : * On the coarsest level of the grid hierarchy, the multigrid algorithm can also
93 : * use a separate "bottom solver". If enabled, the bottom solver replaces the
94 : * pre-smoother on the coarsest grid. The bottom solver can also be any linear
95 : * solver, but it is typically a direct solver that solves the linear problem
96 : * exactly (see `LinearSolver::Actions::BuildMatrix` for a good bottom solver).
97 : *
98 : * The smoother and bottom solver can be used to construct an action list like
99 : * this:
100 : *
101 : * \snippet Test_MultigridAlgorithm.cpp action_list
102 : *
103 : * \par Algorithm overview
104 : * Every iteration of the multigrid algorithm performs a V-cycle over the grid
105 : * hierarchy. One V-cycle consists of first "going down" the grid hierarchy,
106 : * from the finest to successively coarser grids, smoothing on every level
107 : * ("pre-smoothing", can be controlled by the option
108 : * `LinearSolver::multigrid::Tags::EnablePreSmoothing`),
109 : * and then "going up" the grid hierarchy again, smoothing on
110 : * every level again ("post-smoothing"). When going down, the algorithm projects
111 : * the remaining residual of the smoother to the next-coarser grid, setting it
112 : * as the source for the smoother on the coarser grid. When going up again, the
113 : * algorithm projects the solution of the smoother to the next-finer grid,
114 : * adding it to the solution on the finer grid as a correction. The bottom-most
115 : * coarsest grid (the "tip" of the V-cycle) may run the bottom solver instead of
116 : * the pre-smoother. It may also skip the post-smoother, so the result of the
117 : * bottom solver is immediately projected up to the finer grid (controlled by
118 : * the options `LinearSolver::multigrid::Tags::UseBottomSolver` and
119 : * `LinearSolver::multigrid::Tags::EnablePostSmoothingAtBottom`). On the
120 : * top-most finest grid (the "original" grid that represents the overall
121 : * solution) the algorithm applies the smoothing and the corrections from the
122 : * coarser grids directly to the solution fields.
123 : */
124 : template <typename Metavariables, size_t Dim, typename FieldsTag,
125 : typename OptionsGroup, typename ResidualIsMassiveTag,
126 : typename SourceTag =
127 : db::add_tag_prefix<::Tags::FixedSource, FieldsTag>>
128 1 : struct Multigrid {
129 0 : using fields_tag = FieldsTag;
130 0 : using options_group = OptionsGroup;
131 0 : using source_tag = SourceTag;
132 :
133 0 : using operand_tag = FieldsTag;
134 :
135 0 : using smooth_source_tag = source_tag;
136 0 : using smooth_fields_tag = fields_tag;
137 :
138 0 : using component_list = tmpl::list<>;
139 :
140 0 : using observed_reduction_data_tags = observers::make_reduction_data_tags<
141 : tmpl::list<async_solvers::reduction_data>>;
142 :
143 0 : using initialize_element = tmpl::list<
144 : async_solvers::InitializeElement<FieldsTag, OptionsGroup, SourceTag>,
145 : detail::InitializeElement<Dim, FieldsTag, OptionsGroup, SourceTag>>;
146 :
147 0 : using amr_projectors = initialize_element;
148 :
149 0 : using register_element = tmpl::list<
150 : async_solvers::RegisterElement<FieldsTag, OptionsGroup, SourceTag,
151 : ::amr::Tags::IsFinestGrid>,
152 : observers::Actions::RegisterWithObservers<
153 : detail::RegisterWithVolumeObserver<OptionsGroup>>>;
154 :
155 : template <typename ApplyOperatorActions, typename PreSmootherActions,
156 : typename PostSmootherActions,
157 : typename BottomSolverActions = tmpl::list<>,
158 : typename Label = OptionsGroup>
159 0 : using solve = tmpl::list<
160 : async_solvers::PrepareSolve<FieldsTag, OptionsGroup, SourceTag, Label,
161 : ::amr::Tags::IsFinestGrid, false>,
162 : detail::ReceiveResidualFromFinerGrid<Dim, FieldsTag, OptionsGroup,
163 : SourceTag>,
164 : detail::PreparePreSmoothing<
165 : FieldsTag, OptionsGroup, SourceTag,
166 : not std::is_same_v<BottomSolverActions, tmpl::list<>>>,
167 : // No need to apply the linear operator here:
168 : // - On the finest grid, the operator applied to the fields should have
169 : // already been computed at this point, either applied to the initial
170 : // fields before the multigrid solver is invoked, or by the smoother at
171 : // the end of the previous V-cycle.
172 : // - On coarser grids, the initial fields are zero, so the operator
173 : // applied to them is also zero.
174 : PreSmootherActions,
175 : detail::SkipBottomSolver<FieldsTag, OptionsGroup, SourceTag>,
176 : BottomSolverActions,
177 : detail::SkipPostSmoothingAtBottom<FieldsTag, OptionsGroup, SourceTag>,
178 : detail::SendResidualToCoarserGrid<FieldsTag, OptionsGroup,
179 : ResidualIsMassiveTag, SourceTag>,
180 : detail::ReceiveCorrectionFromCoarserGrid<Dim, FieldsTag, OptionsGroup,
181 : SourceTag>,
182 : ApplyOperatorActions, ::Actions::Label<detail::PostSmoothingBeginLabel>,
183 : PostSmootherActions,
184 : detail::SendCorrectionToFinerGrid<FieldsTag, OptionsGroup, SourceTag>,
185 : detail::ObserveVolumeData<FieldsTag, OptionsGroup, SourceTag>,
186 : async_solvers::CompleteStep<FieldsTag, OptionsGroup, SourceTag, Label,
187 : ::amr::Tags::IsFinestGrid, false>>;
188 : };
189 :
190 : } // namespace LinearSolver::multigrid
|
