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