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