Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <blaze/math/DynamicMatrix.h>
7 : #include <blaze/math/DynamicVector.h>
8 : #include <cstddef>
9 : #include <tuple>
10 : #include <utility>
11 :
12 : #include "DataStructures/DataBox/DataBox.hpp"
13 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
14 : #include "IO/Logging/Tags.hpp"
15 : #include "IO/Logging/Verbosity.hpp"
16 : #include "NumericalAlgorithms/Convergence/Tags.hpp"
17 : #include "Parallel/GlobalCache.hpp"
18 : #include "Parallel/Invoke.hpp"
19 : #include "Parallel/Printf/Printf.hpp"
20 : #include "ParallelAlgorithms/LinearSolver/Gmres/Tags/InboxTags.hpp"
21 : #include "ParallelAlgorithms/LinearSolver/Observe.hpp"
22 : #include "ParallelAlgorithms/LinearSolver/Tags.hpp"
23 : #include "Utilities/EqualWithinRoundoff.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::gmres::detail {
36 :
37 : template <typename FieldsTag, typename OptionsGroup, typename BroadcastTarget>
38 : struct InitializeResidualMagnitude {
39 : private:
40 : using fields_tag = FieldsTag;
41 : using residual_magnitude_tag = LinearSolver::Tags::Magnitude<
42 : db::add_tag_prefix<LinearSolver::Tags::Residual, fields_tag>>;
43 : using initial_residual_magnitude_tag =
44 : ::Tags::Initial<residual_magnitude_tag>;
45 : using previous_residual_magnitude_tag =
46 : ::Tags::Previous<residual_magnitude_tag>;
47 : using orthogonalization_history_tag =
48 : LinearSolver::Tags::OrthogonalizationHistory<fields_tag>;
49 :
50 : public:
51 : template <typename ParallelComponent, typename DbTagsList,
52 : typename Metavariables, typename ArrayIndex,
53 : typename DataBox = db::DataBox<DbTagsList>>
54 : static void apply(db::DataBox<DbTagsList>& box,
55 : Parallel::GlobalCache<Metavariables>& cache,
56 : const ArrayIndex& /*array_index*/,
57 : const double residual_magnitude) {
58 : constexpr size_t iteration_id = 0;
59 :
60 : db::mutate<initial_residual_magnitude_tag, previous_residual_magnitude_tag>(
61 : [residual_magnitude](
62 : const gsl::not_null<double*> initial_residual_magnitude,
63 : const gsl::not_null<double*> previous_residual_magnitude) {
64 : *initial_residual_magnitude = residual_magnitude;
65 : *previous_residual_magnitude = residual_magnitude;
66 : },
67 : make_not_null(&box));
68 :
69 : LinearSolver::observe_detail::contribute_to_reduction_observer<
70 : OptionsGroup, ParallelComponent>(iteration_id, residual_magnitude,
71 : cache);
72 :
73 : // Determine whether the linear solver has already converged
74 : Convergence::HasConverged has_converged{
75 : get<Convergence::Tags::Criteria<OptionsGroup>>(box), iteration_id,
76 : residual_magnitude, residual_magnitude};
77 :
78 : // Do some logging
79 : if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(cache) >=
80 : ::Verbosity::Quiet)) {
81 : Parallel::printf("%s initialized with residual: %e\n",
82 : pretty_type::name<OptionsGroup>(), residual_magnitude);
83 : }
84 : if (UNLIKELY(has_converged and get<logging::Tags::Verbosity<OptionsGroup>>(
85 : cache) >= ::Verbosity::Quiet)) {
86 : Parallel::printf("%s has converged without any iterations: %s\n",
87 : pretty_type::name<OptionsGroup>(), has_converged);
88 : }
89 :
90 : Parallel::receive_data<Tags::InitialOrthogonalization<OptionsGroup>>(
91 : Parallel::get_parallel_component<BroadcastTarget>(cache), iteration_id,
92 : // NOLINTNEXTLINE(performance-move-const-arg)
93 : std::make_tuple(residual_magnitude, std::move(has_converged)));
94 : }
95 : };
96 :
97 : template <typename FieldsTag, typename OptionsGroup, typename BroadcastTarget>
98 : struct StoreOrthogonalization {
99 : private:
100 : using fields_tag = FieldsTag;
101 : using residual_magnitude_tag = LinearSolver::Tags::Magnitude<
102 : db::add_tag_prefix<LinearSolver::Tags::Residual, fields_tag>>;
103 : using initial_residual_magnitude_tag =
104 : ::Tags::Initial<residual_magnitude_tag>;
105 : using previous_residual_magnitude_tag =
106 : ::Tags::Previous<residual_magnitude_tag>;
107 : using orthogonalization_history_tag =
108 : LinearSolver::Tags::OrthogonalizationHistory<fields_tag>;
109 :
110 : public:
111 : template <typename ParallelComponent, typename DbTagsList,
112 : typename Metavariables, typename ArrayIndex,
113 : typename DataBox = db::DataBox<DbTagsList>>
114 : static void apply(db::DataBox<DbTagsList>& box,
115 : Parallel::GlobalCache<Metavariables>& cache,
116 : const ArrayIndex& /*array_index*/,
117 : const size_t iteration_id,
118 : const size_t orthogonalization_iteration_id,
119 : const double orthogonalization) {
120 : if (UNLIKELY(orthogonalization_iteration_id == 0)) {
121 : // Append a row and a column to the orthogonalization history. Zero the
122 : // entries that won't be set during the orthogonalization procedure below.
123 : db::mutate<orthogonalization_history_tag>(
124 : [iteration_id](const auto orthogonalization_history) {
125 : orthogonalization_history->resize(iteration_id + 1, iteration_id);
126 : for (size_t j = 0; j < orthogonalization_history->columns() - 1;
127 : ++j) {
128 : (*orthogonalization_history)(
129 : orthogonalization_history->rows() - 1, j) = 0.;
130 : }
131 : },
132 : make_not_null(&box));
133 : }
134 :
135 : // While the orthogonalization procedure is not complete, store the
136 : // orthogonalization, broadcast it back to all elements and return early
137 : if (orthogonalization_iteration_id < iteration_id) {
138 : db::mutate<orthogonalization_history_tag>(
139 : [orthogonalization, iteration_id, orthogonalization_iteration_id](
140 : const auto orthogonalization_history) {
141 : (*orthogonalization_history)(orthogonalization_iteration_id,
142 : iteration_id - 1) = orthogonalization;
143 : },
144 : make_not_null(&box));
145 :
146 : Parallel::receive_data<Tags::Orthogonalization<OptionsGroup>>(
147 : Parallel::get_parallel_component<BroadcastTarget>(cache),
148 : iteration_id, orthogonalization);
149 : return;
150 : }
151 :
152 : // At this point, the orthogonalization procedure is complete.
153 : db::mutate<orthogonalization_history_tag>(
154 : [orthogonalization, iteration_id,
155 : orthogonalization_iteration_id](const auto orthogonalization_history) {
156 : (*orthogonalization_history)(orthogonalization_iteration_id,
157 : iteration_id - 1) =
158 : sqrt(orthogonalization);
159 : },
160 : make_not_null(&box));
161 :
162 : // Perform a QR decomposition of the Hessenberg matrix that was built during
163 : // the orthogonalization
164 : const auto& orthogonalization_history =
165 : get<orthogonalization_history_tag>(box);
166 : const auto num_rows = orthogonalization_iteration_id + 1;
167 : blaze::DynamicMatrix<double> qr_Q;
168 : blaze::DynamicMatrix<double> qr_R;
169 : blaze::qr(orthogonalization_history, qr_Q, qr_R);
170 : // Compute the residual vector from the QR decomposition
171 : blaze::DynamicVector<double> beta(num_rows, 0.);
172 : const double initial_residual_magnitude =
173 : get<initial_residual_magnitude_tag>(box);
174 : beta[0] = initial_residual_magnitude;
175 : blaze::DynamicVector<double> minres =
176 : blaze::inv(qr_R) * blaze::trans(qr_Q) * beta;
177 : const double residual_magnitude =
178 : blaze::length(beta - orthogonalization_history * minres);
179 :
180 : // At this point, the iteration is complete. We proceed with observing,
181 : // logging and checking convergence before broadcasting back to the
182 : // elements.
183 :
184 : LinearSolver::observe_detail::contribute_to_reduction_observer<
185 : OptionsGroup, ParallelComponent>(iteration_id, residual_magnitude,
186 : cache);
187 :
188 : // Determine whether the linear solver has converged.
189 : // GMRES is guaranteed to decrease the residual monotonically, so an
190 : // increase in the residual is an error.
191 : const auto& convergence_criteria =
192 : get<Convergence::Tags::Criteria<OptionsGroup>>(box);
193 : const double previous_residual_magnitude =
194 : get<previous_residual_magnitude_tag>(box);
195 : auto has_converged =
196 : residual_magnitude < previous_residual_magnitude
197 : ? Convergence::HasConverged{convergence_criteria, iteration_id,
198 : residual_magnitude,
199 : initial_residual_magnitude}
200 : : Convergence::HasConverged{
201 : Convergence::Reason::Error,
202 : MakeString{} << std::scientific
203 : << "Residual should decrease monotonically, but "
204 : "increased from "
205 : << previous_residual_magnitude << " to "
206 : << residual_magnitude << ".",
207 : iteration_id};
208 :
209 : db::mutate<previous_residual_magnitude_tag>(
210 : [residual_magnitude](
211 : const gsl::not_null<double*> stored_previous_residual_magnitude) {
212 : *stored_previous_residual_magnitude = residual_magnitude;
213 : },
214 : make_not_null(&box));
215 :
216 : // Do some logging
217 : if (UNLIKELY(get<logging::Tags::Verbosity<OptionsGroup>>(cache) >=
218 : ::Verbosity::Quiet)) {
219 : Parallel::printf("%s(%zu) iteration complete. Remaining residual: %e\n",
220 : pretty_type::name<OptionsGroup>(), iteration_id,
221 : residual_magnitude);
222 : }
223 : if (UNLIKELY(has_converged and get<logging::Tags::Verbosity<OptionsGroup>>(
224 : cache) >= ::Verbosity::Quiet)) {
225 : if (has_converged.reason() == Convergence::Reason::Error) {
226 : Parallel::printf("%s has encountered an error in iteration %zu: %s\n",
227 : pretty_type::name<OptionsGroup>(), iteration_id,
228 : has_converged.error_message());
229 : } else {
230 : Parallel::printf("%s has converged in %zu iterations: %s\n",
231 : pretty_type::name<OptionsGroup>(), iteration_id,
232 : has_converged);
233 : }
234 : }
235 :
236 : Parallel::receive_data<Tags::FinalOrthogonalization<OptionsGroup>>(
237 : Parallel::get_parallel_component<BroadcastTarget>(cache), iteration_id,
238 : std::make_tuple(sqrt(orthogonalization), std::move(minres),
239 : // NOLINTNEXTLINE(performance-move-const-arg)
240 : std::move(has_converged)));
241 : }
242 : };
243 :
244 : } // namespace LinearSolver::gmres::detail
|
