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