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