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$. The preconditioner always begins at an initial
46  * guess of zero. It does not need to compute the operator applied to the
47  * initial guess, since it's zero as well due to the linearity of the operator.
48  *
49  * Note that the operand \f$q\f$ for which \f$A(q)\f$ needs to be computed is
50  * not the field \f$x\f$ we are solving for but
51  * `db::add_tag_prefix<LinearSolver::Tags::Operand, FieldsTag>`. This field is
52  * initially set to the residual \f$q_0 = b - A(x_0)\f$ where \f$x_0\f$ is the
53  * initial value of the `FieldsTag`.
54  *
55  * When the algorithm step is performed after the operator action \f$A(q)\f$ has
56  * been computed and stored in the DataBox, the GMRES algorithm implemented here
57  * will converge the field \f$x\f$ towards the solution and update the operand
58  * \f$q\f$ in the process. This requires reductions over all elements that are
59  * received by a `ResidualMonitor` singleton parallel component, processed, and
60  * then broadcast back to all elements. Since the reductions are performed to
61  * find a vector that is orthogonal to those used in previous steps, the number
62  * of reductions increases linearly with iterations. No restarting mechanism is
63  * currently implemented. The actions are implemented in the `gmres::detail`
64  * namespace and constitute the full algorithm in the following order:
65  * 1. `PerformStep` (on elements): Start an Arnoldi orthogonalization by
66  * computing the inner product between \f$A(q)\f$ and the first of the
67  * previously determined set of orthogonal vectors.
68  * 2. `StoreOrthogonalization` (on `ResidualMonitor`): Keep track of the
69  * computed inner product in a Hessenberg matrix, then broadcast.
70  * 3. `OrthogonalizeOperand` (on elements): Proceed with the Arnoldi
71  * orthogonalization by computing inner products and reducing to
72  * `StoreOrthogonalization` on the `ResidualMonitor` until the new orthogonal
73  * vector is constructed. Then compute its magnitude and reduce.
74  * 4. `StoreOrthogonalization` (on `ResidualMonitor`): Perform a QR
75  * decomposition of the Hessenberg matrix to produce a residual vector.
76  * Broadcast to `NormalizeOperandAndUpdateField` along with a termination
77  * flag if the `Convergence::Tags::Criteria` are met.
78  * 5. `NormalizeOperandAndUpdateField` (on elements): Set the operand \f$q\f$ as
79  * the new orthogonal vector and normalize. Use the residual vector and the set
80  * of orthogonal vectors to determine the solution \f$x\f$.
81  *
82  * \see ConjugateGradient for a linear solver that is more efficient when the
83  * linear operator \f$A\f$ is symmetric.
84  */
85 template <typename Metavariables, typename FieldsTag, typename OptionsGroup,
86  bool Preconditioned,
87  typename SourceTag =
89 struct Gmres {
90  using fields_tag = FieldsTag;
91  using options_group = OptionsGroup;
92  using source_tag = SourceTag;
93  static constexpr bool preconditioned = Preconditioned;
94 
95  /// Apply the linear operator to this tag in each iteration
97  Preconditioned,
102 
103  /// Invoke a linear solver on the `operand_tag` sourced by the
104  /// `preconditioner_source_tag` before applying the operator in each step
107 
108  /*!
109  * \brief The parallel components used by the GMRES linear solver
110  */
111  using component_list = tmpl::list<
112  detail::ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>;
113 
114  using initialize_element =
115  detail::InitializeElement<FieldsTag, OptionsGroup, Preconditioned>;
116 
117  using register_element = tmpl::list<>;
118 
119  using observed_reduction_data_tags = observers::make_reduction_data_tags<
120  tmpl::list<observe_detail::reduction_data>>;
121 
122  template <typename ApplyOperatorActions, typename Label = OptionsGroup>
123  using solve = tmpl::list<
124  detail::PrepareSolve<FieldsTag, OptionsGroup, Preconditioned, Label,
125  SourceTag>,
126  detail::NormalizeInitialOperand<FieldsTag, OptionsGroup, Preconditioned,
127  Label>,
128  detail::PrepareStep<FieldsTag, OptionsGroup, Preconditioned, Label>,
129  ApplyOperatorActions,
130  detail::PerformStep<FieldsTag, OptionsGroup, Preconditioned, Label>,
131  detail::OrthogonalizeOperand<FieldsTag, OptionsGroup, Preconditioned,
132  Label>,
133  detail::NormalizeOperandAndUpdateField<FieldsTag, OptionsGroup,
134  Preconditioned, Label>>;
135 };
136 
137 } // namespace LinearSolver::gmres
LinearSolver::gmres::Gmres
A GMRES solver for nonsymmetric linear systems of equations .
Definition: Gmres.hpp:89
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:112
Tags.hpp
LinearSolver::Tags::Preconditioned
Indicates the Tag is related to preconditioning of the linear solve.
Definition: Tags.hpp:170
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:106
LinearSolver::gmres
Items related to the GMRES linear solver.
Definition: Gmres.hpp:29
Prefixes.hpp
TMPL.hpp