MakeIdentityIfSkipped.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <tuple>
7 #include <utility>
8 
10 #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
11 #include "NumericalAlgorithms/Convergence/Reason.hpp"
12 #include "NumericalAlgorithms/Convergence/Tags.hpp"
13 #include "Utilities/Gsl.hpp"
14 
15 /// \cond
16 namespace tuples {
17 template <typename... Tags>
18 struct TaggedTuple;
19 } // namespace tuples
20 namespace Parallel {
21 template <typename Metavariables>
22 struct GlobalCache;
23 } // namespace Parallel
24 /// \endcond
25 
26 namespace LinearSolver {
27 /// Actions applicable to any parallel linear solver
28 namespace Actions {
29 
30 /*!
31  * \brief Make the iterative linear solve the identity operation on the source
32  * vector if no iterations were performed at all. Useful for disabling a
33  * preconditioner by setting its number of iterations to zero.
34  *
35  * When the linear solve is skipped, i.e. when it performs no iterations because
36  * its number of iterations is set to zero, this action assumes \f$A=1\f$ so
37  * \f$Ax=b\f$ solves to \f$x=b\f$. This is useful when the solver is used as
38  * preconditioner, because then we can disable preconditioning by just not
39  * iterating the preconditioner, i.e. by setting its number of iterations to
40  * zero in the input file.
41  *
42  * To use this action, insert it into the action list just after iterating the
43  * linear solver, i.e. after its `solve` action list:
44  *
45  * \snippet LinearSolverAlgorithmTestHelpers.hpp make_identity_if_skipped
46  *
47  * This action will set the `LinearSolverType::fields_tag` to the
48  * `LinearSolverType::source_tag` whenever the linear solver has converged with
49  * the reason `Convergence::Reason::NumIterations` or
50  * `Convergence::Reason::MaxIterations` without actually having performed any
51  * iterations.
52  *
53  * \par Details:
54  * The standard behaviour of most linear solvers (i.e. when _not_ using this
55  * action) is to keep the fields at their initial guess \f$x=x_0\f$ when it
56  * performs no iterations. That behaviour is not handled well when the solver is
57  * used as preconditioner and its parent always makes it start at \f$x_0=0\f$.
58  * This is a reasonable choice to start the preconditioner but it also means
59  * that the preconditioner doesn't reduce to the identity when it performs no
60  * iterations. Using this action means not iterating the preconditioner at all
61  * essentially disables preconditioning, switching the parent solver to an
62  * unpreconditioned solve with some runtime and memory overhead associated with
63  * initializing the preconditioner.
64  */
65 template <typename LinearSolverType>
67  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
68  typename ArrayIndex, typename ActionList,
69  typename ParallelComponent>
70  static std::tuple<db::DataBox<DbTagsList>&&> apply(
71  db::DataBox<DbTagsList>& box,
72  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
73  const Parallel::GlobalCache<Metavariables>& /*cache*/,
74  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
75  const ParallelComponent* const /*meta*/) noexcept {
76  const auto& has_converged = get<Convergence::Tags::HasConverged<
77  typename LinearSolverType::options_group>>(box);
78  if (has_converged and
79  (has_converged.reason() == Convergence::Reason::NumIterations or
80  has_converged.reason() == Convergence::Reason::MaxIterations) and
81  has_converged.num_iterations() == 0) {
82  db::mutate<typename LinearSolverType::fields_tag>(
83  make_not_null(&box),
84  [](const auto fields, const auto& source) noexcept {
85  *fields = source;
86  },
87  get<typename LinearSolverType::source_tag>(box));
88  }
89  return {std::move(box)};
90  }
91 };
92 
93 } // namespace Actions
94 } // namespace LinearSolver
utility
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:660
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:15
tuple
LinearSolver::Actions::MakeIdentityIfSkipped
Make the iterative linear solve the identity operation on the source vector if no iterations were per...
Definition: MakeIdentityIfSkipped.hpp:66
DataBox.hpp
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
Convergence::Reason::NumIterations
@ NumIterations
Reached the target number of iterations.
LinearSolver
Functionality for solving linear systems of equations.
Definition: ExplicitInverse.hpp:20
Gsl.hpp
make_not_null
gsl::not_null< T * > make_not_null(T *ptr) noexcept
Construct a not_null from a pointer. Often this will be done as an implicit conversion,...
Definition: Gsl.hpp:880
Convergence::Tags::HasConverged
Holds a Convergence::HasConverged flag that signals the iterative algorithm has converged,...
Definition: Tags.hpp:86
Convergence::Reason::MaxIterations
@ MaxIterations
Reached the maximum number of iterations. Can be interpreted as an error condition.
Parallel
Functionality for parallelization.
Definition: ElementReceiveInterpPoints.hpp:13