SpECTRE
v2021.04.06
Documentation
Introduction
Releases
Installation
User Tutorials
Dev Guide
Code of Conduct
Contributing Guide
Code Reference
Topics
Namespaces
Files
Bibliography
View on GitHub
src
ParallelAlgorithms
LinearSolver
Schwarz
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 "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 Possible improvements:
139
* - Specify the number of overlap points as a fraction of the element width
140
* instead of a fixed number. This was shown in \cite Stiller2016b to achieve
141
* scale-independent convergence at the cost of increasing subdomain sizes.
142
* - Subdomain preconditioning (see paragraph on subdomain solves above)
143
*/
144
template
<
typename
FieldsTag,
typename
OptionsGroup,
typename
SubdomainOperator,
145
typename
SourceTag =
146
db::add_tag_prefix<::Tags::FixedSource, FieldsTag>
>
147
struct
Schwarz
{
148
using
operand_tag = FieldsTag;
149
using
fields_tag = FieldsTag;
150
using
source_tag = SourceTag;
151
using
options_group = OptionsGroup;
152
153
using
component_list = tmpl::list<>;
154
using
observed_reduction_data_tags = observers::make_reduction_data_tags<
155
tmpl::list<async_solvers::reduction_data, detail::reduction_data>>;
156
using
subdomain_solver
=
157
detail::subdomain_solver<FieldsTag, SubdomainOperator>
;
158
159
using
initialize_element = tmpl::list<
160
async_solvers::InitializeElement<FieldsTag, OptionsGroup, SourceTag>
,
161
detail::InitializeElement<FieldsTag, OptionsGroup, SubdomainOperator>>;
162
163
using
register_element = tmpl::list<
164
async_solvers::RegisterElement<FieldsTag, OptionsGroup, SourceTag>
,
165
detail::RegisterElement<FieldsTag, OptionsGroup, SourceTag>
>;
166
167
template
<
typename
ApplyOperatorActions,
typename
Label = OptionsGroup>
168
using
solve = tmpl::list<
169
async_solvers::PrepareSolve<FieldsTag, OptionsGroup, SourceTag, Label>
,
170
detail::SendOverlapData<FieldsTag, OptionsGroup, SubdomainOperator>
,
171
detail::ReceiveOverlapData<FieldsTag, OptionsGroup, SubdomainOperator>
,
172
detail::SolveSubdomain<FieldsTag, OptionsGroup, SubdomainOperator>,
173
detail::ReceiveOverlapSolution<FieldsTag, OptionsGroup,
174
SubdomainOperator
>,
175
ApplyOperatorActions,
176
async_solvers::CompleteStep<FieldsTag, OptionsGroup, SourceTag, Label>
>;
177
};
178
179
}
// 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:159
LinearSolver::Schwarz::SubdomainOperator
Abstract base class for the subdomain operator, i.e. the linear operator restricted to an element-cen...
Definition:
SubdomainOperator.hpp:43
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:126
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:147
LinearSolver::async_solvers::PrepareSolve
Definition:
ElementActions.hpp:184
TMPL.hpp
observers::Actions::RegisterWithObservers
Register an observation ID with the observers.
Definition:
RegisterWithObservers.hpp:47
LinearSolver::async_solvers::CompleteStep
Definition:
ElementActions.hpp:238
© Copyright 2017 - 2021
SXS Collaboration
,
Distributed under the
MIT License