Gmres.hpp
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:
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:
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:
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:
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
Invoke a linear solver on the operand_tag sourced by the preconditioner_source_tag before applying th...
Definition: Gmres.hpp:96
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