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