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