Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
7 : #include "DataStructures/DataBox/Prefixes.hpp"
8 : #include "ParallelAlgorithms/LinearSolver/Gmres/ElementActions.hpp"
9 : #include "ParallelAlgorithms/LinearSolver/Gmres/InitializeElement.hpp"
10 : #include "ParallelAlgorithms/LinearSolver/Gmres/ResidualMonitor.hpp"
11 : #include "ParallelAlgorithms/LinearSolver/Tags.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 step of
26 : * the algorithm expects that \f$A(q)\f$ is computed and stored in the DataBox
27 : * as `db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, operand_tag>`.
28 : * To perform a solve, add the `solve` action list to an array parallel
29 : * component. Pass the actions that compute \f$A(q)\f$ as the first template
30 : * parameter to `solve`. The second template parameter allows you to pass any
31 : * further actions you wish to run in each step of the algorithm (such as
32 : * observations). If you add the `solve` action list multiple times, use the
33 : * third template parameter to label each solve with a different type.
34 : *
35 : * This linear solver supports preconditioning. Enable preconditioning by
36 : * setting the `Preconditioned` template parameter to `true`. If you do, run a
37 : * preconditioner (e.g. another parallel linear solver) in each step. The
38 : * preconditioner should approximately solve the linear problem \f$A(q)=b\f$
39 : * where \f$q\f$ is the `operand_tag` and \f$b\f$ is the
40 : * `preconditioner_source_tag`. Make sure the tag
41 : * `db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, operand_tag>`
42 : * is updated with the preconditioned result in each step of the algorithm, i.e.
43 : * that it is \f$A(q)\f$ where \f$q\f$ is the preconditioner's approximate
44 : * solution to \f$A(q)=b\f$. The preconditioner always begins at an initial
45 : * guess of zero. It does not need to compute the operator applied to the
46 : * initial guess, since it's zero as well due to the linearity of the operator.
47 : *
48 : * Note that the operand \f$q\f$ for which \f$A(q)\f$ needs to be computed is
49 : * not the field \f$x\f$ we are solving for but
50 : * `db::add_tag_prefix<LinearSolver::Tags::Operand, FieldsTag>`. This field is
51 : * initially set to the residual \f$q_0 = b - A(x_0)\f$ where \f$x_0\f$ is the
52 : * initial value of the `FieldsTag`.
53 : *
54 : * When the algorithm step is performed after the operator action \f$A(q)\f$ has
55 : * been computed and stored in the DataBox, the GMRES algorithm implemented here
56 : * will converge the field \f$x\f$ towards the solution and update the operand
57 : * \f$q\f$ in the process. This requires reductions over all elements that are
58 : * received by a `ResidualMonitor` singleton parallel component, processed, and
59 : * then broadcast back to all elements. Since the reductions are performed to
60 : * find a vector that is orthogonal to those used in previous steps, the number
61 : * of reductions increases linearly with iterations. No restarting mechanism is
62 : * currently implemented. The actions are implemented in the `gmres::detail`
63 : * namespace and constitute the full algorithm in the following order:
64 : * 1. `PerformStep` (on elements): Start an Arnoldi orthogonalization by
65 : * computing the inner product between \f$A(q)\f$ and the first of the
66 : * previously determined set of orthogonal vectors.
67 : * 2. `StoreOrthogonalization` (on `ResidualMonitor`): Keep track of the
68 : * computed inner product in a Hessenberg matrix, then broadcast.
69 : * 3. `OrthogonalizeOperand` (on elements): Proceed with the Arnoldi
70 : * orthogonalization by computing inner products and reducing to
71 : * `StoreOrthogonalization` on the `ResidualMonitor` until the new orthogonal
72 : * vector is constructed. Then compute its magnitude and reduce.
73 : * 4. `StoreOrthogonalization` (on `ResidualMonitor`): Perform a QR
74 : * decomposition of the Hessenberg matrix to produce a residual vector.
75 : * Broadcast to `NormalizeOperandAndUpdateField` along with a termination
76 : * flag if the `Convergence::Tags::Criteria` are met.
77 : * 5. `NormalizeOperandAndUpdateField` (on elements): Set the operand \f$q\f$ as
78 : * the new orthogonal vector and normalize. Use the residual vector and the set
79 : * of orthogonal vectors to determine the solution \f$x\f$.
80 : *
81 : * \par Array sections
82 : * This linear solver supports running over a subset of the elements in the
83 : * array parallel component (see `Parallel::Section`). Set the
84 : * `ArraySectionIdTag` template parameter to restrict the solver to elements in
85 : * that section. Only a single section must be associated with the
86 : * `ArraySectionIdTag`. The default is `void`, which means running over all
87 : * elements in the array. Note that the actions in the `ApplyOperatorActions`
88 : * list passed to `solve` will _not_ be restricted to run only on section
89 : * elements, so all elements in the array may participate in preconditioning
90 : * (see LinearSolver::multigrid::Multigrid).
91 : *
92 : * \see ConjugateGradient for a linear solver that is more efficient when the
93 : * linear operator \f$A\f$ is symmetric.
94 : */
95 : template <typename Metavariables, typename FieldsTag, typename OptionsGroup,
96 : bool Preconditioned,
97 : typename SourceTag =
98 : db::add_tag_prefix<::Tags::FixedSource, FieldsTag>,
99 : typename ArraySectionIdTag = void>
100 1 : struct Gmres {
101 0 : using fields_tag = FieldsTag;
102 0 : using options_group = OptionsGroup;
103 0 : using source_tag = SourceTag;
104 0 : static constexpr bool preconditioned = Preconditioned;
105 :
106 : /// Apply the linear operator to this tag in each iteration
107 1 : using operand_tag = std::conditional_t<
108 : Preconditioned,
109 : db::add_tag_prefix<
110 : LinearSolver::Tags::Preconditioned,
111 : db::add_tag_prefix<LinearSolver::Tags::Operand, fields_tag>>,
112 : db::add_tag_prefix<LinearSolver::Tags::Operand, fields_tag>>;
113 :
114 : /// Invoke a linear solver on the `operand_tag` sourced by the
115 : /// `preconditioner_source_tag` before applying the operator in each step
116 1 : using preconditioner_source_tag =
117 : db::add_tag_prefix<LinearSolver::Tags::Operand, fields_tag>;
118 :
119 : /*!
120 : * \brief The parallel components used by the GMRES linear solver
121 : */
122 1 : using component_list = tmpl::list<
123 : detail::ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>;
124 :
125 0 : using initialize_element =
126 : detail::InitializeElement<FieldsTag, OptionsGroup, Preconditioned>;
127 :
128 0 : using register_element = tmpl::list<>;
129 :
130 0 : using amr_projectors = initialize_element;
131 :
132 : template <typename ApplyOperatorActions,
133 : typename ObserveActions = tmpl::list<>,
134 : typename Label = OptionsGroup>
135 0 : using solve = tmpl::list<
136 : detail::PrepareSolve<FieldsTag, OptionsGroup, Preconditioned, Label,
137 : SourceTag, ArraySectionIdTag>,
138 : ObserveActions,
139 : detail::NormalizeInitialOperand<FieldsTag, OptionsGroup, Preconditioned,
140 : Label, ArraySectionIdTag>,
141 : detail::PrepareStep<FieldsTag, OptionsGroup, Preconditioned, Label,
142 : ArraySectionIdTag>,
143 : ApplyOperatorActions,
144 : detail::PerformStep<FieldsTag, OptionsGroup, Preconditioned, Label,
145 : ArraySectionIdTag>,
146 : detail::OrthogonalizeOperand<FieldsTag, OptionsGroup, Preconditioned,
147 : Label, ArraySectionIdTag>,
148 : detail::NormalizeOperandAndUpdateField<
149 : FieldsTag, OptionsGroup, Preconditioned, Label, ArraySectionIdTag>,
150 : ObserveActions,
151 : detail::CompleteStep<FieldsTag, OptionsGroup, Preconditioned, Label,
152 : ArraySectionIdTag>>;
153 : };
154 :
155 : } // namespace LinearSolver::gmres
|
