Gmres.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include "DataStructures/DataBox/PrefixHelpers.hpp"
8 #include "IO/Observer/Helpers.hpp"
9 #include "ParallelAlgorithms/LinearSolver/Gmres/ElementActions.hpp"
10 #include "ParallelAlgorithms/LinearSolver/Gmres/InitializeElement.hpp"
11 #include "ParallelAlgorithms/LinearSolver/Gmres/ResidualMonitor.hpp"
12 #include "ParallelAlgorithms/LinearSolver/Observe.hpp"
14 #include "Utilities/TMPL.hpp"
15 
16 /// Items related to the GMRES linear solver
17 ///
18 /// \see `LinearSolver::gmres::Gmres`
19 namespace LinearSolver::gmres {
20 
21 /*!
22  * \ingroup LinearSolverGroup
23  * \brief A GMRES solver for nonsymmetric linear systems of equations
24  * \f$Ax=b\f$.
25  *
26  * \details The only operation we need to supply to the algorithm is the
27  * result of the operation \f$A(p)\f$ (see \ref LinearSolverGroup). Each step of
28  * the algorithm expects that \f$A(q)\f$ is computed and stored in the DataBox
29  * as `db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, operand_tag>`.
30  * To perform a solve, add the `solve` action list to an array parallel
31  * component. Pass the actions that compute \f$A(q)\f$, as well as any further
32  * actions you wish to run in each step of the algorithm, as the first template
33  * parameter to `solve`. If you add the `solve` action list multiple times, use
34  * the second template parameter to label each solve with a different type.
35  *
36  * This linear solver supports preconditioning. Enable preconditioning by
37  * setting the `Preconditioned` template parameter to `true`. If you do, run a
38  * preconditioner (e.g. another parallel linear solver) in each step. The
39  * preconditioner should approximately solve the linear problem \f$A(q)=b\f$
40  * where \f$q\f$ is the `operand_tag` and \f$b\f$ is the
41  * `preconditioner_source_tag`. Make sure the tag
42  * `db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, operand_tag>`
43  * is updated with the preconditioned result in each step of the algorithm, i.e.
44  * that it is \f$A(q)\f$ where \f$q\f$ is the preconditioner's approximate
45  * solution to \f$A(q)=b\f$.
46  *
47  * Note that the operand \f$q\f$ for which \f$A(q)\f$ needs to be computed is
48  * not the field \f$x\f$ we are solving for but
49  * `db::add_tag_prefix<LinearSolver::Tags::Operand, FieldsTag>`. This field is
50  * initially set to the residual \f$q_0 = b - A(x_0)\f$ where \f$x_0\f$ is the
51  * initial value of the `FieldsTag`.
52  *
53  * When the algorithm step is performed after the operator action \f$A(q)\f$ has
54  * been computed and stored in the DataBox, the GMRES algorithm implemented here
55  * will converge the field \f$x\f$ towards the solution and update the operand
56  * \f$q\f$ in the process. This requires reductions over all elements that are
57  * received by a `ResidualMonitor` singleton parallel component, processed, and
58  * then broadcast back to all elements. Since the reductions are performed to
59  * find a vector that is orthogonal to those used in previous steps, the number
60  * of reductions increases linearly with iterations. No restarting mechanism is
61  * currently implemented. The actions are implemented in the `gmres::detail`
62  * namespace and constitute the full algorithm in the following order:
63  * 1. `PerformStep` (on elements): Start an Arnoldi orthogonalization by
64  * computing the inner product between \f$A(q)\f$ and the first of the
65  * previously determined set of orthogonal vectors.
66  * 2. `StoreOrthogonalization` (on `ResidualMonitor`): Keep track of the
67  * computed inner product in a Hessenberg matrix, then broadcast.
68  * 3. `OrthogonalizeOperand` (on elements): Proceed with the Arnoldi
69  * orthogonalization by computing inner products and reducing to
70  * `StoreOrthogonalization` on the `ResidualMonitor` until the new orthogonal
71  * vector is constructed. Then compute its magnitude and reduce.
72  * 4. `StoreOrthogonalization` (on `ResidualMonitor`): Perform a QR
73  * decomposition of the Hessenberg matrix to produce a residual vector.
74  * Broadcast to `NormalizeOperandAndUpdateField` along with a termination
75  * flag if the `Convergence::Tags::Criteria` are met.
76  * 5. `NormalizeOperandAndUpdateField` (on elements): Set the operand \f$q\f$ as
77  * the new orthogonal vector and normalize. Use the residual vector and the set
78  * of orthogonal vectors to determine the solution \f$x\f$.
79  *
80  * \see ConjugateGradient for a linear solver that is more efficient when the
81  * linear operator \f$A\f$ is symmetric.
82  */
83 template <typename Metavariables, typename FieldsTag, typename OptionsGroup,
84  bool Preconditioned,
85  typename SourceTag =
87 struct Gmres {
88  using fields_tag = FieldsTag;
89  using options_group = OptionsGroup;
90  using source_tag = SourceTag;
91  static constexpr bool preconditioned = Preconditioned;
92 
93  /// Apply the linear operator to this tag in each iteration
95  Preconditioned,
100 
101  /// Invoke a linear solver on the `operand_tag` sourced by the
102  /// `preconditioner_source_tag` before applying the operator in each step
105 
106  /*!
107  * \brief The parallel components used by the GMRES linear solver
108  */
109  using component_list = tmpl::list<
110  detail::ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>;
111 
112  using initialize_element =
113  detail::InitializeElement<FieldsTag, OptionsGroup, Preconditioned>;
114 
115  using register_element = tmpl::list<>;
116 
117  using observed_reduction_data_tags = observers::make_reduction_data_tags<
118  tmpl::list<observe_detail::reduction_data>>;
119 
120  template <typename ApplyOperatorActions, typename Label = OptionsGroup>
121  using solve = tmpl::list<
122  detail::PrepareSolve<FieldsTag, OptionsGroup, Preconditioned, Label,
123  SourceTag>,
124  detail::NormalizeInitialOperand<FieldsTag, OptionsGroup, Preconditioned,
125  Label>,
126  detail::PrepareStep<FieldsTag, OptionsGroup, Preconditioned, Label>,
127  ApplyOperatorActions,
128  detail::PerformStep<FieldsTag, OptionsGroup, Preconditioned, Label>,
129  detail::OrthogonalizeOperand<FieldsTag, OptionsGroup, Preconditioned,
130  Label>,
131  detail::NormalizeOperandAndUpdateField<FieldsTag, OptionsGroup,
132  Preconditioned, Label>>;
133 };
134 
135 } // namespace LinearSolver::gmres
LinearSolver::gmres::Gmres
A GMRES solver for nonsymmetric linear systems of equations .
Definition: Gmres.hpp:87
db::add_tag_prefix
typename detail::add_tag_prefix_impl< Prefix, Tag, Args... >::type add_tag_prefix
Definition: PrefixHelpers.hpp:51
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:110
Tags.hpp
LinearSolver::Tags::Preconditioned
Indicates the Tag is related to preconditioning of the linear solve.
Definition: Tags.hpp:171
std::conditional_t
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:104
LinearSolver::gmres
Items related to the GMRES linear solver.
Definition: Gmres.hpp:29
Prefixes.hpp
TMPL.hpp