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