Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <optional>
       7             : #include <tuple>
       8             : #include <utility>
       9             : 
      10             : #include "DataStructures/DataBox/DataBox.hpp"
      11             : #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
      12             : #include "NumericalAlgorithms/Convergence/Reason.hpp"
      13             : #include "NumericalAlgorithms/Convergence/Tags.hpp"
      14             : #include "Parallel/AlgorithmExecution.hpp"
      15             : #include "ParallelAlgorithms/Actions/Goto.hpp"
      16             : #include "Utilities/Gsl.hpp"
      17             : 
      18             : /// \cond
      19             : namespace tuples {
      20             : template <typename... Tags>
      21             : struct TaggedTuple;
      22             : }  // namespace tuples
      23             : namespace Parallel {
      24             : template <typename Metavariables>
      25             : struct GlobalCache;
      26             : }  // namespace Parallel
      27             : /// \endcond
      28             : 
      29             : namespace LinearSolver {
      30             : /// Actions applicable to any parallel linear solver
      31             : namespace Actions {
      32             : 
      33             : /*!
      34             :  * \brief Make the iterative linear solve the identity operation on the source
      35             :  * vector if no iterations were performed at all. Useful for disabling a
      36             :  * preconditioner by setting its number of iterations to zero.
      37             :  *
      38             :  * When the linear solve is skipped, i.e. when it performs no iterations because
      39             :  * its number of iterations is set to zero, this action assumes \f$A=1\f$ so
      40             :  * \f$Ax=b\f$ solves to \f$x=b\f$. This is useful when the solver is used as
      41             :  * preconditioner, because then we can disable preconditioning by just not
      42             :  * iterating the preconditioner, i.e. by setting its number of iterations to
      43             :  * zero in the input file.
      44             :  *
      45             :  * To use this action, insert it into the action list just after iterating the
      46             :  * linear solver, i.e. after its `solve` action list:
      47             :  *
      48             :  * \snippet LinearSolverAlgorithmTestHelpers.hpp make_identity_if_skipped
      49             :  *
      50             :  * This action will set the `LinearSolverType::fields_tag` to the
      51             :  * `LinearSolverType::source_tag` whenever the linear solver has converged with
      52             :  * the reason `Convergence::Reason::NumIterations` or
      53             :  * `Convergence::Reason::MaxIterations` without actually having performed any
      54             :  * iterations.
      55             :  *
      56             :  * To run additional actions after this action has triggered, i.e. when the
      57             :  * linear solver is skipped, place them after this action and follow them by an
      58             :  * `::Actions::Label<ProceedLabel>`, where `ProceedLabel` is a type used for
      59             :  * identification. Pass the `ProceedLabel` as the second template parameter to
      60             :  * this action. Then, the actions between this action and the label will run
      61             :  * only when the linear solver is skipped. This is useful to set DataBox tags
      62             :  * that are usually updated by the linear solver.
      63             :  *
      64             :  * \par Details:
      65             :  * The standard behaviour of most linear solvers (i.e. when _not_ using this
      66             :  * action) is to keep the fields at their initial guess \f$x=x_0\f$ when it
      67             :  * performs no iterations. That behaviour is not handled well when the solver is
      68             :  * used as preconditioner and its parent always makes it start at \f$x_0=0\f$.
      69             :  * This is a reasonable choice to start the preconditioner but it also means
      70             :  * that the preconditioner doesn't reduce to the identity when it performs no
      71             :  * iterations. Using this action means not iterating the preconditioner at all
      72             :  * essentially disables preconditioning, switching the parent solver to an
      73             :  * unpreconditioned solve with some runtime and memory overhead associated with
      74             :  * initializing the preconditioner.
      75             :  */
      76             : template <typename LinearSolverType, typename ProceedLabel = void>
      77           1 : struct MakeIdentityIfSkipped {
      78             :   template <typename DbTagsList, typename... InboxTags, typename Metavariables,
      79             :             typename ArrayIndex, typename ActionList,
      80             :             typename ParallelComponent>
      81           0 :   static Parallel::iterable_action_return_t apply(
      82             :       db::DataBox<DbTagsList>& box,
      83             :       const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
      84             :       const Parallel::GlobalCache<Metavariables>& /*cache*/,
      85             :       const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
      86             :       const ParallelComponent* const /*meta*/) {
      87             :     constexpr size_t this_action_index =
      88             :         tmpl::index_of<ActionList, MakeIdentityIfSkipped>::value;
      89             :     constexpr size_t proceed_action_index =
      90             :         tmpl::conditional_t<
      91             :             std::is_same_v<ProceedLabel, void>, tmpl::size_t<this_action_index>,
      92             :             tmpl::index_of<ActionList, ::Actions::Label<ProceedLabel>>>::value +
      93             :         1;
      94             : 
      95             :     const auto& has_converged = get<Convergence::Tags::HasConverged<
      96             :         typename LinearSolverType::options_group>>(box);
      97             :     if (has_converged and
      98             :         (has_converged.reason() == Convergence::Reason::NumIterations or
      99             :          has_converged.reason() == Convergence::Reason::MaxIterations) and
     100             :         has_converged.num_iterations() == 0) {
     101             :       db::mutate<typename LinearSolverType::fields_tag>(
     102             :           [](const auto fields, const auto& source) { *fields = source; },
     103             :           make_not_null(&box), get<typename LinearSolverType::source_tag>(box));
     104             :       return {Parallel::AlgorithmExecution::Continue, std::nullopt};
     105             :     }
     106             :     return {Parallel::AlgorithmExecution::Continue, proceed_action_index};
     107             :   }
     108             : };
     109             : 
     110             : namespace detail {
     111             : template <typename Tag>
     112             : struct ProceedLabel {};
     113             : }  // namespace detail
     114             : 
     115             : /// Run `MakeIdentityIfSkipped`, and also run the `BuildOperatorActions` if
     116             : /// the linear solver is skipped. See `MakeIdentityIfSkipped` for details.
     117             : template <typename LinearSolverType,
     118             :           typename BuildOperatorActions = tmpl::list<>,
     119             :           typename Label = typename LinearSolverType::options_group>
     120           1 : using make_identity_if_skipped = tmpl::conditional_t<
     121             :     std::is_same_v<BuildOperatorActions, tmpl::list<>>,
     122             :     MakeIdentityIfSkipped<LinearSolverType>,
     123             :     tmpl::list<
     124             :         MakeIdentityIfSkipped<LinearSolverType, detail::ProceedLabel<Label>>,
     125             :         BuildOperatorActions, ::Actions::Label<detail::ProceedLabel<Label>>>>;
     126             : 
     127             : }  // namespace Actions
     128             : }  // namespace LinearSolver