Line data    Source code

       1            0 : // Distributed under the MIT License.
       2              : // See LICENSE.txt for details.
       3              : 
       4              : #pragma once
       5              : 
       6              : #include <cstddef>
       7              : #include <optional>
       8              : #include <tuple>
       9              : #include <utility>
      10              : 
      11              : #include "DataStructures/DataBox/DataBox.hpp"
      12              : #include "DataStructures/DataBox/PrefixHelpers.hpp"
      13              : #include "DataStructures/DataBox/Prefixes.hpp"
      14              : #include "IO/Observer/Helpers.hpp"
      15              : #include "Parallel/AlgorithmExecution.hpp"
      16              : #include "ParallelAlgorithms/LinearSolver/AsynchronousSolvers/ElementActions.hpp"
      17              : #include "ParallelAlgorithms/LinearSolver/Richardson/Tags.hpp"
      18              : #include "ParallelAlgorithms/LinearSolver/Tags.hpp"
      19              : #include "Utilities/TMPL.hpp"
      20              : 
      21              : /// \cond
      22              : namespace tuples {
      23              : template <typename...>
      24              : class TaggedTuple;
      25              : }  // namespace tuples
      26              : namespace Parallel {
      27              : template <typename Metavariables>
      28              : struct GlobalCache;
      29              : }  // namespace Parallel
      30              : /// \endcond
      31              : 
      32              : /// Items related to the %Richardson linear solver
      33              : ///
      34              : /// \see `LinearSolver::Richardson::Richardson`
      35            1 : namespace LinearSolver::Richardson {
      36              : 
      37              : namespace detail {
      38              : 
      39              : template <typename FieldsTag, typename OptionsGroup, typename SourceTag>
      40              : struct UpdateFields {
      41              :  private:
      42              :   using residual_tag =
      43              :       db::add_tag_prefix<LinearSolver::Tags::Residual, FieldsTag>;
      44              : 
      45              :  public:
      46              :   using const_global_cache_tags =
      47              :       tmpl::list<Tags::RelaxationParameter<OptionsGroup>>;
      48              : 
      49              :   template <typename DbTagsList, typename... InboxTags, typename Metavariables,
      50              :             typename ArrayIndex, typename ActionList,
      51              :             typename ParallelComponent>
      52              :   static Parallel::iterable_action_return_t apply(
      53              :       db::DataBox<DbTagsList>& box,
      54              :       const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
      55              :       const Parallel::GlobalCache<Metavariables>& /*cache*/,
      56              :       const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
      57              :       const ParallelComponent* const /*meta*/) {
      58              :     // Update the solution fields according to the Richardson scheme
      59              :     db::mutate<FieldsTag>(
      60              :         [](const auto fields, const auto& residual,
      61              :            const double relaxation_parameter) {
      62              :           *fields += relaxation_parameter * residual;
      63              :         },
      64              :         make_not_null(&box), get<residual_tag>(box),
      65              :         get<Tags::RelaxationParameter<OptionsGroup>>(box));
      66              :     return {Parallel::AlgorithmExecution::Continue, std::nullopt};
      67              :   }
      68              : };
      69              : 
      70              : }  // namespace detail
      71              : 
      72              : /*!
      73              :  * \ingroup LinearSolverGroup
      74              :  * \brief A simple %Richardson scheme for solving a system of linear equations
      75              :  * \f$Ax=b\f$
      76              :  *
      77              :  * \warning This linear solver is useful only for basic preconditioning of
      78              :  * another linear solver or for testing purposes. See
      79              :  * `LinearSolver::cg::ConjugateGradient` or `LinearSolver::gmres::Gmres` for
      80              :  * more useful general-purpose linear solvers.
      81              :  *
      82              :  * In each step the solution is updated from its initial state \f$x_0\f$ as
      83              :  *
      84              :  * \f[
      85              :  * x_{k+1} = x_k + \omega \left(b - Ax\right)
      86              :  * \f]
      87              :  *
      88              :  * where \f$\omega\f$ is a _relaxation parameter_ that weights the residual.
      89              :  *
      90              :  * The scheme converges if the spectral radius (i.e. the largest absolute
      91              :  * eigenvalue) of the iteration operator \f$G=1-\omega A\f$ is smaller than one.
      92              :  * For symmetric positive definite (SPD) matrices \f$A\f$ with largest
      93              :  * eigenvalue \f$\lambda_\mathrm{max}\f$ and smallest eigenvalue
      94              :  * \f$\lambda_\mathrm{min}\f$ choose
      95              :  *
      96              :  * \f[
      97              :  * \omega_\mathrm{SPD,optimal} = \frac{2}{\lambda_\mathrm{max} +
      98              :  * \lambda_\mathrm{min}}
      99              :  * \f]
     100              :  *
     101              :  * for optimal convergence.
     102              :  *
     103              :  * \par Array sections
     104              :  * This linear solver requires no synchronization between elements, so it runs
     105              :  * on all elements in the array parallel component. Partitioning of the elements
     106              :  * in sections is only relevant for observing residual norms. Pass the section
     107              :  * ID tag for the `ArraySectionIdTag` template parameter if residual norms
     108              :  * should be computed over a section. Pass `void` (default) to compute residual
     109              :  * norms over all elements in the array.
     110              :  */
     111              : template <typename FieldsTag, typename OptionsGroup,
     112              :           typename SourceTag =
     113              :               db::add_tag_prefix<::Tags::FixedSource, FieldsTag>,
     114              :           typename ArraySectionIdTag = void>
     115            1 : struct Richardson {
     116            0 :   using fields_tag = FieldsTag;
     117            0 :   using options_group = OptionsGroup;
     118            0 :   using source_tag = SourceTag;
     119            0 :   using operand_tag = fields_tag;
     120            0 :   using component_list = tmpl::list<>;
     121            0 :   using observed_reduction_data_tags = observers::make_reduction_data_tags<
     122              :       tmpl::list<async_solvers::reduction_data>>;
     123            0 :   using initialize_element =
     124              :       async_solvers::InitializeElement<FieldsTag, OptionsGroup, SourceTag>;
     125            0 :   using amr_projectors = initialize_element;
     126            0 :   using register_element =
     127              :       async_solvers::RegisterElement<FieldsTag, OptionsGroup, SourceTag,
     128              :                                      ArraySectionIdTag>;
     129              :   template <typename ApplyOperatorActions, typename Label = OptionsGroup>
     130            0 :   using solve =
     131              :       tmpl::list<async_solvers::PrepareSolve<FieldsTag, OptionsGroup, SourceTag,
     132              :                                              Label, ArraySectionIdTag>,
     133              :                  detail::UpdateFields<FieldsTag, OptionsGroup, SourceTag>,
     134              :                  ApplyOperatorActions,
     135              :                  async_solvers::CompleteStep<FieldsTag, OptionsGroup, SourceTag,
     136              :                                              Label, ArraySectionIdTag>>;
     137              : };
     138              : }  // namespace LinearSolver::Richardson