2 // See LICENSE.txt for details.
3
4 #pragma once
5
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
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...