11 #include "DataStructures/DataBox/PrefixHelpers.hpp"
12 #include "NumericalAlgorithms/Convergence/Tags.hpp"
15 #include "Parallel/Invoke.hpp"
16 #include "Parallel/Reduction.hpp"
17 #include "ParallelAlgorithms/LinearSolver/ConjugateGradient/ResidualMonitorActions.hpp"
18 #include "ParallelAlgorithms/LinearSolver/ConjugateGradient/Tags/InboxTags.hpp"
20 #include "Utilities/Functional.hpp"
30 template <
typename...>
33 namespace LinearSolver::cg::detail {
34 template <
typename Metavariables,
typename FieldsTag,
typename OptionsGroup>
35 struct ResidualMonitor;
36 template <
typename FieldsTag,
typename OptionsGroup,
typename Label>
41 namespace LinearSolver::cg::detail {
43 template <
typename FieldsTag,
typename OptionsGroup,
typename Label,
47 using fields_tag = FieldsTag;
48 using source_tag = SourceTag;
49 using operator_applied_to_fields_tag =
57 template <
typename DbTagsList,
typename... InboxTags,
typename Metavariables,
58 typename ArrayIndex,
typename ActionList,
59 typename ParallelComponent>
61 db::DataBox<DbTagsList>& box,
64 const ArrayIndex& array_index,
const ActionList ,
65 const ParallelComponent*
const ) noexcept {
66 db::mutate<Convergence::Tags::IterationId<OptionsGroup>, operand_tag,
70 const auto residual,
const auto& source,
71 const auto& operator_applied_to_fields) noexcept {
73 *operand = source - operator_applied_to_fields;
76 get<source_tag>(box), get<operator_applied_to_fields_tag>(box));
80 const auto& residual = get<residual_tag>(box);
82 InitializeResidual<FieldsTag, OptionsGroup, ParallelComponent>>(
83 Parallel::ReductionData<
86 Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
88 ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache));
90 return {std::move(box)};
94 template <
typename FieldsTag,
typename OptionsGroup,
typename Label>
95 struct InitializeHasConverged {
96 using inbox_tags = tmpl::list<Tags::InitialHasConverged<OptionsGroup>>;
98 template <
typename DbTags,
typename... InboxTags,
typename Metavariables,
100 static bool is_ready(
const db::DataBox<DbTags>& box,
103 const ArrayIndex& ) noexcept {
104 const auto& inbox = get<Tags::InitialHasConverged<OptionsGroup>>(inboxes);
106 box)) != inbox.end();
109 template <
typename DbTagsList,
typename... InboxTags,
typename Metavariables,
110 typename ArrayIndex,
typename ActionList,
111 typename ParallelComponent>
115 const ArrayIndex& ,
const ActionList ,
116 const ParallelComponent*
const ) noexcept {
117 auto has_converged = std::move(
118 tuples::get<Tags::InitialHasConverged<OptionsGroup>>(inboxes)
122 db::mutate<Convergence::Tags::HasConverged<OptionsGroup>>(
125 local_has_converged) noexcept {
126 *local_has_converged = std::move(has_converged);
130 constexpr
size_t step_end_index =
131 tmpl::index_of<ActionList,
132 UpdateOperand<FieldsTag, OptionsGroup, Label>>::value;
133 constexpr
size_t this_action_index =
134 tmpl::index_of<ActionList, InitializeHasConverged>::value;
135 return {std::move(box),
false,
136 has_converged ? (step_end_index + 1) : (this_action_index + 1)};
140 template <
typename FieldsTag,
typename OptionsGroup,
typename Label>
142 template <
typename DbTagsList,
typename... InboxTags,
typename Metavariables,
143 typename ArrayIndex,
typename ActionList,
144 typename ParallelComponent>
146 db::DataBox<DbTagsList>& box,
149 const ArrayIndex& array_index,
const ActionList ,
150 const ParallelComponent*
const ) noexcept {
151 using fields_tag = FieldsTag;
160 const double local_conj_grad_inner_product =
161 inner_product(get<operand_tag>(box), get<operator_tag>(box));
164 ComputeAlpha<FieldsTag, OptionsGroup, ParallelComponent>>(
165 Parallel::ReductionData<
168 get<Convergence::Tags::IterationId<OptionsGroup>>(box),
169 local_conj_grad_inner_product},
170 Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
172 ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache));
174 return {std::move(box)};
178 template <
typename FieldsTag,
typename OptionsGroup,
typename Label>
179 struct UpdateFieldValues {
181 using fields_tag = FieldsTag;
190 using inbox_tags = tmpl::list<Tags::Alpha<OptionsGroup>>;
192 template <
typename DbTags,
typename... InboxTags,
typename Metavariables,
194 static bool is_ready(
const db::DataBox<DbTags>& box,
197 const ArrayIndex& ) noexcept {
198 const auto& inbox = get<Tags::Alpha<OptionsGroup>>(inboxes);
200 box)) != inbox.end();
203 template <
typename DbTagsList,
typename... InboxTags,
typename Metavariables,
204 typename ArrayIndex,
typename ActionList,
205 typename ParallelComponent>
209 const ArrayIndex& array_index,
const ActionList ,
210 const ParallelComponent*
const ) noexcept {
211 const double alpha = std::move(
216 db::mutate<residual_tag, fields_tag>(
218 [alpha](
const auto residual,
const auto fields,
const auto& operand,
219 const auto& operator_applied_to_operand) noexcept {
220 *fields += alpha * operand;
221 *residual -= alpha * operator_applied_to_operand;
223 get<operand_tag>(box), get<operator_tag>(box));
226 const auto& residual = get<residual_tag>(box);
227 const double local_residual_magnitude_square =
231 UpdateResidual<FieldsTag, OptionsGroup, ParallelComponent>>(
232 Parallel::ReductionData<
235 get<Convergence::Tags::IterationId<OptionsGroup>>(box),
236 local_residual_magnitude_square},
237 Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
239 ResidualMonitor<Metavariables, FieldsTag, OptionsGroup>>(cache));
241 return {std::move(box)};
245 template <
typename FieldsTag,
typename OptionsGroup,
typename Label>
246 struct UpdateOperand {
248 using fields_tag = FieldsTag;
256 tmpl::list<Tags::ResidualRatioAndHasConverged<OptionsGroup>>;
258 template <
typename DbTags,
typename... InboxTags,
typename Metavariables,
260 static bool is_ready(
const db::DataBox<DbTags>& box,
263 const ArrayIndex& ) noexcept {
265 get<Tags::ResidualRatioAndHasConverged<OptionsGroup>>(inboxes);
267 box)) != inbox.end();
270 template <
typename DbTagsList,
typename... InboxTags,
typename Metavariables,
271 typename ArrayIndex,
typename ActionList,
272 typename ParallelComponent>
276 const ArrayIndex& ,
const ActionList ,
277 const ParallelComponent*
const ) noexcept {
278 auto received_data = std::move(
279 tuples::get<Tags::ResidualRatioAndHasConverged<OptionsGroup>>(inboxes)
282 const double res_ratio = get<0>(received_data);
283 auto& has_converged = get<1>(received_data);
285 db::mutate<operand_tag, Convergence::Tags::IterationId<OptionsGroup>,
288 [res_ratio, &has_converged](
291 const auto& residual) noexcept {
292 *operand = residual + res_ratio * *operand;
294 *local_has_converged = std::move(has_converged);
296 get<residual_tag>(box));
299 constexpr
size_t this_action_index =
300 tmpl::index_of<ActionList, UpdateOperand>::value;
301 constexpr
size_t prepare_step_index =
302 tmpl::index_of<ActionList, InitializeHasConverged<
303 FieldsTag, OptionsGroup, Label>>::value +
305 return {std::move(box),
false,
306 get<Convergence::Tags::HasConverged<OptionsGroup>>(box)
307 ? (this_action_index + 1)
308 : prepare_step_index};