ConjugateGradient.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include "NumericalAlgorithms/LinearSolver/ConjugateGradient/ElementActions.hpp"
7 #include "NumericalAlgorithms/LinearSolver/ConjugateGradient/InitializeElement.hpp"
8 #include "NumericalAlgorithms/LinearSolver/ConjugateGradient/ResidualMonitor.hpp"
9 #include "Utilities/TMPL.hpp"
10 
11 namespace LinearSolver {
12 
13 /*!
14  * \ingroup LinearSolverGroup
15  * \brief A conjugate gradient solver for linear systems of equations \f$Ax=b\f$
16  * where the operator \f$A\f$ is symmetric.
17  *
18  * \details The only operation we need to supply to the algorithm is the
19  * result of the operation \f$A(p)\f$ (see \ref LinearSolverGroup) that in the
20  * case of the conjugate gradient algorithm must be symmetric. Each invocation
21  * of the `perform_step` action expects that \f$A(p)\f$ has been computed in a
22  * preceding action and stored in the DataBox as
23  * %db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo,
24  * db::add_tag_prefix<LinearSolver::Tags::Operand, typename
25  * Metavariables::system::fields_tag>>.
26  *
27  * Note that the operand \f$p\f$ for which \f$A(p)\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$r = 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(p)\f$ has been computed and stored in the DataBox, the conjugate
36  * gradient algorithm implemented here will converge the field \f$x\f$ towards
37  * the solution and update the operand \f$p\f$ in the process. This requires
38  * two reductions over all elements that are received by a `ResidualMonitor`
39  * singleton parallel component, processed, and then broadcast back to all
40  * elements. The actions are implemented in the `cg_detail` namespace and
41  * constitute the full algorithm in the following order:
42  * 1. `PerformStep` (on elements): Compute the inner product \f$\langle p,
43  * A(p)\rangle\f$ and reduce.
44  * 2. `ComputeAlpha` (on `ResidualMonitor`): Compute
45  * \f$\alpha=\frac{r^2}{\langle p, A(p)\rangle}\f$ and broadcast.
46  * 3. `UpdateFieldValues` (on elements): Update \f$x\f$ and \f$r\f$, then
47  * compute the inner product \f$\langle r, r\rangle\f$ and reduce to find the
48  * new \f$r^2\f$.
49  * 4. `UpdateResidual` (on `ResidualMonitor`): Store the new \f$r^2\f$ and
50  * broadcast the ratio of the new and old \f$r^2\f$, as well as a termination
51  * flag if the `LinearSolver::Tags::ConvergenceCriteria` are met.
52  * 5. `UpdateOperand` (on elements): Update \f$p\f$.
53  *
54  * \see Gmres for a linear solver that can invert nonsymmetric operators
55  * \f$A\f$.
56  */
57 template <typename Metavariables>
59  /*!
60  * \brief The parallel components used by the conjugate gradient linear solver
61  *
62  * Uses:
63  * - System:
64  * * `fields_tag`
65  */
66  using component_list = tmpl::list<cg_detail::ResidualMonitor<Metavariables>>;
67 
68  /*!
69  * \brief Initialize the tags used by the conjugate gradient linear solver
70  *
71  * Uses:
72  * - System:
73  * * `fields_tag`
74  * - ConstGlobalCache: nothing
75  *
76  * With:
77  * - `operand_tag` =
78  * `db::add_tag_prefix<LinearSolver::Tags::Operand, fields_tag>`
79  * - `operator_tag` =
80  * `db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, operand_tag>`
81  * - `residual_tag` =
82  * `db::add_tag_prefix<LinearSolver::Tags::Residual, fields_tag>`
83  *
84  * DataBox changes:
85  * - Adds:
86  * * `LinearSolver::Tags::IterationId`
87  * * `Tags::Next<LinearSolver::Tags::IterationId>`
88  * * `residual_tag`
89  * * `LinearSolver::Tags::HasConverged`
90  * - Removes: nothing
91  * - Modifies:
92  * * `operand_tag`
93  *
94  * \note The `operand_tag` must already be present in the DataBox and is set
95  * to its initial value here. It is typically added to the DataBox by the
96  * system, which uses it to compute the `operator_tag` in each step. Also the
97  * `operator_tag` is typically added to the DataBox by the system, but does
98  * not need to be initialized until it is computed for the first time in the
99  * first step of the algorithm.
100  */
101  using tags = cg_detail::InitializeElement<Metavariables>;
102 
103  // Compile-time interface for observers
104  using observed_reduction_data_tags = observers::make_reduction_data_tags<
105  tmpl::list<observe_detail::reduction_data>>;
106 
107  /*!
108  * \brief Perform an iteration of the conjugate gradient linear solver
109  *
110  * Uses:
111  * - System:
112  * * `fields_tag`
113  * - ConstGlobalCache: nothing
114  *
115  * With:
116  * - `operand_tag` =
117  * `db::add_tag_prefix<LinearSolver::Tags::Operand, fields_tag>`
118  * - `residual_tag` =
119  * `db::add_tag_prefix<LinearSolver::Tags::Residual, fields_tag>`
120  *
121  * DataBox changes:
122  * - Adds: nothing
123  * - Removes: nothing
124  * - Modifies:
125  * * `LinearSolver::Tags::IterationId`
126  * * `Tags::Next<LinearSolver::Tags::IterationId>`
127  * * `fields_tag`
128  * * `operand_tag`
129  * * `residual_tag`
130  * * `LinearSolver::Tags::HasConverged`
131  */
132  using perform_step = cg_detail::PerformStep;
133 };
134 
135 } // namespace LinearSolver
Functionality for solving linear systems of equations.
Definition: TerminateIfConverged.hpp:22
A conjugate gradient solver for linear systems of equations where the operator is symmetric...
Definition: ConjugateGradient.hpp:58
cg_detail::InitializeElement< Metavariables > tags
Initialize the tags used by the conjugate gradient linear solver.
Definition: ConjugateGradient.hpp:101
tmpl::list< cg_detail::ResidualMonitor< Metavariables > > component_list
The parallel components used by the conjugate gradient linear solver.
Definition: ConjugateGradient.hpp:66
Wraps the template metaprogramming library used (brigand)
cg_detail::PerformStep perform_step
Perform an iteration of the conjugate gradient linear solver.
Definition: ConjugateGradient.hpp:132