Richardson.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 #include <tuple>
8 #include <utility>
9 
11 #include "DataStructures/DataBox/PrefixHelpers.hpp"
13 #include "IO/Observer/Helpers.hpp"
14 #include "ParallelAlgorithms/LinearSolver/AsynchronousSolvers/ElementActions.hpp"
15 #include "ParallelAlgorithms/LinearSolver/Richardson/Tags.hpp"
17 #include "Utilities/TMPL.hpp"
18 
19 /// \cond
20 namespace tuples {
21 template <typename...>
22 class TaggedTuple;
23 } // namespace tuples
24 namespace Parallel {
25 template <typename Metavariables>
26 struct GlobalCache;
27 } // namespace Parallel
28 /// \endcond
29 
30 /// Items related to the %Richardson linear solver
31 ///
32 /// \see `LinearSolver::Richardson::Richardson`
34 
35 namespace detail {
36 
37 template <typename FieldsTag, typename OptionsGroup, typename SourceTag>
38 struct UpdateFields {
39  private:
40  using residual_tag =
42 
43  public:
44  using const_global_cache_tags =
45  tmpl::list<Tags::RelaxationParameter<OptionsGroup>>;
46 
47  template <typename DbTagsList, typename... InboxTags, typename Metavariables,
48  typename ArrayIndex, typename ActionList,
49  typename ParallelComponent>
50  static std::tuple<db::DataBox<DbTagsList>&&> apply(
51  db::DataBox<DbTagsList>& box,
52  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
53  const Parallel::GlobalCache<Metavariables>& /*cache*/,
54  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
55  const ParallelComponent* const /*meta*/) noexcept {
56  // Update the solution fields according to the Richardson scheme
57  db::mutate<FieldsTag>(
58  make_not_null(&box),
59  [](const auto fields, const auto& residual,
60  const double relaxation_parameter) noexcept {
61  *fields += relaxation_parameter * residual;
62  },
63  get<residual_tag>(box),
64  get<Tags::RelaxationParameter<OptionsGroup>>(box));
65  return {std::move(box)};
66  }
67 };
68 
69 } // namespace detail
70 
71 /*!
72  * \ingroup LinearSolverGroup
73  * \brief A simple %Richardson scheme for solving a system of linear equations
74  * \f$Ax=b\f$
75  *
76  * \warning This linear solver is useful only for basic preconditioning of
77  * another linear solver or for testing purposes. See
78  * `LinearSolver::cg::ConjugateGradient` or `LinearSolver::gmres::Gmres` for
79  * more useful general-purpose linear solvers.
80  *
81  * In each step the solution is updated from its initial state \f$x_0\f$ as
82  *
83  * \f[
84  * x_{k+1} = x_k + \omega \left(b - Ax\right)
85  * \f]
86  *
87  * where \f$\omega\f$ is a _relaxation parameter_ that weights the residual.
88  *
89  * The scheme converges if the spectral radius (i.e. the largest absolute
90  * eigenvalue) of the iteration operator \f$G=1-\omega A\f$ is smaller than one.
91  * For symmetric positive definite (SPD) matrices \f$A\f$ with largest
92  * eigenvalue \f$\lambda_\mathrm{max}\f$ and smallest eigenvalue
93  * \f$\lambda_\mathrm{min}\f$ choose
94  *
95  * \f[
96  * \omega_\mathrm{SPD,optimal} = \frac{2}{\lambda_\mathrm{max} +
97  * \lambda_\mathrm{min}}
98  * \f]
99  *
100  * for optimal convergence.
101  */
102 template <typename FieldsTag, typename OptionsGroup,
103  typename SourceTag =
105 struct Richardson {
106  using fields_tag = FieldsTag;
107  using options_group = OptionsGroup;
108  using source_tag = SourceTag;
109  using operand_tag = fields_tag;
110  using component_list = tmpl::list<>;
111  using observed_reduction_data_tags = observers::make_reduction_data_tags<
112  tmpl::list<async_solvers::reduction_data>>;
113  using initialize_element =
115  using register_element =
117  template <typename ApplyOperatorActions, typename Label = OptionsGroup>
118  using solve = tmpl::list<
120  detail::UpdateFields<FieldsTag, OptionsGroup, SourceTag>,
121  ApplyOperatorActions,
123 };
124 } // namespace LinearSolver::Richardson
LinearSolver::Richardson::Richardson
A simple Richardson scheme for solving a system of linear equations .
Definition: Richardson.hpp:105
utility
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:15
db::add_tag_prefix
typename detail::add_tag_prefix_impl< Prefix, Tag, Args... >::type add_tag_prefix
Definition: PrefixHelpers.hpp:51
tuple
LinearSolver::async_solvers::InitializeElement
Definition: ElementActions.hpp:102
Tags.hpp
DataBox.hpp
cstddef
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
LinearSolver::Richardson
Items related to the Richardson linear solver.
Definition: Richardson.hpp:33
LinearSolver::async_solvers::PrepareSolve
Definition: ElementActions.hpp:160
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
Prefixes.hpp
Parallel
Contains functions that forward to Charm++ parallel functions.
Definition: ElementReceiveInterpPoints.hpp:13
TMPL.hpp
observers::Actions::RegisterWithObservers
Register an observation ID with the observers.
Definition: RegisterWithObservers.hpp:39
LinearSolver::async_solvers::CompleteStep
Definition: ElementActions.hpp:214