Gmres.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include "IO/Observer/Actions/RegisterWithObservers.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, operand_tag>`.
29  *
30  * This linear solver supports preconditioning. Enable preconditioning by
31  * setting the `Preconditioned` template parameter to `true`. If you do, run a
32  * preconditioner (e.g. another parallel linear solver) in each step before
33  * invoking `perform_step`. The preconditioner should approximately solve the
34  * linear problem \f$A(q)=b\f$ where \f$q\f$ is the `operand_tag` and \f$b\f$ is
35  * the `preconditioner_source_tag`. Make sure the tag
36  * `db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, operand_tag>`
37  * is up-to-date with the preconditioned result before invoking `perform_step`,
38  * i.e. that it is \f$A(q)\f$ where \f$q\f$ is the preconditioner's approximate
39  * solution to \f$A(q)=b\f$.
40  *
41  * Note that the operand \f$q\f$ for which \f$A(q)\f$ needs to be computed is
42  * not the field \f$x\f$ we are solving for but
43  * `db::add_tag_prefix<LinearSolver::Tags::Operand, FieldsTag>`. This field is
44  * initially set to the residual \f$q_0 = b - A(x_0)\f$ where \f$x_0\f$ is the
45  * initial value of the `FieldsTag`.
46  *
47  * When the `perform_step` action is invoked after the operator action
48  * \f$A(q)\f$ has been computed and stored in the DataBox, the GMRES algorithm
49  * implemented here will converge the field \f$x\f$ towards the solution and
50  * update the operand \f$q\f$ in the process. This requires reductions over
51  * all elements that are received by a `ResidualMonitor` singleton parallel
52  * component, processed, and then broadcast back to all elements. Since the
53  * reductions are performed to find a vector that is orthogonal to those used in
54  * previous steps, the number of reductions increases linearly with iterations.
55  * No restarting mechanism is currently implemented. The actions are implemented
56  * in the `gmres::detail` namespace and constitute the full algorithm in the
57  * following order:
58  * 1. `PerformStep` (on elements): Start an Arnoldi orthogonalization by
59  * computing the inner product between \f$A(q)\f$ and the first of the
60  * previously determined set of orthogonal vectors.
61  * 2. `StoreOrthogonalization` (on `ResidualMonitor`): Keep track of the
62  * computed inner product in a Hessenberg matrix, then broadcast.
63  * 3. `OrthogonalizeOperand` (on elements): Proceed with the Arnoldi
64  * orthogonalization by computing inner products and reducing to
65  * `StoreOrthogonalization` on the `ResidualMonitor` until the new orthogonal
66  * vector is constructed. Then compute its magnitude and reduce.
67  * 4. `StoreFinalOrthogonalization` (on `ResidualMonitor`): Perform a QR
68  * decomposition of the Hessenberg matrix to produce a residual vector.
69  * Broadcast to `NormalizeOperandAndUpdateField` along with a termination
70  * flag if the `LinearSolver::Tags::ConvergenceCriteria` are met.
71  * 5. `NormalizeOperandAndUpdateField` (on elements): Set the operand \f$q\f$ as
72  * the new orthogonal vector and normalize. Use the residual vector and the set
73  * of orthogonal vectors to determine the solution \f$x\f$.
74  *
75  * \see ConjugateGradient for a linear solver that is more efficient when the
76  * linear operator \f$A\f$ is symmetric.
77  */
78 template <typename Metavariables, typename FieldsTag, typename OptionsGroup,
79  bool Preconditioned>
80 struct Gmres {
81  using fields_tag = FieldsTag;
82  using options_group = OptionsGroup;
83  static constexpr bool preconditioned = Preconditioned;
84 
85  /// Apply the linear operator to this tag in each iteration
87  Preconditioned,
92 
93  /// Invoke a linear solver on the `operand_tag` sourced by the
94  /// `preconditioner_source_tag` before applying the operator in each step
97 
98  /*!
99  * \brief The parallel components used by the GMRES linear solver
100  */
101  using component_list = tmpl::list<
102  detail::ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>;
103 
104  /*!
105  * \brief Initialize the tags used by the GMRES linear solver.
106  *
107  * Since we have not started iterating yet, we initialize the state _before_
108  * the first iteration. So `LinearSolver::Tags::IterationId` is undefined at
109  * this point and `Tags::Next<LinearSolver::Tags::IterationId>` is the initial
110  * step number. Invoke `prepare_step` to advance the state to the first
111  * iteration.
112  *
113  * \warning This action involves a blocking reduction, so it is a global
114  * synchronization point.
115  *
116  * With:
117  * - `initial_fields_tag` =
118  * `db::add_tag_prefix<LinearSolver::Tags::Initial, FieldsTag>`
119  * - `operand_tag` =
120  * `db::add_tag_prefix<LinearSolver::Tags::Operand, FieldsTag>`
121  * - `operator_tag` =
122  * `db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, operand_tag>`
123  * - `orthogonalization_iteration_id_tag` =
124  * `LinearSolver::Tags::Orthogonalization<LinearSolver::Tags::IterationId>`
125  * - `basis_history_tag` =
126  * `LinearSolver::Tags::KrylovSubspaceBasis<FieldsTag>`
127  *
128  * DataBox changes:
129  * - Adds:
130  * * `LinearSolver::Tags::IterationId`
131  * * `Tags::Next<LinearSolver::Tags::IterationId>`
132  * * `initial_fields_tag`
133  * * `orthogonalization_iteration_id_tag`
134  * * `basis_history_tag`
135  * * `LinearSolver::Tags::HasConverged`
136  * - Removes: nothing
137  * - Modifies:
138  * * `operand_tag`
139  *
140  * \note The `operand_tag` must already be present in the DataBox and is set
141  * to its initial value here. It is typically added to the DataBox by the
142  * system, which uses it to compute the `operator_tag` in each step. Also the
143  * `operator_tag` is typically added to the DataBox by the system, but does
144  * not need to be initialized until it is computed for the first time in the
145  * first step of the algorithm.
146  */
147  using initialize_element =
148  detail::InitializeElement<FieldsTag, OptionsGroup, Preconditioned>;
149 
150  using register_element = tmpl::list<>;
151 
152  /*!
153  * \brief Reset the linear solver to its initial state.
154  *
155  * Uses:
156  * - System:
157  * * `fields_tag`
158  *
159  * With:
160  * - `initial_fields_tag` =
161  * `db::add_tag_prefix<LinearSolver::Tags::Initial, fields_tag>`
162  * - `operand_tag` =
163  * `db::add_tag_prefix<LinearSolver::Tags::Operand, fields_tag>`
164  * - `orthogonalization_iteration_id_tag` =
165  * `LinearSolver::Tags::Orthogonalization<LinearSolver::Tags::IterationId>`
166  * - `basis_history_tag` =
167  * `LinearSolver::Tags::KrylovSubspaceBasis<fields_tag>`
168  *
169  * DataBox changes:
170  * - Adds: nothing
171  * - Removes: nothing
172  * - Modifies:
173  * * `LinearSolver::Tags::IterationId`
174  * * `initial_fields_tag`
175  * * `orthogonalization_iteration_id_tag`
176  * * `basis_history_tag`
177  * * `LinearSolver::Tags::HasConverged`
178  * * `operand_tag`
179  *
180  * \see `initialize_element`
181  */
182  using prepare_solve =
183  detail::PrepareSolve<FieldsTag, OptionsGroup, Preconditioned>;
184 
185  // Compile-time interface for observers
186  using observed_reduction_data_tags = observers::make_reduction_data_tags<
187  tmpl::list<observe_detail::reduction_data>>;
188 
189  /*!
190  * \brief Advance the linear solver to the next iteration.
191  *
192  * DataBox changes:
193  * - Adds: nothing
194  * - Removes: nothing
195  * - Modifies:
196  * * `LinearSolver::Tags::IterationId`
197  * * `Tags::Next<LinearSolver::Tags::IterationId>`
198  * * `orthogonalization_iteration_id_tag`
199  */
200  using prepare_step =
201  detail::PrepareStep<FieldsTag, OptionsGroup, Preconditioned>;
202 
203  /*!
204  * \brief Perform an iteration of the GMRES linear solver.
205  *
206  * \warning This action involves a blocking reduction, so it is a global
207  * synchronization point.
208  *
209  * With:
210  * - `operand_tag` =
211  * `db::add_tag_prefix<LinearSolver::Tags::Operand, FieldsTag>`
212  * - `orthogonalization_iteration_id_tag` =
213  * LinearSolver::Tags::Orthogonalization<LinearSolver::Tags::IterationId>`
214  * - `basis_history_tag` =
215  * `LinearSolver::Tags::KrylovSubspaceBasis<FieldsTag>`
216  *
217  * DataBox changes:
218  * - Adds: nothing
219  * - Removes: nothing
220  * - Modifies:
221  * * `FieldsTag`
222  * * `operand_tag`
223  * * `orthogonalization_iteration_id_tag`
224  * * `basis_history_tag`
225  * * `LinearSolver::Tags::HasConverged`
226  */
227  using perform_step =
228  detail::PerformStep<FieldsTag, OptionsGroup, Preconditioned>;
229 };
230 
231 } // namespace LinearSolver::gmres
LinearSolver::gmres::Gmres::prepare_solve
detail::PrepareSolve< FieldsTag, OptionsGroup, Preconditioned > prepare_solve
Reset the linear solver to its initial state.
Definition: Gmres.hpp:183
LinearSolver::gmres::Gmres
A GMRES solver for nonsymmetric linear systems of equations .
Definition: Gmres.hpp:80
LinearSolver::gmres::Gmres::preconditioner_source_tag
db::add_tag_prefix< LinearSolver::Tags::Operand, fields_tag > preconditioner_source_tag
Invoke a linear solver on the operand_tag sourced by the preconditioner_source_tag before applying th...
Definition: Gmres.hpp:96
db::add_tag_prefix
typename detail::add_tag_prefix_impl< Prefix, Tag, Args... >::type add_tag_prefix
Definition: PrefixHelpers.hpp:52
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:102
LinearSolver::gmres::Gmres::prepare_step
detail::PrepareStep< FieldsTag, OptionsGroup, Preconditioned > prepare_step
Advance the linear solver to the next iteration.
Definition: Gmres.hpp:201
LinearSolver::Tags::Preconditioned
Indicates the Tag is related to preconditioning of the linear solve.
Definition: Tags.hpp:224
LinearSolver::gmres::Gmres::perform_step
detail::PerformStep< FieldsTag, OptionsGroup, Preconditioned > perform_step
Perform an iteration of the GMRES linear solver.
Definition: Gmres.hpp:228
std::conditional_t
LinearSolver::gmres
Items related to the GMRES linear solver.
Definition: Gmres.cpp:16
LinearSolver::gmres::Gmres::initialize_element
detail::InitializeElement< FieldsTag, OptionsGroup, Preconditioned > initialize_element
Initialize the tags used by the GMRES linear solver.
Definition: Gmres.hpp:148
TMPL.hpp