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