Gmres.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 "NumericalAlgorithms/LinearSolver/Gmres/ElementActions.hpp"
8 #include "NumericalAlgorithms/LinearSolver/Gmres/InitializeElement.hpp"
9 #include "NumericalAlgorithms/LinearSolver/Gmres/ResidualMonitor.hpp"
10 #include "Utilities/TMPL.hpp"
11 
12 namespace LinearSolver {
13 
14 /*!
15  * \ingroup LinearSolverGroup
16  * \brief A GMRES solver for nonsymmetric linear systems of equations
17  * \f$Ax=b\f$.
18  *
19  * \details The only operation we need to supply to the algorithm is the
20  * result of the operation \f$A(p)\f$ (see \ref LinearSolverGroup). Each
21  * invocation of the `perform_step` action expects that \f$A(q)\f$ has been
22  * computed in a preceding action and stored in the DataBox as
23  * %db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo,
24  * db::add_tag_prefix<LinearSolver::Tags::Operand, typename
25  * Metavariables::system::fields_tag>>.
26  *
27  * Note that the operand \f$q\f$ for which \f$A(q)\f$ needs to be computed is
28  * not the field \f$x\f$ we are solving for but
29  * `db::add_tag_prefix<LinearSolver::Tags::Operand, typename
30  * Metavariables::system::fields_tag>`. This field is initially set to the
31  * residual \f$q_0 = b - A(x_0)\f$ where \f$x_0\f$ is the initial value of the
32  * `Metavariables::system::fields_tag`.
33  *
34  * When the `perform_step` action is invoked after the operator action
35  * \f$A(q)\f$ has been computed and stored in the DataBox, the GMRES algorithm
36  * implemented here will converge the field \f$x\f$ towards the solution and
37  * update the operand \f$q\f$ in the process. This requires reductions over
38  * all elements that are received by a `ResidualMonitor` singleton parallel
39  * component, processed, and then broadcast back to all elements. Since the
40  * reductions are performed to find a vector that is orthogonal to those used in
41  * previous steps, the number of reductions increases linearly with iterations.
42  * No restarting mechanism is currently implemented. The actions are implemented
43  * in the `gmres_detail` namespace and constitute the full algorithm in the
44  * following order:
45  * 1. `PerformStep` (on elements): Start an Arnoldi orthogonalization by
46  * computing the inner product between \f$A(q)\f$ and the first of the
47  * previously determined set of orthogonal vectors.
48  * 2. `StoreOrthogonalization` (on `ResidualMonitor`): Keep track of the
49  * computed inner product in a Hessenberg matrix, then broadcast.
50  * 3. `OrthogonalizeOperand` (on elements): Proceed with the Arnoldi
51  * orthogonalization by computing inner products and reducing to
52  * `StoreOrthogonalization` on the `ResidualMonitor` until the new orthogonal
53  * vector is constructed. Then compute its magnitude and reduce.
54  * 4. `StoreFinalOrthogonalization` (on `ResidualMonitor`): Perform a QR
55  * decomposition of the Hessenberg matrix to produce a residual vector.
56  * Broadcast to `NormalizeOperandAndUpdateField` along with a termination
57  * flag if the `LinearSolver::Tags::ConvergenceCriteria` are met.
58  * 5. `NormalizeOperandAndUpdateField` (on elements): Set the operand \f$q\f$ as
59  * the new orthogonal vector and normalize. Use the residual vector and the set
60  * of orthogonal vectors to determine the solution \f$x\f$.
61  *
62  * \see ConjugateGradient for a linear solver that is more efficient when the
63  * linear operator \f$A\f$ is symmetric.
64  */
65 template <typename Metavariables>
66 struct Gmres {
67  /*!
68  * \brief The parallel components used by the GMRES linear solver
69  *
70  * Uses:
71  * - System:
72  * * `fields_tag`
73  */
74  using component_list =
75  tmpl::list<gmres_detail::ResidualMonitor<Metavariables>>;
76 
77  /*!
78  * \brief Initialize the tags used by the GMRES linear solver
79  *
80  * Uses:
81  * - System:
82  * * `fields_tag`
83  * - ConstGlobalCache: nothing
84  *
85  * With:
86  * - `initial_fields_tag` =
87  * `db::add_tag_prefix<LinearSolver::Tags::Initial, fields_tag>`
88  * - `operand_tag` =
89  * `db::add_tag_prefix<LinearSolver::Tags::Operand, fields_tag>`
90  * - `operator_tag` =
91  * `db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, operand_tag>`
92  * - `orthogonalization_iteration_id_tag` =
93  * `db::add_tag_prefix<LinearSolver::Tags::Orthogonalization,
94  * LinearSolver::Tags::IterationId>`
95  * - `basis_history_tag` =
96  * `LinearSolver::Tags::KrylovSubspaceBasis<fields_tag>`
97  *
98  * DataBox changes:
99  * - Adds:
100  * * `LinearSolver::Tags::IterationId`
101  * * `Tags::Next<LinearSolver::Tags::IterationId>`
102  * * `initial_fields_tag`
103  * * `orthogonalization_iteration_id_tag`
104  * * `basis_history_tag`
105  * * `LinearSolver::Tags::HasConverged`
106  * - Removes: nothing
107  * - Modifies:
108  * * `operand_tag`
109  *
110  * \note The `operand_tag` must already be present in the DataBox and is set
111  * to its initial value here. It is typically added to the DataBox by the
112  * system, which uses it to compute the `operator_tag` in each step. Also the
113  * `operator_tag` is typically added to the DataBox by the system, but does
114  * not need to be initialized until it is computed for the first time in the
115  * first step of the algorithm.
116  */
117  using tags = gmres_detail::InitializeElement<Metavariables>;
118 
119  // Compile-time interface for observers
120  using observed_reduction_data_tags = observers::make_reduction_data_tags<
121  tmpl::list<observe_detail::reduction_data>>;
122 
123  /*!
124  * \brief Perform an iteration of the GMRES linear solver
125  *
126  * Uses:
127  * - System:
128  * * `fields_tag`
129  * - ConstGlobalCache: nothing
130  *
131  * With:
132  * - `operand_tag` =
133  * `db::add_tag_prefix<LinearSolver::Tags::Operand, fields_tag>`
134  * - `orthogonalization_iteration_id_tag` =
135  * `db::add_tag_prefix<LinearSolver::Tags::Orthogonalization,
136  * LinearSolver::Tags::IterationId>`
137  * - `basis_history_tag` =
138  * `LinearSolver::Tags::KrylovSubspaceBasis<fields_tag>`
139  *
140  * DataBox changes:
141  * - Adds: nothing
142  * - Removes: nothing
143  * - Modifies:
144  * * `LinearSolver::Tags::IterationId`
145  * * `Tags::Next<LinearSolver::Tags::IterationId>`
146  * * `fields_tag`
147  * * `operand_tag`
148  * * `orthogonalization_iteration_id_tag`
149  * * `basis_history_tag`
150  * * `LinearSolver::Tags::HasConverged`
151  */
152  using perform_step = gmres_detail::PerformStep;
153 };
154 
155 } // namespace LinearSolver
A GMRES solver for nonsymmetric linear systems of equations .
Definition: Gmres.hpp:66
Functionality for solving linear systems of equations.
Definition: TerminateIfConverged.hpp:22
gmres_detail::PerformStep perform_step
Perform an iteration of the GMRES linear solver.
Definition: Gmres.hpp:152
gmres_detail::InitializeElement< Metavariables > tags
Initialize the tags used by the GMRES linear solver.
Definition: Gmres.hpp:117
Wraps the template metaprogramming library used (brigand)
tmpl::list< gmres_detail::ResidualMonitor< Metavariables > > component_list
The parallel components used by the GMRES linear solver.
Definition: Gmres.hpp:75