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/ConjugateGradient/ElementActions.hpp"
       9              : #include "ParallelAlgorithms/LinearSolver/ConjugateGradient/InitializeElement.hpp"
      10              : #include "ParallelAlgorithms/LinearSolver/ConjugateGradient/ResidualMonitor.hpp"
      11              : #include "Utilities/TMPL.hpp"
      12              : 
      13              : /// Items related to the conjugate gradient linear solver
      14              : ///
      15              : /// \see `LinearSolver::cg::ConjugateGradient`
      16            1 : namespace LinearSolver::cg {
      17              : 
      18              : /*!
      19              :  * \ingroup LinearSolverGroup
      20              :  * \brief A conjugate gradient solver for linear systems of equations \f$Ax=b\f$
      21              :  * where the operator \f$A\f$ is symmetric.
      22              :  *
      23              :  * \details The only operation we need to supply to the algorithm is the
      24              :  * result of the operation \f$A(p)\f$ (see \ref LinearSolverGroup) that in the
      25              :  * case of the conjugate gradient algorithm must be symmetric. Each step of the
      26              :  * algorithm expects that \f$A(p)\f$ is computed and stored in the DataBox as
      27              :  * %db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo,
      28              :  * db::add_tag_prefix<LinearSolver::Tags::Operand, FieldsTag>>. To perform a
      29              :  * solve, add the `solve` action list to an array parallel component. Pass the
      30              :  * actions that compute \f$A(q)\f$, as well as any further actions you wish to
      31              :  * run in each step of the algorithm, as the first template parameter to
      32              :  * `solve`. If you add the `solve` action list multiple times, use the second
      33              :  * template parameter to label each solve with a different type.
      34              :  *
      35              :  * Note that the operand \f$p\f$ for which \f$A(p)\f$ needs to be computed is
      36              :  * not the field \f$x\f$ we are solving for but
      37              :  * `db::add_tag_prefix<LinearSolver::Tags::Operand, FieldsTag>`. This field is
      38              :  * initially set to the residual \f$r = b - A(x_0)\f$ where \f$x_0\f$ is the
      39              :  * initial value of the `FieldsTag`.
      40              :  *
      41              :  * When the algorithm step is performed after the operator action \f$A(p)\f$ has
      42              :  * been computed and stored in the DataBox, the conjugate gradient algorithm
      43              :  * implemented here will converge the field \f$x\f$ towards the solution and
      44              :  * update the operand \f$p\f$ in the process. This requires two reductions over
      45              :  * all elements that are received by a `ResidualMonitor` singleton parallel
      46              :  * component, processed, and then broadcast back to all elements. The actions
      47              :  * are implemented in the `cg::detail` namespace and constitute the full
      48              :  * algorithm in the following order:
      49              :  * 1. `PerformStep` (on elements): Compute the inner product \f$\langle p,
      50              :  * A(p)\rangle\f$ and reduce.
      51              :  * 2. `ComputeAlpha` (on `ResidualMonitor`): Compute
      52              :  * \f$\alpha=\frac{r^2}{\langle p, A(p)\rangle}\f$ and broadcast.
      53              :  * 3. `UpdateFieldValues` (on elements): Update \f$x\f$ and \f$r\f$, then
      54              :  * compute the inner product \f$\langle r, r\rangle\f$ and reduce to find the
      55              :  * new \f$r^2\f$.
      56              :  * 4. `UpdateResidual` (on `ResidualMonitor`): Store the new \f$r^2\f$ and
      57              :  * broadcast the ratio of the new and old \f$r^2\f$, as well as a termination
      58              :  * flag if the `Convergence::Tags::Criteria` are met.
      59              :  * 5. `UpdateOperand` (on elements): Update \f$p\f$.
      60              :  *
      61              :  * \see Gmres for a linear solver that can invert nonsymmetric operators
      62              :  * \f$A\f$.
      63              :  */
      64              : template <typename Metavariables, typename FieldsTag, typename OptionsGroup,
      65              :           typename SourceTag =
      66              :               db::add_tag_prefix<::Tags::FixedSource, FieldsTag>>
      67            1 : struct ConjugateGradient {
      68            0 :   using fields_tag = FieldsTag;
      69            0 :   using options_group = OptionsGroup;
      70            0 :   using source_tag = SourceTag;
      71              : 
      72              :   /// Apply the linear operator to this tag in each iteration
      73            1 :   using operand_tag =
      74              :       db::add_tag_prefix<LinearSolver::Tags::Operand, fields_tag>;
      75              : 
      76              :   /*!
      77              :    * \brief The parallel components used by the conjugate gradient linear solver
      78              :    */
      79            1 :   using component_list = tmpl::list<
      80              :       detail::ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>;
      81              : 
      82            0 :   using initialize_element = detail::InitializeElement<FieldsTag, OptionsGroup>;
      83              : 
      84            0 :   using register_element = tmpl::list<>;
      85              : 
      86              :   template <typename ApplyOperatorActions, typename Label = OptionsGroup>
      87            0 :   using solve = tmpl::list<
      88              :       detail::PrepareSolve<FieldsTag, OptionsGroup, Label, SourceTag>,
      89              :       detail::InitializeHasConverged<FieldsTag, OptionsGroup, Label>,
      90              :       ApplyOperatorActions, detail::PerformStep<FieldsTag, OptionsGroup, Label>,
      91              :       detail::UpdateFieldValues<FieldsTag, OptionsGroup, Label>,
      92              :       detail::UpdateOperand<FieldsTag, OptionsGroup, Label>>;
      93              : };
      94              : 
      95              : }  // namespace LinearSolver::cg