Schwarz.hpp
1 // 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/Protocols.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::protocols::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, we use
120  * `LinearSolver::Serial::Gmres` to solve subdomains. It supports
121  * preconditioning, and adding useful subdomain preconditioners will be the
122  * subject of future work. Note that the choice of subdomain solver (and, by
123  * extension, the choice of subdomain preconditioner) affects only the
124  * _performance_ of the Schwarz solver, not its convergence or parallelization
125  * properties (assuming the subdomain solutions it produces are sufficiently
126  * precise).
127  *
128  * \par Weighting:
129  * Once the subdomain solutions \f$\delta x_s\f$ have been found they must be
130  * combined where they have multiple values, i.e. on overlap regions of the
131  * subdomains. We use an additive approach where we combine the subdomain
132  * solutions as a weighted sum, which has the advantage over "multiplicative"
133  * Schwarz methods that the subdomains decouple and can be solved in parallel.
134  * See `LinearSolver::Schwarz::element_weight` and Section 3.1 of
135  * \cite Stiller2016a for details. Note that we must account for the missing
136  * corner- and edge-neighbors when constructing the weights. See
137  * `LinearSolver::Schwarz::intruding_weight` for a discussion.
138  *
139  * \par Possible improvements:
140  * - Specify the number of overlap points as a fraction of the element width
141  * instead of a fixed number. This was shown in \cite Stiller2016b to achieve
142  * scale-independent convergence at the cost of increasing subdomain sizes.
143  * - Subdomain preconditioning (see paragraph on subdomain solves above)
144  */
145 template <typename FieldsTag, typename OptionsGroup, typename SubdomainOperator,
146  typename SourceTag =
148 struct Schwarz {
149  using operand_tag = FieldsTag;
150  using fields_tag = FieldsTag;
151  using source_tag = SourceTag;
152  using options_group = OptionsGroup;
153  static_assert(
154  tt::assert_conforms_to<SubdomainOperator, protocols::SubdomainOperator>);
155 
156  using component_list = tmpl::list<>;
157  using observed_reduction_data_tags = observers::make_reduction_data_tags<
158  tmpl::list<async_solvers::reduction_data, detail::reduction_data>>;
159  using subdomain_solver =
161 
162  using initialize_element = tmpl::list<
164  detail::InitializeElement<FieldsTag, OptionsGroup, SubdomainOperator>>;
165 
166  using register_element = tmpl::list<
169 
170  template <typename ApplyOperatorActions, typename Label = OptionsGroup>
171  using solve = tmpl::list<
175  detail::SolveSubdomain<FieldsTag, OptionsGroup, SubdomainOperator>,
176  detail::ReceiveOverlapSolution<FieldsTag, OptionsGroup,
177  SubdomainOperator>,
178  ApplyOperatorActions,
180 };
181 
182 } // namespace LinearSolver::Schwarz
LinearSolver::Serial::LinearSolver
Base class for serial linear solvers that supports factory-creation.
Definition: LinearSolver.hpp:37
LinearSolver::Schwarz::Actions::ReceiveOverlapFields
Receive data from regions of this element's subdomain that overlap with other elements.
Definition: CommunicateOverlapFields.hpp:160
db::add_tag_prefix
typename detail::add_tag_prefix_impl< Prefix, Tag, Args... >::type add_tag_prefix
Definition: PrefixHelpers.hpp:51
LinearSolver::async_solvers::InitializeElement
Definition: ElementActions.hpp:102
LinearSolver::Schwarz
Items related to the Schwarz linear solver.
Definition: CommunicateOverlapFields.hpp:36
LinearSolver::Schwarz::Actions::SendOverlapFields
Send data on regions that overlap with other subdomains to their corresponding elements.
Definition: CommunicateOverlapFields.hpp:68
LinearSolver::Schwarz::Schwarz
An additive Schwarz subdomain solver for linear systems of equations .
Definition: Schwarz.hpp:148
LinearSolver::async_solvers::PrepareSolve
Definition: ElementActions.hpp:160
TMPL.hpp
observers::Actions::RegisterWithObservers
Register an observation ID with the observers.
Definition: RegisterWithObservers.hpp:39
LinearSolver::async_solvers::CompleteStep
Definition: ElementActions.hpp:214