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 register_element = 126 : async_solvers::RegisterElement<FieldsTag, OptionsGroup, SourceTag, 127 : ArraySectionIdTag>; 128 : template <typename ApplyOperatorActions, typename Label = OptionsGroup> 129 0 : using solve = 130 : tmpl::list<async_solvers::PrepareSolve<FieldsTag, OptionsGroup, SourceTag, 131 : Label, ArraySectionIdTag>, 132 : detail::UpdateFields<FieldsTag, OptionsGroup, SourceTag>, 133 : ApplyOperatorActions, 134 : async_solvers::CompleteStep<FieldsTag, OptionsGroup, SourceTag, 135 : Label, ArraySectionIdTag>>; 136 : }; 137 : } // namespace LinearSolver::Richardson