Gmres.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include "IO/Observer/Actions.hpp"
7 #include "IO/Observer/Helpers.hpp"
8 #include "ParallelAlgorithms/Initialization/MergeIntoDataBox.hpp"
9 #include "ParallelAlgorithms/LinearSolver/Gmres/ElementActions.hpp"
10 #include "ParallelAlgorithms/LinearSolver/Gmres/InitializeElement.hpp"
11 #include "ParallelAlgorithms/LinearSolver/Gmres/ResidualMonitor.hpp"
12 #include "Utilities/TMPL.hpp"
13 
14 /// Items related to the GMRES linear solver
15 ///
16 /// \see `LinearSolver::gmres::Gmres`
17 namespace LinearSolver::gmres {
18 
19 /*!
20  * \ingroup LinearSolverGroup
21  * \brief A GMRES solver for nonsymmetric linear systems of equations
22  * \f$Ax=b\f$.
23  *
24  * \details The only operation we need to supply to the algorithm is the
25  * result of the operation \f$A(p)\f$ (see \ref LinearSolverGroup). Each
26  * invocation of the `perform_step` action expects that \f$A(q)\f$ has been
27  * computed in a preceding action and stored in the DataBox as
28  * %db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo,
29  * db::add_tag_prefix<LinearSolver::Tags::Operand, FieldsTag>>.
30  *
31  * Note that the operand \f$q\f$ for which \f$A(q)\f$ needs to be computed is
32  * not the field \f$x\f$ we are solving for but
33  * `db::add_tag_prefix<LinearSolver::Tags::Operand, FieldsTag>`. This field is
34  * initially set to the residual \f$q_0 = b - A(x_0)\f$ where \f$x_0\f$ is the
35  * initial value of the `FieldsTag`.
36  *
37  * When the `perform_step` action is invoked after the operator action
38  * \f$A(q)\f$ has been computed and stored in the DataBox, the GMRES algorithm
39  * implemented here will converge the field \f$x\f$ towards the solution and
40  * update the operand \f$q\f$ in the process. This requires reductions over
41  * all elements that are received by a `ResidualMonitor` singleton parallel
42  * component, processed, and then broadcast back to all elements. Since the
43  * reductions are performed to find a vector that is orthogonal to those used in
44  * previous steps, the number of reductions increases linearly with iterations.
45  * No restarting mechanism is currently implemented. The actions are implemented
46  * in the `gmres::detail` namespace and constitute the full algorithm in the
47  * following order:
48  * 1. `PerformStep` (on elements): Start an Arnoldi orthogonalization by
49  * computing the inner product between \f$A(q)\f$ and the first of the
50  * previously determined set of orthogonal vectors.
51  * 2. `StoreOrthogonalization` (on `ResidualMonitor`): Keep track of the
52  * computed inner product in a Hessenberg matrix, then broadcast.
53  * 3. `OrthogonalizeOperand` (on elements): Proceed with the Arnoldi
54  * orthogonalization by computing inner products and reducing to
55  * `StoreOrthogonalization` on the `ResidualMonitor` until the new orthogonal
56  * vector is constructed. Then compute its magnitude and reduce.
57  * 4. `StoreFinalOrthogonalization` (on `ResidualMonitor`): Perform a QR
58  * decomposition of the Hessenberg matrix to produce a residual vector.
59  * Broadcast to `NormalizeOperandAndUpdateField` along with a termination
60  * flag if the `LinearSolver::Tags::ConvergenceCriteria` are met.
61  * 5. `NormalizeOperandAndUpdateField` (on elements): Set the operand \f$q\f$ as
62  * the new orthogonal vector and normalize. Use the residual vector and the set
63  * of orthogonal vectors to determine the solution \f$x\f$.
64  *
65  * \see ConjugateGradient for a linear solver that is more efficient when the
66  * linear operator \f$A\f$ is symmetric.
67  */
68 template <typename Metavariables, typename FieldsTag, typename OptionsGroup>
69 struct Gmres {
70  using fields_tag = FieldsTag;
71  using options_group = OptionsGroup;
72 
73  /// Apply the linear operator to this tag in each iteration
74  using operand_tag =
76 
77  /*!
78  * \brief The parallel components used by the GMRES linear solver
79  */
80  using component_list = tmpl::list<
81  detail::ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>;
82 
83  /*!
84  * \brief Initialize the tags used by the GMRES linear solver.
85  *
86  * Since we have not started iterating yet, we initialize the state _before_
87  * the first iteration. So `LinearSolver::Tags::IterationId` is undefined at
88  * this point and `Tags::Next<LinearSolver::Tags::IterationId>` is the initial
89  * step number. Invoke `prepare_step` to advance the state to the first
90  * iteration.
91  *
92  * \warning This action involves a blocking reduction, so it is a global
93  * synchronization point.
94  *
95  * With:
96  * - `initial_fields_tag` =
97  * `db::add_tag_prefix<LinearSolver::Tags::Initial, FieldsTag>`
98  * - `operand_tag` =
99  * `db::add_tag_prefix<LinearSolver::Tags::Operand, FieldsTag>`
100  * - `operator_tag` =
101  * `db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, operand_tag>`
102  * - `orthogonalization_iteration_id_tag` =
103  * `db::add_tag_prefix<LinearSolver::Tags::Orthogonalization,
104  * LinearSolver::Tags::IterationId>`
105  * - `basis_history_tag` =
106  * `LinearSolver::Tags::KrylovSubspaceBasis<FieldsTag>`
107  *
108  * DataBox changes:
109  * - Adds:
110  * * `LinearSolver::Tags::IterationId`
111  * * `Tags::Next<LinearSolver::Tags::IterationId>`
112  * * `initial_fields_tag`
113  * * `orthogonalization_iteration_id_tag`
114  * * `basis_history_tag`
115  * * `LinearSolver::Tags::HasConverged`
116  * - Removes: nothing
117  * - Modifies:
118  * * `operand_tag`
119  *
120  * \note The `operand_tag` must already be present in the DataBox and is set
121  * to its initial value here. It is typically added to the DataBox by the
122  * system, which uses it to compute the `operator_tag` in each step. Also the
123  * `operator_tag` is typically added to the DataBox by the system, but does
124  * not need to be initialized until it is computed for the first time in the
125  * first step of the algorithm.
126  */
127  using initialize_element = detail::InitializeElement<FieldsTag, OptionsGroup>;
128 
129  using register_element =
131 
132  /*!
133  * \brief Reset the linear solver to its initial state.
134  *
135  * Uses:
136  * - System:
137  * * `fields_tag`
138  *
139  * With:
140  * - `initial_fields_tag` =
141  * `db::add_tag_prefix<LinearSolver::Tags::Initial, fields_tag>`
142  * - `operand_tag` =
143  * `db::add_tag_prefix<LinearSolver::Tags::Operand, fields_tag>`
144  * - `orthogonalization_iteration_id_tag` =
145  * `db::add_tag_prefix<LinearSolver::Tags::Orthogonalization,
146  * LinearSolver::Tags::IterationId>`
147  * - `basis_history_tag` =
148  * `LinearSolver::Tags::KrylovSubspaceBasis<fields_tag>`
149  *
150  * DataBox changes:
151  * - Adds: nothing
152  * - Removes: nothing
153  * - Modifies:
154  * * `LinearSolver::Tags::IterationId`
155  * * `initial_fields_tag`
156  * * `orthogonalization_iteration_id_tag`
157  * * `basis_history_tag`
158  * * `LinearSolver::Tags::HasConverged`
159  * * `operand_tag`
160  *
161  * \see `initialize_element`
162  */
163  using prepare_solve = detail::PrepareSolve<FieldsTag, OptionsGroup>;
164 
165  // Compile-time interface for observers
166  using observed_reduction_data_tags = observers::make_reduction_data_tags<
167  tmpl::list<observe_detail::reduction_data>>;
168 
169  /*!
170  * \brief Advance the linear solver to the next iteration.
171  *
172  * DataBox changes:
173  * - Adds: nothing
174  * - Removes: nothing
175  * - Modifies:
176  * * `LinearSolver::Tags::IterationId`
177  * * `Tags::Next<LinearSolver::Tags::IterationId>`
178  * * `orthogonalization_iteration_id_tag`
179  */
180  using prepare_step = detail::PrepareStep<FieldsTag, OptionsGroup>;
181 
182  /*!
183  * \brief Perform an iteration of the GMRES linear solver.
184  *
185  * \warning This action involves a blocking reduction, so it is a global
186  * synchronization point.
187  *
188  * With:
189  * - `operand_tag` =
190  * `db::add_tag_prefix<LinearSolver::Tags::Operand, FieldsTag>`
191  * - `orthogonalization_iteration_id_tag` =
192  * `db::add_tag_prefix<LinearSolver::Tags::Orthogonalization,
193  * LinearSolver::Tags::IterationId>`
194  * - `basis_history_tag` =
195  * `LinearSolver::Tags::KrylovSubspaceBasis<FieldsTag>`
196  *
197  * DataBox changes:
198  * - Adds: nothing
199  * - Removes: nothing
200  * - Modifies:
201  * * `FieldsTag`
202  * * `operand_tag`
203  * * `orthogonalization_iteration_id_tag`
204  * * `basis_history_tag`
205  * * `LinearSolver::Tags::HasConverged`
206  */
207  using perform_step = detail::PerformStep<FieldsTag, OptionsGroup>;
208 };
209 
210 } // namespace LinearSolver::gmres
LinearSolver::gmres::Gmres::prepare_solve
detail::PrepareSolve< FieldsTag, OptionsGroup > prepare_solve
Reset the linear solver to its initial state.
Definition: Gmres.hpp:163
LinearSolver::gmres::Gmres
A GMRES solver for nonsymmetric linear systems of equations .
Definition: Gmres.hpp:69
LinearSolver::gmres::Gmres::initialize_element
detail::InitializeElement< FieldsTag, OptionsGroup > initialize_element
Initialize the tags used by the GMRES linear solver.
Definition: Gmres.hpp:127
LinearSolver::gmres::Gmres::perform_step
detail::PerformStep< FieldsTag, OptionsGroup > perform_step
Perform an iteration of the GMRES linear solver.
Definition: Gmres.hpp:207
LinearSolver::gmres::Gmres::operand_tag
db::add_tag_prefix< LinearSolver::Tags::Operand, fields_tag > operand_tag
Apply the linear operator to this tag in each iteration.
Definition: Gmres.hpp:75
db::add_tag_prefix
Prefix< DataBox_detail::dispatch_add_tag_prefix_impl< Prefix, Tag, Args... >, Args... > add_tag_prefix
Definition: PrefixHelpers.hpp:145
LinearSolver::gmres::Gmres::prepare_step
detail::PrepareStep< FieldsTag, OptionsGroup > prepare_step
Advance the linear solver to the next iteration.
Definition: Gmres.hpp:180
LinearSolver::gmres::Gmres::component_list
tmpl::list< detail::ResidualMonitor< Metavariables, FieldsTag, OptionsGroup > > component_list
The parallel components used by the GMRES linear solver.
Definition: Gmres.hpp:81
LinearSolver::gmres
Items related to the GMRES linear solver.
Definition: Gmres.cpp:16
TMPL.hpp
observers::Actions::RegisterWithObservers
Registers itself with the local observer parallel component so the observer knows to expect data from...
Definition: Actions.hpp:232