Line data Source code
1 1 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : /// \file
5 : /// Functionality to build the explicit matrix representation of the linear
6 : /// operator column-by-column, and to solve the linear problem directly with it.
7 :
8 : #pragma once
9 :
10 : #include <blaze/math/Submatrix.h>
11 : #include <blaze/math/lapack/getrf.h>
12 : #include <blaze/math/lapack/getrs.h>
13 : #include <cstddef>
14 : #include <map>
15 : #include <optional>
16 : #include <utility>
17 : #include <vector>
18 :
19 : #include "DataStructures/CompressedMatrix.hpp"
20 : #include "DataStructures/DataBox/DataBox.hpp"
21 : #include "DataStructures/DataBox/Tag.hpp"
22 : #include "DataStructures/TaggedTuple.hpp"
23 : #include "Domain/Structure/ElementId.hpp"
24 : #include "Domain/Tags.hpp"
25 : #include "IO/H5/TensorData.hpp"
26 : #include "IO/Logging/Tags.hpp"
27 : #include "IO/Logging/Verbosity.hpp"
28 : #include "IO/Observer/Actions/RegisterWithObservers.hpp"
29 : #include "IO/Observer/GetSectionObservationKey.hpp"
30 : #include "IO/Observer/VolumeActions.hpp"
31 : #include "NumericalAlgorithms/Convergence/Tags.hpp"
32 : #include "Parallel/AlgorithmExecution.hpp"
33 : #include "Parallel/GetSection.hpp"
34 : #include "Parallel/Phase.hpp"
35 : #include "Parallel/Printf/Printf.hpp"
36 : #include "Parallel/Reduction.hpp"
37 : #include "Parallel/Section.hpp"
38 : #include "Parallel/Tags/Section.hpp"
39 : #include "Parallel/TypeTraits.hpp"
40 : #include "ParallelAlgorithms/Actions/Goto.hpp"
41 : #include "ParallelAlgorithms/Amr/Protocols/Projector.hpp"
42 : #include "ParallelAlgorithms/LinearSolver/Tags.hpp"
43 : #include "Utilities/EqualWithinRoundoff.hpp"
44 : #include "Utilities/Functional.hpp"
45 : #include "Utilities/GetOutput.hpp"
46 : #include "Utilities/ProtocolHelpers.hpp"
47 : #include "Utilities/TMPL.hpp"
48 :
49 : namespace LinearSolver {
50 0 : namespace OptionTags {
51 :
52 : /// Options for building the explicit matrix representation of the linear
53 : /// operator.
54 1 : struct BuildMatrixOptionsGroup {
55 0 : static std::string name() { return "BuildMatrix"; }
56 0 : static constexpr Options::String help = {
57 : "Options for building the explicit matrix representation of the linear "
58 : "operator. This is done by applying the linear operator to unit "
59 : "vectors. The matrix can be written to disk for analysis, or used to "
60 : "solve small problems directly, e.g. on the coarsest multigrid level."};
61 : };
62 :
63 : /// Subfile name in the volume data H5 files where the matrix will be stored.
64 1 : struct MatrixSubfileName {
65 0 : using type = Options::Auto<std::string, Options::AutoLabel::None>;
66 0 : using group = BuildMatrixOptionsGroup;
67 0 : static constexpr Options::String help = {
68 : "Subfile name in the volume data H5 files where the matrix will be "
69 : "stored. Each observation in the subfile is a column of the matrix. The "
70 : "row index is the order of elements defined by the ElementId in the "
71 : "volume data, by the order of tensor components encoded in the name of "
72 : "the components, and by the contiguous ordering of grid points for each "
73 : "component."};
74 : };
75 :
76 : /// Option for enabling direct solve of the linear problem.
77 1 : struct EnableDirectSolve {
78 0 : using type = bool;
79 0 : using group = BuildMatrixOptionsGroup;
80 0 : static constexpr Options::String help = {
81 : "Solve the linear problem directly with an LU factorization of the "
82 : "assembled matrix. The factorization is reused for subsequent solves "
83 : "until the matrix changes. The dense factorization needs memory "
84 : "quadratic in the number of points, so this can be unfeasible if the "
85 : "linear problem is too big."};
86 : };
87 :
88 : /// Option for skipping resets of the built matrix.
89 1 : struct SkipResets {
90 0 : using type = bool;
91 0 : using group = BuildMatrixOptionsGroup;
92 0 : static constexpr Options::String help = {
93 : "Skip resets of the built matrix. This only has an effect in cases "
94 : "where the operator changes, e.g. between nonlinear-solver iterations. "
95 : "Skipping resets avoids expensive re-building of the matrix, but comes "
96 : "at the cost of less accurate preconditioning and thus potentially more "
97 : "preconditioned iterations. Whether or not this helps convergence "
98 : "overall is highly problem-dependent."};
99 : };
100 :
101 : } // namespace OptionTags
102 :
103 1 : namespace Tags {
104 :
105 : /// Subfile name in the volume data H5 files where the matrix will be stored.
106 1 : struct MatrixSubfileName : db::SimpleTag {
107 0 : using type = std::optional<std::string>;
108 0 : using option_tags = tmpl::list<OptionTags::MatrixSubfileName>;
109 0 : static constexpr bool pass_metavariables = false;
110 0 : static type create_from_options(const type& value) { return value; }
111 : };
112 :
113 : /// Size of the matrix: number of grid points times number of variables
114 1 : struct TotalNumPoints : db::SimpleTag {
115 0 : using type = size_t;
116 : };
117 :
118 : /// Index of the first point in this element
119 1 : struct LocalFirstIndex : db::SimpleTag {
120 0 : using type = size_t;
121 : };
122 :
123 : /// Matrix representation of the linear operator. Stored as sparse matrix. The
124 : /// full matrix has size `total_num_points` by `total_num_points`. Each element
125 : /// only stores the rows corresponding to its grid points (starting at
126 : /// `local_first_index`), but all columns.
127 : template <typename ValueType>
128 1 : struct Matrix : db::SimpleTag {
129 0 : using type = blaze::CompressedMatrix<ValueType>;
130 : };
131 :
132 : /// Whether the matrix on this element was (re)built since it was last sent to
133 : /// the singleton that inverts the full matrix. The singleton only factorizes
134 : /// the full matrix again if this is `true` on any element.
135 1 : struct MatrixChanged : db::SimpleTag {
136 0 : using type = bool;
137 : };
138 :
139 : /// LU factorization of the full matrix, held by the singleton that solves the
140 : /// linear problem directly. It is not serialized, but rebuilt from the matrix
141 : /// slices that the elements send with every solve.
142 : template <typename ValueType>
143 1 : struct LuFactorization : db::SimpleTag {
144 0 : struct type {
145 0 : blaze::DynamicMatrix<ValueType, blaze::columnMajor> lu{};
146 0 : std::vector<blaze::blas_int_t> pivots{};
147 0 : void pup(PUP::er& /*p*/) {}
148 : };
149 : };
150 :
151 : /// Option for enabling direct solve of the linear problem.
152 1 : struct EnableDirectSolve : db::SimpleTag {
153 0 : using type = bool;
154 0 : using option_tags = tmpl::list<OptionTags::EnableDirectSolve>;
155 0 : static constexpr bool pass_metavariables = false;
156 0 : static type create_from_options(const type& value) { return value; }
157 : };
158 :
159 : /// Option for skipping resets of the built matrix.
160 1 : struct SkipResets : db::SimpleTag {
161 0 : using type = bool;
162 0 : using option_tags = tmpl::list<OptionTags::SkipResets>;
163 0 : static constexpr bool pass_metavariables = false;
164 0 : static type create_from_options(const type& value) { return value; }
165 : };
166 :
167 : } // namespace Tags
168 :
169 1 : namespace Actions {
170 :
171 : namespace detail {
172 :
173 : template <typename FieldsTag, typename FixedSourcesTag, typename OperandTag,
174 : typename OperatorAppliedToOperandTag, typename CoordsTag,
175 : typename ArraySectionIdTag>
176 : struct BuildMatrixMetavars {
177 : using iteration_id_tag = Convergence::Tags::IterationId<
178 : LinearSolver::OptionTags::BuildMatrixOptionsGroup>;
179 : using fields_tag = FieldsTag;
180 : using fixed_sources_tag = FixedSourcesTag;
181 : using operand_tag = OperandTag;
182 : using operator_applied_to_operand_tag = OperatorAppliedToOperandTag;
183 : using coords_tag = CoordsTag;
184 : using array_section_id_tag = ArraySectionIdTag;
185 :
186 : using value_type = typename OperatorAppliedToOperandTag::type::value_type;
187 :
188 : struct end_label {};
189 : };
190 :
191 : /// \brief The total number of grid points (size of the matrix) and the index of
192 : /// the first grid point in this element (the offset into the matrix
193 : /// corresponding to this element). The `num_points_per_element` should hold
194 : /// the total number of degrees of freedom for each element, so the number of
195 : /// variables times the number of grid points.
196 : template <size_t Dim>
197 : std::pair<size_t, size_t> total_num_points_and_local_first_index(
198 : const ElementId<Dim>& element_id,
199 : const std::map<ElementId<Dim>, size_t>& num_points_per_element);
200 :
201 : /// \brief The index of the '1' of the unit vector in this element, or
202 : /// std::nullopt if the '1' is in another element.
203 : ///
204 : /// \param iteration_id enumerates all grid points
205 : /// \param local_first_index the index of the first grid point in this element
206 : /// \param local_num_points the number of grid points in this element
207 : std::optional<size_t> local_unit_vector_index(size_t iteration_id,
208 : size_t local_first_index,
209 : size_t local_num_points);
210 :
211 : /// \brief Observe matrix column as volume data.
212 : ///
213 : /// This is the format of the volume data:
214 : ///
215 : /// - Columns of the matrix are enumerated by the observation ID
216 : /// - Rows are enumerated by the order of elements defined by the ElementId
217 : /// in the volume data, by the order of tensor components encoded in the
218 : /// name of the components, and by the contiguous ordering of grid points
219 : /// for each component.
220 : /// - We also observe coordinates so the data can be plotted as volume data.
221 : ///
222 : /// The matrix can be reconstructed from this data in Python with the function
223 : /// `spectre.Elliptic.ReadH5.read_matrix`.
224 : template <typename ParallelComponent, typename OperatorAppliedToOperandTags,
225 : size_t Dim, typename CoordsFrame, typename Metavariables>
226 : void observe_matrix_column(
227 : const size_t column,
228 : const Variables<OperatorAppliedToOperandTags>& operator_applied_to_operand,
229 : const ElementId<Dim>& element_id, const Mesh<Dim>& mesh,
230 : const tnsr::I<DataVector, Dim, CoordsFrame>& coords,
231 : const std::string& subfile_name, const std::string& section_observation_key,
232 : Parallel::GlobalCache<Metavariables>& cache) {
233 : std::vector<TensorComponent> observe_components{};
234 : observe_components.reserve(
235 : Dim +
236 : Variables<
237 : OperatorAppliedToOperandTags>::number_of_independent_components);
238 : for (size_t i = 0; i < Dim; ++i) {
239 : observe_components.emplace_back(
240 : get_output(CoordsFrame{}) + "Coordinates" + coords.component_suffix(i),
241 : coords.get(i));
242 : }
243 : size_t component_i = 0;
244 : tmpl::for_each<OperatorAppliedToOperandTags>([&operator_applied_to_operand,
245 : &observe_components,
246 : &component_i](auto tag_v) {
247 : using tag = tmpl::type_from<decltype(tag_v)>;
248 : const auto& tensor = get<tag>(operator_applied_to_operand);
249 : using TensorType = std::decay_t<decltype(tensor)>;
250 : using VectorType = typename TensorType::type;
251 : using ValueType = typename VectorType::value_type;
252 : for (size_t i = 0; i < tensor.size(); ++i) {
253 : if constexpr (std::is_same_v<ValueType, std::complex<double>>) {
254 : observe_components.emplace_back(
255 : "Re(Variable_" + std::to_string(component_i) + ")",
256 : real(tensor[i]));
257 : observe_components.emplace_back(
258 : "Im(Variable_" + std::to_string(component_i) + ")",
259 : imag(tensor[i]));
260 : } else {
261 : observe_components.emplace_back(
262 : "Variable_" + std::to_string(component_i), tensor[i]);
263 : }
264 : ++component_i;
265 : }
266 : });
267 :
268 : observers::contribute_volume_data<
269 : not Parallel::is_nodegroup_v<ParallelComponent>>(
270 : cache,
271 : observers::ObservationId(static_cast<double>(column),
272 : subfile_name + section_observation_key),
273 : "/" + subfile_name,
274 : Parallel::make_array_component_id<ParallelComponent>(element_id),
275 : ElementVolumeData{element_id, std::move(observe_components), mesh});
276 : }
277 :
278 : /// \brief Register with the volume observer
279 : template <typename ArraySectionIdTag>
280 : struct RegisterWithVolumeObserver {
281 : template <typename ParallelComponent, typename DbTagsList,
282 : typename ArrayIndex>
283 : static std::pair<observers::TypeOfObservation, observers::ObservationKey>
284 : register_info(const db::DataBox<DbTagsList>& box,
285 : const ArrayIndex& /*array_index*/) {
286 : // Get the observation key, or "Unused" if the element does not belong
287 : // to a section with this tag. In the latter case, no observations will
288 : // ever be contributed.
289 : const std::optional<std::string> section_observation_key =
290 : observers::get_section_observation_key<ArraySectionIdTag>(box);
291 : ASSERT(section_observation_key != "Unused",
292 : "The identifier 'Unused' is reserved to indicate that no "
293 : "observations with this key will be contributed. Use a different "
294 : "key, or change the identifier 'Unused' to something else.");
295 : const auto& subfile_name = get<Tags::MatrixSubfileName>(box);
296 : return {
297 : observers::TypeOfObservation::Volume,
298 : observers::ObservationKey(subfile_name.value_or("Unused") +
299 : section_observation_key.value_or("Unused"))};
300 : }
301 : };
302 :
303 : } // namespace detail
304 :
305 : /// \cond
306 : template <typename BuildMatrixMetavars>
307 : struct PrepareBuildMatrix;
308 : template <typename BuildMatrixMetavars>
309 : struct AssembleFullMatrix;
310 : /// \endcond
311 :
312 : /// Dispatch global reduction to get the size of the matrix
313 : template <typename BuildMatrixMetavars>
314 1 : struct CollectTotalNumPoints {
315 : private:
316 0 : using IterationIdTag = typename BuildMatrixMetavars::iteration_id_tag;
317 0 : using OperandTag = typename BuildMatrixMetavars::operand_tag;
318 0 : using ArraySectionIdTag = typename BuildMatrixMetavars::array_section_id_tag;
319 0 : using value_type = typename BuildMatrixMetavars::value_type;
320 :
321 : public:
322 0 : using simple_tags =
323 : tmpl::list<Tags::TotalNumPoints, Tags::LocalFirstIndex, IterationIdTag,
324 : OperandTag, Tags::Matrix<value_type>, Tags::MatrixChanged>;
325 0 : using compute_tags = tmpl::list<>;
326 0 : using const_global_cache_tags =
327 : tmpl::list<logging::Tags::Verbosity<OptionTags::BuildMatrixOptionsGroup>>;
328 :
329 : template <typename DbTags, typename... InboxTags, typename Metavariables,
330 : size_t Dim, typename ActionList, typename ParallelComponent>
331 0 : static Parallel::iterable_action_return_t apply(
332 : db::DataBox<DbTags>& box,
333 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
334 : Parallel::GlobalCache<Metavariables>& cache,
335 : const ElementId<Dim>& element_id, const ActionList /*meta*/,
336 : const ParallelComponent* const /*meta*/) {
337 : // Skip everything on elements that are not part of the section
338 : if constexpr (not std::is_same_v<ArraySectionIdTag, void>) {
339 : if (not db::get<Parallel::Tags::Section<ParallelComponent,
340 : ArraySectionIdTag>>(box)
341 : .has_value()) {
342 : constexpr size_t last_action_index = tmpl::index_of<
343 : ActionList,
344 : ::Actions::Label<typename BuildMatrixMetavars::end_label>>::value;
345 : return {Parallel::AlgorithmExecution::Continue, last_action_index + 1};
346 : }
347 : }
348 : if (db::get<Tags::TotalNumPoints>(box) != 0) {
349 : // We have built the matrix already, so we can skip ahead
350 : constexpr size_t assemble_matrix_index =
351 : tmpl::index_of<ActionList,
352 : AssembleFullMatrix<BuildMatrixMetavars>>::value;
353 : return {Parallel::AlgorithmExecution::Continue, assemble_matrix_index};
354 : }
355 : db::mutate<Tags::TotalNumPoints, OperandTag>(
356 : [](const auto total_num_points, const auto operand,
357 : const size_t num_points) {
358 : // Set num points to zero until we know the total number of points
359 : *total_num_points = 0;
360 : // Size operand and fill with zeros
361 : operand->initialize(num_points, 0.);
362 : },
363 : make_not_null(&box),
364 : get<domain::Tags::Mesh<Dim>>(box).number_of_grid_points());
365 : // Collect the total number of grid points
366 : auto& section = Parallel::get_section<ParallelComponent, ArraySectionIdTag>(
367 : make_not_null(&box));
368 : Parallel::contribute_to_reduction<PrepareBuildMatrix<BuildMatrixMetavars>>(
369 : Parallel::ReductionData<
370 : Parallel::ReductionDatum<size_t, funcl::AssertEqual<>>,
371 : Parallel::ReductionDatum<std::map<ElementId<Dim>, size_t>,
372 : funcl::Merge<>>>{
373 : element_id.grid_index(),
374 : std::map<ElementId<Dim>, size_t>{
375 : std::make_pair(element_id, get<OperandTag>(box).size())}},
376 : Parallel::get_parallel_component<ParallelComponent>(cache)[element_id],
377 : Parallel::get_parallel_component<ParallelComponent>(cache),
378 : make_not_null(§ion));
379 : // Pause the algorithm for now. The reduction will be broadcast to the next
380 : // action which is responsible for restarting the algorithm.
381 : return {Parallel::AlgorithmExecution::Pause, std::nullopt};
382 : }
383 : };
384 :
385 : /// Receive the reduction and initialize the algorithm
386 : template <typename BuildMatrixMetavars>
387 1 : struct PrepareBuildMatrix {
388 : private:
389 0 : using IterationIdTag = typename BuildMatrixMetavars::iteration_id_tag;
390 0 : using OperandTag = typename BuildMatrixMetavars::operand_tag;
391 0 : using value_type = typename BuildMatrixMetavars::value_type;
392 :
393 : public:
394 : template <typename ParallelComponent, typename DbTagsList,
395 : typename Metavariables, size_t Dim>
396 0 : static void apply(
397 : db::DataBox<DbTagsList>& box, Parallel::GlobalCache<Metavariables>& cache,
398 : const ElementId<Dim>& element_id, const size_t grid_index,
399 : const std::map<ElementId<Dim>, size_t>& num_points_per_element) {
400 : if (grid_index != element_id.grid_index()) {
401 : return;
402 : }
403 : const auto [total_num_points, local_first_index] =
404 : detail::total_num_points_and_local_first_index(element_id,
405 : num_points_per_element);
406 : if (get<logging::Tags::Verbosity<OptionTags::BuildMatrixOptionsGroup>>(
407 : box) >= Verbosity::Quiet and
408 : local_first_index == 0) {
409 : Parallel::printf(
410 : "Building explicit matrix representation of size %zu x %zu.\n",
411 : total_num_points, total_num_points);
412 : }
413 : db::mutate<Tags::TotalNumPoints, Tags::LocalFirstIndex, IterationIdTag,
414 : Tags::Matrix<value_type>, Tags::MatrixChanged>(
415 : [captured_total_num_points = total_num_points,
416 : captured_local_first_index = local_first_index,
417 : local_num_points = num_points_per_element.at(element_id)](
418 : const auto stored_total_num_points,
419 : const auto stored_local_first_index, const auto iteration_id,
420 : const gsl::not_null<blaze::CompressedMatrix<value_type>*> matrix,
421 : const gsl::not_null<bool*> matrix_changed) {
422 : *stored_total_num_points = captured_total_num_points;
423 : *stored_local_first_index = captured_local_first_index;
424 : *iteration_id = 0;
425 : matrix->clear();
426 : matrix->resize(local_num_points, captured_total_num_points);
427 : *matrix_changed = true;
428 : },
429 : make_not_null(&box));
430 : // Proceed with algorithm
431 : Parallel::get_parallel_component<ParallelComponent>(cache)[element_id]
432 : .perform_algorithm(true);
433 : }
434 : };
435 :
436 : /// Set the operand to a unit vector with a '1' at the current grid point.
437 : /// Applying the operator to this operand gives a column of the matrix. We jump
438 : /// back to this until we have iterated over all grid points.
439 : template <typename BuildMatrixMetavars>
440 1 : struct SetUnitVector {
441 : private:
442 0 : using IterationIdTag = typename BuildMatrixMetavars::iteration_id_tag;
443 0 : using OperandTag = typename BuildMatrixMetavars::operand_tag;
444 :
445 : public:
446 : template <typename DbTags, typename... InboxTags, typename Metavariables,
447 : typename ArrayIndex, typename ActionList,
448 : typename ParallelComponent>
449 0 : static Parallel::iterable_action_return_t apply(
450 : db::DataBox<DbTags>& box,
451 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
452 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
453 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
454 : const ParallelComponent* const /*meta*/) {
455 : const std::optional<size_t> local_unit_vector_index =
456 : detail::local_unit_vector_index(get<IterationIdTag>(box),
457 : get<Tags::LocalFirstIndex>(box),
458 : get<OperandTag>(box).size());
459 : if (local_unit_vector_index.has_value()) {
460 : db::mutate<OperandTag>(
461 : [&local_unit_vector_index](const auto operand) {
462 : operand->data()[*local_unit_vector_index] = 1.;
463 : },
464 : make_not_null(&box));
465 : }
466 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
467 : }
468 : };
469 :
470 : // --- Linear operator will be applied to the unit vector here ---
471 :
472 : /// Write result out to disk, reset operand back to zero, and keep iterating
473 : template <typename BuildMatrixMetavars>
474 1 : struct StoreMatrixColumn {
475 : private:
476 0 : using IterationIdTag = typename BuildMatrixMetavars::iteration_id_tag;
477 0 : using OperandTag = typename BuildMatrixMetavars::operand_tag;
478 0 : using OperatorAppliedToOperandTag =
479 : typename BuildMatrixMetavars::operator_applied_to_operand_tag;
480 0 : using CoordsTag = typename BuildMatrixMetavars::coords_tag;
481 0 : using ArraySectionIdTag = typename BuildMatrixMetavars::array_section_id_tag;
482 0 : using value_type = typename BuildMatrixMetavars::value_type;
483 :
484 : public:
485 0 : using const_global_cache_tags = tmpl::list<Tags::MatrixSubfileName>;
486 :
487 : template <typename DbTags, typename... InboxTags, typename Metavariables,
488 : size_t Dim, typename ActionList, typename ParallelComponent>
489 0 : static Parallel::iterable_action_return_t apply(
490 : db::DataBox<DbTags>& box,
491 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
492 : Parallel::GlobalCache<Metavariables>& cache,
493 : const ElementId<Dim>& element_id, const ActionList /*meta*/,
494 : const ParallelComponent* const /*meta*/) {
495 : const size_t iteration_id = get<IterationIdTag>(box);
496 : const size_t local_first_index = get<Tags::LocalFirstIndex>(box);
497 : if (get<logging::Tags::Verbosity<OptionTags::BuildMatrixOptionsGroup>>(
498 : box) >= Verbosity::Verbose and
499 : local_first_index == 0 and (iteration_id + 1) % 100 == 0) {
500 : Parallel::printf("Column %zu / %zu done.\n", iteration_id + 1,
501 : get<Tags::TotalNumPoints>(box));
502 : }
503 : // This is the result of applying the linear operator to the unit vector. It
504 : // is a column of the operator matrix.
505 : const auto& operator_applied_to_operand =
506 : get<OperatorAppliedToOperandTag>(box);
507 : // Store in sparse matrix
508 : db::mutate<Tags::Matrix<value_type>>(
509 : [&operator_applied_to_operand, &iteration_id](
510 : const gsl::not_null<blaze::CompressedMatrix<value_type>*> matrix) {
511 : for (size_t i = 0; i < operator_applied_to_operand.size(); ++i) {
512 : const auto value = operator_applied_to_operand.data()[i];
513 : if (not equal_within_roundoff(value, 0.)) {
514 : matrix->insert(i, iteration_id, value);
515 : }
516 : }
517 : },
518 : make_not_null(&box));
519 : // Write it out to disk
520 : if (const auto& subfile_name = get<Tags::MatrixSubfileName>(box);
521 : subfile_name.has_value()) {
522 : detail::observe_matrix_column<ParallelComponent>(
523 : iteration_id, operator_applied_to_operand, element_id,
524 : get<domain::Tags::Mesh<Dim>>(box), get<CoordsTag>(box), *subfile_name,
525 : *observers::get_section_observation_key<ArraySectionIdTag>(box),
526 : cache);
527 : }
528 : // Reset operand to zero
529 : const std::optional<size_t> local_unit_vector_index =
530 : detail::local_unit_vector_index(iteration_id, local_first_index,
531 : get<OperandTag>(box).size());
532 : if (local_unit_vector_index.has_value()) {
533 : db::mutate<OperandTag>(
534 : [&local_unit_vector_index](const auto operand) {
535 : operand->data()[*local_unit_vector_index] = 0.;
536 : },
537 : make_not_null(&box));
538 : }
539 : // Keep iterating
540 : db::mutate<IterationIdTag>(
541 : [](const auto local_iteration_id) { ++(*local_iteration_id); },
542 : make_not_null(&box));
543 : if (get<IterationIdTag>(box) < get<Tags::TotalNumPoints>(box)) {
544 : constexpr size_t set_unit_vector_index =
545 : tmpl::index_of<ActionList, SetUnitVector<BuildMatrixMetavars>>::value;
546 : return {Parallel::AlgorithmExecution::Continue, set_unit_vector_index};
547 : }
548 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
549 : }
550 : };
551 :
552 : /// Adds the LU factorization of the full matrix to the singleton's DataBox
553 : template <typename BuildMatrixMetavars>
554 1 : struct InitializeSingleton {
555 0 : using value_type = typename BuildMatrixMetavars::value_type;
556 0 : using simple_tags = tmpl::list<Tags::LuFactorization<value_type>>;
557 :
558 : template <typename DbTags, typename... InboxTags, typename Metavariables,
559 : typename ArrayIndex, typename ActionList,
560 : typename ParallelComponent>
561 0 : static Parallel::iterable_action_return_t apply(
562 : db::DataBox<DbTags>& /*box*/,
563 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
564 : Parallel::GlobalCache<Metavariables>& /*cache*/,
565 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
566 : const ParallelComponent* const /*meta*/) {
567 : return {Parallel::AlgorithmExecution::Pause, std::nullopt};
568 : }
569 : };
570 :
571 : template <typename Metavariables, typename BuildMatrixMetavars>
572 0 : struct BuildMatrixSingleton {
573 0 : using chare_type = Parallel::Algorithms::Singleton;
574 0 : static constexpr bool checkpoint_data = true;
575 0 : using const_global_cache_tags = tmpl::list<>;
576 0 : using metavariables = Metavariables;
577 0 : using phase_dependent_action_list = tmpl::list<Parallel::PhaseActions<
578 : Parallel::Phase::Initialization,
579 : tmpl::list<InitializeSingleton<BuildMatrixMetavars>>>>;
580 0 : using simple_tags_from_options = Parallel::get_simple_tags_from_options<
581 : Parallel::get_initialization_actions_list<phase_dependent_action_list>>;
582 :
583 0 : static void execute_next_phase(
584 : const Parallel::Phase next_phase,
585 : Parallel::CProxy_GlobalCache<Metavariables>& global_cache) {
586 : auto& local_cache = *Parallel::local_branch(global_cache);
587 : Parallel::get_parallel_component<BuildMatrixSingleton>(local_cache)
588 : .start_phase(next_phase);
589 : }
590 : };
591 :
592 : template <typename ElementArrayComponent, typename BuildMatrixMetavars>
593 : struct InvertMatrix;
594 :
595 : template <typename BuildMatrixMetavars>
596 0 : struct AssembleFullMatrix {
597 : private:
598 0 : using value_type = typename BuildMatrixMetavars::value_type;
599 0 : using ArraySectionIdTag = typename BuildMatrixMetavars::array_section_id_tag;
600 0 : using FixedSourcesTag = typename BuildMatrixMetavars::fixed_sources_tag;
601 0 : using ReductionType = std::pair<blaze::CompressedMatrix<value_type>,
602 : typename FixedSourcesTag::type>;
603 :
604 : public:
605 0 : using const_global_cache_tags = tmpl::list<Tags::EnableDirectSolve>;
606 :
607 : template <typename DbTags, typename... InboxTags, typename Metavariables,
608 : size_t Dim, typename ActionList, typename ParallelComponent>
609 0 : static Parallel::iterable_action_return_t apply(
610 : db::DataBox<DbTags>& box,
611 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
612 : Parallel::GlobalCache<Metavariables>& cache,
613 : const ElementId<Dim>& element_id, const ActionList /*meta*/,
614 : const ParallelComponent* const /*meta*/) {
615 : if (not db::get<Tags::EnableDirectSolve>(box)) {
616 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
617 : }
618 : auto& section = Parallel::get_section<ParallelComponent, ArraySectionIdTag>(
619 : make_not_null(&box));
620 : Parallel::contribute_to_reduction<
621 : InvertMatrix<ParallelComponent, BuildMatrixMetavars>>(
622 : Parallel::ReductionData<
623 : Parallel::ReductionDatum<size_t, funcl::AssertEqual<>>,
624 : Parallel::ReductionDatum<bool, funcl::Or<>>,
625 : Parallel::ReductionDatum<std::map<ElementId<Dim>, ReductionType>,
626 : funcl::Merge<>>>{
627 : element_id.grid_index(), get<Tags::MatrixChanged>(box),
628 : std::map<ElementId<Dim>, ReductionType>{std::make_pair(
629 : element_id, ReductionType{get<Tags::Matrix<value_type>>(box),
630 : get<FixedSourcesTag>(box)})}},
631 : Parallel::get_parallel_component<ParallelComponent>(cache)[element_id],
632 : Parallel::get_parallel_component<
633 : BuildMatrixSingleton<Metavariables, BuildMatrixMetavars>>(cache),
634 : make_not_null(§ion));
635 : db::mutate<Tags::MatrixChanged>(
636 : [](const gsl::not_null<bool*> matrix_changed) {
637 : *matrix_changed = false;
638 : },
639 : make_not_null(&box));
640 : // Pause the algorithm for now. The reduction will be received by the next
641 : // action which is responsible for restarting the algorithm.
642 : return {Parallel::AlgorithmExecution::Pause, std::nullopt};
643 : }
644 : };
645 :
646 : template <typename BuildMatrixMetavars>
647 : struct StoreSolution;
648 :
649 : template <typename ElementArrayComponent, typename BuildMatrixMetavars>
650 0 : struct InvertMatrix {
651 : private:
652 0 : using value_type = typename BuildMatrixMetavars::value_type;
653 0 : using FixedSourcesTag = typename BuildMatrixMetavars::fixed_sources_tag;
654 0 : using ReductionType = std::pair<blaze::CompressedMatrix<value_type>,
655 : typename FixedSourcesTag::type>;
656 :
657 : public:
658 : template <typename ParallelComponent, typename DbTagsList,
659 : typename Metavariables, typename ArrayIndex, size_t Dim>
660 0 : static void apply(db::DataBox<DbTagsList>& box,
661 : Parallel::GlobalCache<Metavariables>& cache,
662 : const ArrayIndex& /*array_index*/, const size_t grid_index,
663 : const bool matrix_changed,
664 : const std::map<ElementId<Dim>, ReductionType>&
665 : matrix_and_sources_slices) {
666 : // Assemble and factorize the full matrix if it changed since we last
667 : // factorized it. The factorization is not serialized, so it is also rebuilt
668 : // after a restart. Possible performance optimization: Don't send the matrix
669 : // slices again if they haven't changed.
670 : const size_t total_num_points =
671 : matrix_and_sources_slices.begin()->second.first.columns();
672 : auto& [lu, pivots] =
673 : db::get_mutable_reference<Tags::LuFactorization<value_type>>(
674 : make_not_null(&box));
675 : const bool factorize = matrix_changed or lu.rows() == 0;
676 : if (factorize) {
677 : lu.resize(total_num_points, total_num_points);
678 : pivots.resize(total_num_points);
679 : }
680 : ASSERT(lu.rows() == total_num_points,
681 : "The factorized matrix has size "
682 : << lu.rows() << " but the source has size " << total_num_points
683 : << ". The matrix changed but was not factorized again.");
684 : blaze::DynamicVector<value_type> solution(total_num_points);
685 : size_t i = 0;
686 : for (const auto& [element_id, slice] : matrix_and_sources_slices) {
687 : const auto& [matrix_slice, source_slice] = slice;
688 : ASSERT(matrix_slice.rows() == source_slice.size(),
689 : "Matrix and source slice have different sizes.");
690 : const size_t slice_length = matrix_slice.rows();
691 : if (factorize) {
692 : blaze::submatrix(lu, i, 0, slice_length, total_num_points) =
693 : matrix_slice;
694 : }
695 : std::copy_n(source_slice.data(), slice_length, solution.data() + i);
696 : i += slice_length;
697 : }
698 : ASSERT(i == total_num_points, "The matrix is not fully assembled. Expected "
699 : << total_num_points << " rows, got "
700 : << i);
701 : if (factorize) {
702 : if (get<logging::Tags::Verbosity<OptionTags::BuildMatrixOptionsGroup>>(
703 : box) >= Verbosity::Quiet) {
704 : Parallel::printf("Factorizing matrix of size %zu x %zu\n",
705 : total_num_points, total_num_points);
706 : }
707 : const auto size = static_cast<blaze::blas_int_t>(total_num_points);
708 : blaze::blas_int_t info = 0;
709 : blaze::getrf(size, size, lu.data(),
710 : static_cast<blaze::blas_int_t>(lu.spacing()), pivots.data(),
711 : &info);
712 : if (info != 0) {
713 : ERROR("LU factorization of the matrix of size "
714 : << total_num_points << " x " << total_num_points
715 : << " failed with info = " << info
716 : << ". A positive value means the matrix is singular.");
717 : }
718 : }
719 : // Solve the equations Ax = b using the LU factorization
720 : blaze::getrs(lu, solution, 'N', pivots.data());
721 : // Broadcast the solution to the elements
722 : Parallel::simple_action<StoreSolution<BuildMatrixMetavars>>(
723 : Parallel::get_parallel_component<ElementArrayComponent>(cache),
724 : grid_index, solution);
725 : }
726 : };
727 :
728 : template <typename BuildMatrixMetavars>
729 0 : struct StoreSolution {
730 : private:
731 0 : using value_type = typename BuildMatrixMetavars::value_type;
732 0 : using FieldsTag = typename BuildMatrixMetavars::fields_tag;
733 0 : using FixedSourcesTag = typename BuildMatrixMetavars::fixed_sources_tag;
734 :
735 : public:
736 : template <typename ParallelComponent, typename DbTagsList,
737 : typename Metavariables, size_t Dim>
738 0 : static void apply(db::DataBox<DbTagsList>& box,
739 : Parallel::GlobalCache<Metavariables>& cache,
740 : const ElementId<Dim>& element_id, const size_t grid_index,
741 : const blaze::DynamicVector<value_type>& solution) {
742 : if (grid_index != element_id.grid_index()) {
743 : return;
744 : }
745 : const size_t local_first_index = db::get<Tags::LocalFirstIndex>(box);
746 : db::mutate<
747 : FieldsTag,
748 : db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, FieldsTag>>(
749 : [&solution, &local_first_index](const auto fields,
750 : const auto operator_applied_to_fields,
751 : const auto& fixed_sources) {
752 : std::copy_n(
753 : solution.begin() + static_cast<ptrdiff_t>(local_first_index),
754 : fields->size(), fields->data());
755 : // The linear problem is solved now, so Ax = b
756 : *operator_applied_to_fields = fixed_sources;
757 : },
758 : make_not_null(&box), db::get<FixedSourcesTag>(box));
759 : // Proceed with algorithm
760 : Parallel::get_parallel_component<ParallelComponent>(cache)[element_id]
761 : .perform_algorithm(true);
762 : }
763 : };
764 :
765 : template <typename BuildMatrixMetavars>
766 0 : struct ProjectBuildMatrix : tt::ConformsTo<::amr::protocols::Projector> {
767 : private:
768 0 : using IterationIdTag = typename BuildMatrixMetavars::iteration_id_tag;
769 0 : using OperandTag = typename BuildMatrixMetavars::operand_tag;
770 0 : using value_type = typename BuildMatrixMetavars::value_type;
771 :
772 : public:
773 0 : using return_tags = tmpl::list<Tags::TotalNumPoints, Tags::Matrix<value_type>,
774 : Tags::LocalFirstIndex, IterationIdTag,
775 : OperandTag, Tags::MatrixChanged>;
776 0 : using argument_tags = tmpl::list<>;
777 :
778 : template <typename... AmrData>
779 0 : static void apply(
780 : const gsl::not_null<size_t*> total_num_points,
781 : const gsl::not_null<blaze::CompressedMatrix<value_type>*> matrix,
782 : const AmrData&... /*amr_data*/) {
783 : // Reset the built matrix when AMR changes the grid.
784 : // In case AMR is configured to keep coarse grids then the coarse-grid
785 : // elements don't run projectors during AMR, so they are not reset and just
786 : // keep the built matrix. This is good because the coarse grids are
787 : // unchanged and so the matrix is still valid (and can be used as bottom
788 : // solver in multigrid).
789 : *total_num_points = 0;
790 : matrix->clear();
791 : }
792 : };
793 :
794 : /// Dispatch global reduction to get the size of the matrix
795 : template <typename BuildMatrixMetavars>
796 1 : struct ResetBuiltMatrix {
797 : private:
798 0 : using value_type = typename BuildMatrixMetavars::value_type;
799 0 : using ArraySectionIdTag = typename BuildMatrixMetavars::array_section_id_tag;
800 :
801 : public:
802 0 : using const_global_cache_tags = tmpl::list<Tags::SkipResets>;
803 :
804 : template <typename DbTags, typename... InboxTags, typename Metavariables,
805 : size_t Dim, typename ActionList, typename ParallelComponent>
806 0 : static Parallel::iterable_action_return_t apply(
807 : db::DataBox<DbTags>& box,
808 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
809 : Parallel::GlobalCache<Metavariables>& /*cache*/,
810 : const ElementId<Dim>& /*array_index*/, const ActionList /*meta*/,
811 : const ParallelComponent* const /*meta*/) {
812 : // Skip on elements that are not part of the section
813 : if constexpr (not std::is_same_v<ArraySectionIdTag, void>) {
814 : if (not db::get<Parallel::Tags::Section<ParallelComponent,
815 : ArraySectionIdTag>>(box)
816 : .has_value()) {
817 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
818 : }
819 : }
820 : if (db::get<Tags::SkipResets>(box)) {
821 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
822 : }
823 : db::mutate<Tags::TotalNumPoints, Tags::Matrix<value_type>>(
824 : [](const auto total_num_points, const auto matrix) {
825 : *total_num_points = 0;
826 : matrix->clear();
827 : },
828 : make_not_null(&box));
829 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
830 : }
831 : };
832 :
833 : /*!
834 : * \brief Build the explicit matrix representation of the linear operator and
835 : * optionally solve the problem directly with it.
836 : *
837 : * The matrix is built by applying the linear operator to unit vectors, so it
838 : * takes as many operator applications as there are points. Each element stores
839 : * the rows of the matrix that correspond to its points. The matrix is built
840 : * once and kept until it is reset (see `reset_actions`) or AMR changes the
841 : * grid.
842 : *
843 : * If direct solve is enabled, the elements send their rows of the matrix and
844 : * the source to a singleton. The singleton assembles the full matrix, computes
845 : * its LU factorization, and broadcasts the solution back to the elements. The
846 : * LU factorization is reused for subsequent solves until the matrix changes,
847 : * so these solves only cost two triangular solves on the singleton. This is
848 : * useful as bottom solver in a multigrid setting, where the coarsest grid is
849 : * small and the same linear operator is inverted in every V-cycle.
850 : *
851 : * Add the `actions` to the action list to build the matrix. The
852 : * `ApplyOperatorActions` template parameter are the actions that apply the
853 : * linear operator to the `OperandTag`. Also add the `amr_projectors` to the
854 : * list of AMR projectors and the `register_actions`
855 : *
856 : * \tparam FieldsTag The solution will be stored here (if direct solve is
857 : * enabled).
858 : * \tparam FixedSourcesTag The source `b` in the problem `Ax = b`. Only used
859 : * if direct solve is enabled.
860 : * \tparam OperandTag Will be set to unit vectors, with the '1' moving through
861 : * all points over the course of the iteration.
862 : * \tparam OperatorAppliedToOperandTag Where the `ApplyOperatorActions` store
863 : * the result of applying the linear operator to the `OperandTag`.
864 : * \tparam CoordsTag The tag of the coordinates observed alongside the matrix
865 : * for volume data visualization.
866 : * \tparam ArraySectionIdTag Can identify a subset of elements that this
867 : * algorithm should run over, e.g. in a multigrid setting.
868 : */
869 : template <typename FieldsTag, typename FixedSourcesTag, typename OperandTag,
870 : typename OperatorAppliedToOperandTag, typename CoordsTag,
871 : typename ArraySectionIdTag = void>
872 1 : struct BuildMatrix {
873 : private:
874 0 : using BuildMatrixMetavars =
875 : detail::BuildMatrixMetavars<FieldsTag, FixedSourcesTag, OperandTag,
876 : OperatorAppliedToOperandTag, CoordsTag,
877 : ArraySectionIdTag>;
878 : static_assert(
879 : std::is_same_v<FieldsTag, OperandTag>,
880 : "The operand and the fields tags must be the same. This is just so that "
881 : "at the end of the algorithm the operator is applied to the solution "
882 : "fields. This restriction can be lifted if needed by copying the fields "
883 : "into the operand and back.");
884 :
885 : public:
886 : template <typename Metavariables>
887 0 : using component_list =
888 : tmpl::list<BuildMatrixSingleton<Metavariables, BuildMatrixMetavars>>;
889 :
890 : template <typename ApplyOperatorActions>
891 0 : using actions =
892 : tmpl::list<CollectTotalNumPoints<BuildMatrixMetavars>,
893 : // PrepareBuildMatrix is called on reduction broadcast
894 : SetUnitVector<BuildMatrixMetavars>, ApplyOperatorActions,
895 : StoreMatrixColumn<BuildMatrixMetavars>,
896 : // Algorithm iterates until matrix is complete, then proceeds
897 : // below
898 : AssembleFullMatrix<BuildMatrixMetavars>,
899 : // Apply operator to the solution. Note that FieldsTag and
900 : // OperandTag must be the same for this (see assert above).
901 : // Note also that this operator application is not needed in
902 : // most cases, as the linear problem is solved exactly and
903 : // therefore Ax = b is set in 'StoreSolution'. However, this
904 : // operator application is needed if the operator changes
905 : // between nonlinear solver iterations but we skip the reset,
906 : // so the matrix represents the old operator. It's just one
907 : // more operator application, so we keep it always enabled for
908 : // now.
909 : ApplyOperatorActions,
910 : ::Actions::Label<typename BuildMatrixMetavars::end_label>>;
911 :
912 0 : using amr_projectors = tmpl::list<ProjectBuildMatrix<BuildMatrixMetavars>>;
913 :
914 : /// Add to the register phase to enable observations
915 1 : using register_actions = tmpl::list<observers::Actions::RegisterWithObservers<
916 : detail::RegisterWithVolumeObserver<ArraySectionIdTag>>>;
917 :
918 : /// Add to the action list to reset the matrix
919 1 : using reset_actions = tmpl::list<ResetBuiltMatrix<BuildMatrixMetavars>>;
920 : };
921 :
922 : } // namespace Actions
923 : } // namespace LinearSolver
|
