Gmres.hpp
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
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:
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:
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