11 #include "DataStructures/DataBox/PrefixHelpers.hpp"
12 #include "Informer/Tags.hpp"
13 #include "Informer/Verbosity.hpp"
14 #include "NumericalAlgorithms/Convergence/HasConverged.hpp"
15 #include "NumericalAlgorithms/Convergence/Tags.hpp"
18 #include "Parallel/Invoke.hpp"
20 #include "ParallelAlgorithms/LinearSolver/ConjugateGradient/Tags/InboxTags.hpp"
21 #include "ParallelAlgorithms/LinearSolver/Observe.hpp"
23 #include "Utilities/EqualWithinRoundoff.hpp"
24 #include "Utilities/Functional.hpp"
30 template <
typename...>
35 namespace LinearSolver::cg::detail {
37 template <
typename FieldsTag,
typename OptionsGroup,
typename BroadcastTarget>
38 struct InitializeResidual {
40 using fields_tag = FieldsTag;
43 using initial_residual_magnitude_tag =
48 template <
typename ParallelComponent,
typename DbTagsList,
49 typename Metavariables,
typename ArrayIndex,
50 typename DataBox = db::DataBox<DbTagsList>,
53 static void apply(db::DataBox<DbTagsList>& box,
56 const double residual_square) noexcept {
57 constexpr
size_t iteration_id = 0;
58 const double residual_magnitude =
sqrt(residual_square);
60 db::mutate<residual_square_tag, initial_residual_magnitude_tag>(
62 [residual_square, residual_magnitude](
65 *local_residual_square = residual_square;
66 *initial_residual_magnitude = residual_magnitude;
69 LinearSolver::observe_detail::contribute_to_reduction_observer<
70 OptionsGroup>(iteration_id, residual_magnitude, cache);
74 get<Convergence::Tags::Criteria<OptionsGroup>>(box), iteration_id,
75 residual_magnitude, residual_magnitude};
79 ::Verbosity::Verbose)) {
81 Options::name<OptionsGroup>() +
82 "' initialized with residual: %e\n",
86 cache) >= ::Verbosity::Quiet)) {
88 Options::name<OptionsGroup>() +
89 "' has converged without any iterations: %s\n",
93 Parallel::receive_data<Tags::InitialHasConverged<OptionsGroup>>(
94 Parallel::get_parallel_component<BroadcastTarget>(cache), iteration_id,
96 std::move(has_converged));
100 template <
typename FieldsTag,
typename OptionsGroup,
typename BroadcastTarget>
101 struct ComputeAlpha {
103 using fields_tag = FieldsTag;
108 template <
typename ParallelComponent,
typename DbTagsList,
109 typename Metavariables,
typename ArrayIndex,
110 typename DataBox = db::DataBox<DbTagsList>,
113 static void apply(db::DataBox<DbTagsList>& box,
116 const size_t iteration_id,
117 const double conj_grad_inner_product) noexcept {
118 Parallel::receive_data<Tags::Alpha<OptionsGroup>>(
119 Parallel::get_parallel_component<BroadcastTarget>(cache), iteration_id,
120 get<residual_square_tag>(box) / conj_grad_inner_product);
124 template <
typename FieldsTag,
typename OptionsGroup,
typename BroadcastTarget>
125 struct UpdateResidual {
127 using fields_tag = FieldsTag;
130 using initial_residual_magnitude_tag =
135 template <
typename ParallelComponent,
typename DbTagsList,
136 typename Metavariables,
typename ArrayIndex,
137 typename DataBox = db::DataBox<DbTagsList>,
140 static void apply(db::DataBox<DbTagsList>& box,
143 const size_t iteration_id,
144 const double residual_square) noexcept {
146 const double res_ratio = residual_square / get<residual_square_tag>(box);
148 db::mutate<residual_square_tag>(
152 *local_residual_square = residual_square;
159 const size_t completed_iterations = iteration_id + 1;
160 const double residual_magnitude =
sqrt(residual_square);
161 LinearSolver::observe_detail::contribute_to_reduction_observer<
162 OptionsGroup>(completed_iterations, residual_magnitude, cache);
166 get<Convergence::Tags::Criteria<OptionsGroup>>(box),
167 completed_iterations, residual_magnitude,
168 get<initial_residual_magnitude_tag>(box)};
172 ::Verbosity::Verbose)) {
174 Options::name<OptionsGroup>() +
175 "' iteration %zu done. Remaining residual: %e\n",
176 completed_iterations, residual_magnitude);
179 cache) >= ::Verbosity::Quiet)) {
181 Options::name<OptionsGroup>() +
182 "' has converged in %zu iterations: %s\n",
183 completed_iterations, has_converged);
186 Parallel::receive_data<Tags::ResidualRatioAndHasConverged<OptionsGroup>>(
187 Parallel::get_parallel_component<BroadcastTarget>(cache), iteration_id,
189 std::make_tuple(res_ratio, std::move(has_converged)));