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