Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include "IO/Observer/Helpers.hpp"
7 : #include "ParallelAlgorithms/LinearSolver/AsynchronousSolvers/ElementActions.hpp"
8 : #include "ParallelAlgorithms/LinearSolver/Schwarz/ElementActions.hpp"
9 : #include "ParallelAlgorithms/LinearSolver/Schwarz/ObserveVolumeData.hpp"
10 : #include "Utilities/ProtocolHelpers.hpp"
11 : #include "Utilities/TMPL.hpp"
12 :
13 : /// Items related to the %Schwarz linear solver
14 : ///
15 : /// \see `LinearSolver::Schwarz::Schwarz`
16 : namespace LinearSolver::Schwarz {
17 :
18 : /*!
19 : * \ingroup LinearSolverGroup
20 : * \brief An additive Schwarz subdomain solver for linear systems of equations
21 : * \f$Ax=b\f$.
22 : *
23 : * This Schwarz-type linear solver works by solving many sub-problems in
24 : * parallel and combining their solutions as a weighted sum to converge towards
25 : * the global solution. Each sub-problem is the restriction of the global
26 : * problem to an element-centered subdomain, which consists of a central element
27 : * and an overlap region with its neighbors. The decomposition into independent
28 : * sub-problems makes this linear solver very parallelizable, avoiding global
29 : * synchronization points altogether. It is commonly used as a preconditioner to
30 : * Krylov-type solvers such as `LinearSolver::gmres::Gmres` or
31 : * `LinearSolver::cg::ConjugateGradient`, or as the "smoother" for a Multigrid
32 : * solver (which may in turn precondition a Krylov-type solver).
33 : *
34 : * This linear solver relies on an implementation of the global linear operator
35 : * \f$A(x)\f$ and its restriction to a subdomain \f$A_{ss}(x)\f$. Each step of
36 : * the algorithm expects that \f$A(x)\f$ is computed and stored in the DataBox
37 : * as `db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, operand_tag>`.
38 : * To perform a solve, add the `solve` action list to an array parallel
39 : * component. Pass the actions that compute \f$A(x)\f$, as well as any further
40 : * actions you wish to run in each step of the algorithm, as the first template
41 : * parameter to `solve`. If you add the `solve` action list multiple times, use
42 : * the second template parameter to label each solve with a different type.
43 : * Pass the subdomain operator as the `SubdomainOperator` template parameter to
44 : * this class (not to the `solve` action list template). See the paragraph below
45 : * for information on implementing a subdomain operator.
46 : *
47 : * \par Subdomain geometry:
48 : * The image below gives an overview of the structure of an element-centered
49 : * subdomain:
50 : *
51 : * \image html subdomain_structure.svg "Fig. 1: An element-centered subdomain" width=600px
52 : *
53 : * Fig. 1 shows part of a 2-dimensional computational domain. The domain is
54 : * composed of elements (light gray) that each carry a Legendre-Gauss-Lobatto
55 : * mesh of grid points (black dots). The diagonally shaded region to the left
56 : * illustrates an external domain boundary. The lines between neighboring
57 : * element faces illustrate mortar meshes, which are relevant when implementing
58 : * a subdomain operator in a DG context but play no role in the Schwarz
59 : * algorithm. The dashed line gives an example of an element-centered subdomain
60 : * with a maximum of 2 grid points of overlap into neighboring elements. For
61 : * details on how the number of overlap points is chosen see
62 : * `LinearSolver::Schwarz::overlap_extent`. Note that the subdomain does not
63 : * extend into corner- or edge-neighbors, which is different to both
64 : * \cite Stiller2016a and \cite Vincent2019qpd. The reason we don't include
65 : * such diagonal couplings is that in a DG context information only propagates
66 : * across faces, as noted in \cite Stiller2016a. By eliminating the corner-
67 : * and edge-neighbors we significantly reduce the complexity of the subdomain
68 : * geometry and potentially also the communication costs. However, the
69 : * advantages and disadvantages of this choice have yet to be carefully
70 : * evaluated.
71 : *
72 : * \par Subdomain operator:
73 : * The Schwarz subdomain solver relies on a restriction of the global linear
74 : * problem \f$Ax=b\f$ to the subdomains. The subdomain operator
75 : * \f$A_{ss}=R_s A R_s^T\f$, where \f$R_s\f$ is the restriction operator,
76 : * essentially assumes that its operand is zero everywhere but within the
77 : * subdomain (see Section 3.1 in \cite Stiller2016a). Therefore it can be
78 : * evaluated entirely on an element-centered subdomain with no need to
79 : * communicate with neighbors within each subdomain-operator application, as
80 : * opposed to the global linear operator that typically requires
81 : * nearest-neighbor communications. See
82 : * `LinearSolver::Schwarz::SubdomainOperator` for details on how to
83 : * implement a subdomain operator for your problem.
84 : *
85 : * \par Algorithm overview:
86 : * In each iteration, the Schwarz algorithm computes the residual \f$r=b-Ax\f$,
87 : * restricts it to all subdomains as \f$r_s=R_s r\f$ and communicates data on
88 : * overlap regions with neighboring elements. Once an element has received all
89 : * data on its subdomain it solves the sub-problem \f$A_{ss}\delta x_s=r_s\f$
90 : * for the correction \f$\delta x_s\f$, where \f$\delta x_s\f$ and \f$r_s\f$
91 : * have data only on the points of the subdomain and \f$A_{ss}\f$ is the
92 : * subdomain operator. Since all elements perform such a subdomain-solve, we end
93 : * up with a subdomain solution \f$\delta x_s\f$ on every subdomain, and the
94 : * solutions overlap. Therefore, the algorithm communicates subdomain solutions
95 : * on overlap regions with neighboring elements and adds them to the solution
96 : * field \f$x\f$ as a weighted sum.
97 : *
98 : * \par Applying the global linear operator:
99 : * In order to compute the residual \f$r=b-Ax\f$ that will be restricted to the
100 : * subdomains to serve as source for the subdomain solves, we must apply the
101 : * global linear operator \f$A\f$ to the solution field \f$x\f$ once per Schwarz
102 : * iteration. The global linear operator typically introduces nearest-neighbor
103 : * couplings between elements but no global synchronization point (as opposed to
104 : * Krylov-type solvers such as `LinearSolver::gmres::Gmres` that require global
105 : * inner products in each iteration). The algorithm assumes that the action list
106 : * passed to `solve` applies the global linear operator, as described above.
107 : *
108 : * \par Subdomain solves:
109 : * Each subdomain solve is local to an element, since the data on the subdomain
110 : * has been made available through communication with the neighboring elements.
111 : * We can now use any means to solve the subdomain problem that's appropriate
112 : * for the subdomain operator \f$A_{ss}\f$. For example, if the subdomain
113 : * operator was available as an explicit matrix we could invert it directly.
114 : * Since it is typically a matrix-free operator a common approach to solve the
115 : * subdomain problem is by means of an iterative Krylov-type method, such as
116 : * GMRES or Conjugate Gradient, ideally with an appropriate preconditioner (yes,
117 : * this would be preconditioned Krylov-methods solving the subdomains of the
118 : * Schwarz solver, which might in turn precondition a global Krylov-solver -
119 : * it's preconditioners all the way down). Currently, the linear solvers listed
120 : * in `LinearSolver::Serial` are available to solve subdomains. They include
121 : * Krylov-type methods that support nested preconditioning with any other linear
122 : * solver. Additional problem-specific subdomain preconditioners can be added
123 : * with the `SubdomainPreconditioners` template parameter. Note that the choice
124 : * of subdomain solver (and, by extension, the choice of subdomain
125 : * preconditioner) affects only the _performance_ of the Schwarz solver, not its
126 : * convergence or parallelization properties (assuming the subdomain solutions
127 : * it produces are sufficiently precise).
128 : *
129 : * \par Weighting:
130 : * Once the subdomain solutions \f$\delta x_s\f$ have been found they must be
131 : * combined where they have multiple values, i.e. on overlap regions of the
132 : * subdomains. We use an additive approach where we combine the subdomain
133 : * solutions as a weighted sum, which has the advantage over "multiplicative"
134 : * Schwarz methods that the subdomains decouple and can be solved in parallel.
135 : * See `LinearSolver::Schwarz::element_weight` and Section 3.1 of
136 : * \cite Stiller2016a for details. Note that we must account for the missing
137 : * corner- and edge-neighbors when constructing the weights. See
138 : * `LinearSolver::Schwarz::intruding_weight` for a discussion.
139 : *
140 : * \par Array sections
141 : * This linear solver requires no synchronization between elements, so it runs
142 : * on all elements in the array parallel component. Partitioning of the elements
143 : * in sections is only relevant for observing residual norms. Pass the section
144 : * ID tag for the `ArraySectionIdTag` template parameter if residual norms
145 : * should be computed over a section. Pass `void` (default) to compute residual
146 : * norms over all elements in the array.
147 : *
148 : * \par Possible improvements:
149 : * - Specify the number of overlap points as a fraction of the element width
150 : * instead of a fixed number. This was shown in \cite Stiller2016b to achieve
151 : * scale-independent convergence at the cost of increasing subdomain sizes.
152 : */
153 : template <typename FieldsTag, typename OptionsGroup, typename SubdomainOperator,
154 : typename SubdomainPreconditioners = tmpl::list<>,
155 : typename SourceTag =
156 : db::add_tag_prefix<::Tags::FixedSource, FieldsTag>,
157 : typename ArraySectionIdTag = void>
158 1 : struct Schwarz {
159 0 : using operand_tag = FieldsTag;
160 0 : using fields_tag = FieldsTag;
161 0 : using source_tag = SourceTag;
162 0 : using options_group = OptionsGroup;
163 :
164 0 : using component_list = tmpl::list<>;
165 0 : using observed_reduction_data_tags = observers::make_reduction_data_tags<
166 : tmpl::list<async_solvers::reduction_data, detail::reduction_data>>;
167 0 : using subdomain_solver =
168 : detail::subdomain_solver<FieldsTag, SubdomainOperator,
169 : SubdomainPreconditioners>;
170 :
171 0 : using initialize_element = tmpl::list<
172 : async_solvers::InitializeElement<FieldsTag, OptionsGroup, SourceTag>,
173 : detail::InitializeElement<FieldsTag, OptionsGroup, SubdomainOperator,
174 : SubdomainPreconditioners>>;
175 :
176 0 : using amr_projectors = initialize_element;
177 :
178 0 : using register_element = tmpl::list<
179 : async_solvers::RegisterElement<FieldsTag, OptionsGroup, SourceTag,
180 : ArraySectionIdTag>,
181 : detail::RegisterElement<FieldsTag, OptionsGroup, SourceTag,
182 : ArraySectionIdTag>,
183 : observers::Actions::RegisterWithObservers<
184 : detail::RegisterWithVolumeObserver<OptionsGroup, ArraySectionIdTag>>>;
185 :
186 : template <typename ApplyOperatorActions, typename Label = OptionsGroup>
187 0 : using solve = tmpl::list<
188 : async_solvers::PrepareSolve<FieldsTag, OptionsGroup, SourceTag, Label,
189 : ArraySectionIdTag>,
190 : detail::SendOverlapData<FieldsTag, OptionsGroup, SubdomainOperator>,
191 : detail::SolveSubdomain<FieldsTag, OptionsGroup, SubdomainOperator,
192 : SubdomainPreconditioners, ArraySectionIdTag>,
193 : detail::ReceiveOverlapSolution<FieldsTag, OptionsGroup,
194 : SubdomainOperator>,
195 : detail::ObserveVolumeData<FieldsTag, OptionsGroup, SubdomainOperator,
196 : ArraySectionIdTag>,
197 : ApplyOperatorActions,
198 : async_solvers::CompleteStep<FieldsTag, OptionsGroup, SourceTag, Label,
199 : ArraySectionIdTag>>;
200 : };
201 :
202 : } // namespace LinearSolver::Schwarz
|
