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
|
