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 <tuple>
       8             : #include <utility>
       9             : 
      10             : #include "DataStructures/DataBox/DataBox.hpp"
      11             : #include "DataStructures/DataBox/PrefixHelpers.hpp"
      12             : #include "IO/Logging/Tags.hpp"
      13             : #include "IO/Logging/Verbosity.hpp"
      14             : #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
      15             : #include "NumericalAlgorithms/Convergence/Tags.hpp"
      16             : #include "Parallel/GlobalCache.hpp"
      17             : #include "Parallel/Invoke.hpp"
      18             : #include "Parallel/Printf/Printf.hpp"
      19             : #include "ParallelAlgorithms/LinearSolver/ConjugateGradient/Tags/InboxTags.hpp"
      20             : #include "ParallelAlgorithms/LinearSolver/Observe.hpp"
      21             : #include "ParallelAlgorithms/LinearSolver/Tags.hpp"
      22             : #include "Utilities/EqualWithinRoundoff.hpp"
      23             : #include "Utilities/Functional.hpp"
      24             : #include "Utilities/Gsl.hpp"
      25             : #include "Utilities/PrettyType.hpp"
      26             : #include "Utilities/Requires.hpp"
      27             : 
      28             : /// \cond
      29             : namespace tuples {
      30             : template <typename...>
      31             : class TaggedTuple;
      32             : }  // namespace tuples
      33             : /// \endcond
      34             : 
      35             : namespace LinearSolver::cg::detail {
      36             : 
      37             : template <typename FieldsTag, typename OptionsGroup, typename BroadcastTarget>
      38             : struct InitializeResidual {
      39             :  private:
      40             :   using fields_tag = FieldsTag;
      41             :   using residual_square_tag = LinearSolver::Tags::MagnitudeSquare<
      42             :       db::add_tag_prefix<LinearSolver::Tags::Residual, fields_tag>>;
      43             :   using initial_residual_magnitude_tag =
      44             :       ::Tags::Initial<LinearSolver::Tags::Magnitude<
      45             :           db::add_tag_prefix<LinearSolver::Tags::Residual, fields_tag>>>;
      46             : 
      47             :  public:
      48             :   template <typename ParallelComponent, typename DbTagsList,
      49             :             typename Metavariables, typename ArrayIndex,
      50             :             typename DataBox = db::DataBox<DbTagsList>>
      51             :   static void apply(db::DataBox<DbTagsList>& box,
      52             :                     Parallel::GlobalCache<Metavariables>& cache,
      53             :                     const ArrayIndex& /*array_index*/,
      54             :                     const double residual_square) {
      55             :     constexpr size_t iteration_id = 0;
      56             :     const double residual_magnitude = sqrt(residual_square);
      57             : 
      58             :     db::mutate<residual_square_tag, initial_residual_magnitude_tag>(
      59             :         [residual_square, residual_magnitude](
      60             :             const gsl::not_null<double*> local_residual_square,
      61             :             const gsl::not_null<double*> initial_residual_magnitude) {
      62             :           *local_residual_square = residual_square;
      63             :           *initial_residual_magnitude = residual_magnitude;
      64             :         },
      65             :         make_not_null(&box));
      66             : 
      67             :     LinearSolver::observe_detail::contribute_to_reduction_observer<
      68             :         OptionsGroup, ParallelComponent>(iteration_id, residual_magnitude,
      69             :                                          cache);
      70             : 
      71             :     // Determine whether the linear solver has converged
      72             :     Convergence::HasConverged has_converged{
      73             :         get<Convergence::Tags::Criteria<OptionsGroup>>(box), iteration_id,
      74             :         residual_magnitude, residual_magnitude};
      75             : 
      76             :     // Do some logging
      77             :     if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(cache) >=
      78             :                  ::Verbosity::Quiet)) {
      79             :       Parallel::printf("%s initialized with residual: %e\n",
      80             :                        pretty_type::name<OptionsGroup>(), residual_magnitude);
      81             :     }
      82             :     if (UNLIKELY(has_converged and get<logging::Tags::Verbosity<OptionsGroup>>(
      83             :                                        cache) >= ::Verbosity::Quiet)) {
      84             :       Parallel::printf("%s has converged without any iterations: %s\n",
      85             :                        pretty_type::name<OptionsGroup>(), has_converged);
      86             :     }
      87             : 
      88             :     Parallel::receive_data<Tags::InitialHasConverged<OptionsGroup>>(
      89             :         Parallel::get_parallel_component<BroadcastTarget>(cache), iteration_id,
      90             :         // NOLINTNEXTLINE(performance-move-const-arg)
      91             :         std::move(has_converged));
      92             :   }
      93             : };
      94             : 
      95             : template <typename FieldsTag, typename OptionsGroup, typename BroadcastTarget>
      96             : struct ComputeAlpha {
      97             :  private:
      98             :   using fields_tag = FieldsTag;
      99             :   using residual_square_tag = LinearSolver::Tags::MagnitudeSquare<
     100             :       db::add_tag_prefix<LinearSolver::Tags::Residual, fields_tag>>;
     101             : 
     102             :  public:
     103             :   template <typename ParallelComponent, typename DbTagsList,
     104             :             typename Metavariables, typename ArrayIndex,
     105             :             typename DataBox = db::DataBox<DbTagsList>>
     106             :   static void apply(db::DataBox<DbTagsList>& box,
     107             :                     Parallel::GlobalCache<Metavariables>& cache,
     108             :                     const ArrayIndex& /*array_index*/,
     109             :                     const size_t iteration_id,
     110             :                     const double conj_grad_inner_product) {
     111             :     Parallel::receive_data<Tags::Alpha<OptionsGroup>>(
     112             :         Parallel::get_parallel_component<BroadcastTarget>(cache), iteration_id,
     113             :         get<residual_square_tag>(box) / conj_grad_inner_product);
     114             :   }
     115             : };
     116             : 
     117             : template <typename FieldsTag, typename OptionsGroup, typename BroadcastTarget>
     118             : struct UpdateResidual {
     119             :  private:
     120             :   using fields_tag = FieldsTag;
     121             :   using residual_square_tag = LinearSolver::Tags::MagnitudeSquare<
     122             :       db::add_tag_prefix<LinearSolver::Tags::Residual, fields_tag>>;
     123             :   using initial_residual_magnitude_tag =
     124             :       ::Tags::Initial<LinearSolver::Tags::Magnitude<
     125             :           db::add_tag_prefix<LinearSolver::Tags::Residual, fields_tag>>>;
     126             : 
     127             :  public:
     128             :   template <typename ParallelComponent, typename DbTagsList,
     129             :             typename Metavariables, typename ArrayIndex,
     130             :             typename DataBox = db::DataBox<DbTagsList>>
     131             :   static void apply(db::DataBox<DbTagsList>& box,
     132             :                     Parallel::GlobalCache<Metavariables>& cache,
     133             :                     const ArrayIndex& /*array_index*/,
     134             :                     const size_t iteration_id, const double residual_square) {
     135             :     // Compute the residual ratio before mutating the DataBox
     136             :     const double res_ratio = residual_square / get<residual_square_tag>(box);
     137             : 
     138             :     db::mutate<residual_square_tag>(
     139             :         [residual_square](const gsl::not_null<double*> local_residual_square) {
     140             :           *local_residual_square = residual_square;
     141             :         },
     142             :         make_not_null(&box));
     143             : 
     144             :     // At this point, the iteration is complete. We proceed with observing,
     145             :     // logging and checking convergence before broadcasting back to the
     146             :     // elements.
     147             : 
     148             :     const size_t completed_iterations = iteration_id + 1;
     149             :     const double residual_magnitude = sqrt(residual_square);
     150             :     LinearSolver::observe_detail::contribute_to_reduction_observer<
     151             :         OptionsGroup, ParallelComponent>(completed_iterations,
     152             :                                          residual_magnitude, cache);
     153             : 
     154             :     // Determine whether the linear solver has converged
     155             :     Convergence::HasConverged has_converged{
     156             :         get<Convergence::Tags::Criteria<OptionsGroup>>(box),
     157             :         completed_iterations, residual_magnitude,
     158             :         get<initial_residual_magnitude_tag>(box)};
     159             : 
     160             :     // Do some logging
     161             :     if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(cache) >=
     162             :                  ::Verbosity::Quiet)) {
     163             :       Parallel::printf("%s(%zu) iteration complete. Remaining residual: %e\n",
     164             :                        pretty_type::name<OptionsGroup>(), completed_iterations,
     165             :                        residual_magnitude);
     166             :     }
     167             :     if (UNLIKELY(has_converged and get<logging::Tags::Verbosity<OptionsGroup>>(
     168             :                                        cache) >= ::Verbosity::Quiet)) {
     169             :       Parallel::printf("%s has converged in %zu iterations: %s\n",
     170             :                        pretty_type::name<OptionsGroup>(), completed_iterations,
     171             :                        has_converged);
     172             :     }
     173             : 
     174             :     Parallel::receive_data<Tags::ResidualRatioAndHasConverged<OptionsGroup>>(
     175             :         Parallel::get_parallel_component<BroadcastTarget>(cache), iteration_id,
     176             :         // NOLINTNEXTLINE(performance-move-const-arg)
     177             :         std::make_tuple(res_ratio, std::move(has_converged)));
     178             :   }
     179             : };
     180             : 
     181             : }  // namespace LinearSolver::cg::detail