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/ObserverComponent.hpp"
29 : #include "IO/Observer/VolumeActions.hpp"
30 : #include "NumericalAlgorithms/Convergence/Tags.hpp"
31 : #include "Parallel/AlgorithmExecution.hpp"
32 : #include "Parallel/GetSection.hpp"
33 : #include "Parallel/Phase.hpp"
34 : #include "Parallel/Printf/Printf.hpp"
35 : #include "Parallel/Reduction.hpp"
36 : #include "Parallel/Section.hpp"
37 : #include "Parallel/Tags/Section.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 : auto& local_observer = *Parallel::local_branch(
245 : Parallel::get_parallel_component<observers::Observer<Metavariables>>(
246 : cache));
247 : Parallel::simple_action<observers::Actions::ContributeVolumeData>(
248 : local_observer,
249 : observers::ObservationId(static_cast<double>(column),
250 : subfile_name + section_observation_key),
251 : "/" + subfile_name,
252 : Parallel::make_array_component_id<ParallelComponent>(element_id),
253 : ElementVolumeData{element_id, std::move(observe_components), mesh});
254 : }
255 :
256 : /// \brief Register with the volume observer
257 : template <typename ArraySectionIdTag>
258 : struct RegisterWithVolumeObserver {
259 : template <typename ParallelComponent, typename DbTagsList,
260 : typename ArrayIndex>
261 : static std::pair<observers::TypeOfObservation, observers::ObservationKey>
262 : register_info(const db::DataBox<DbTagsList>& box,
263 : const ArrayIndex& /*array_index*/) {
264 : // Get the observation key, or "Unused" if the element does not belong
265 : // to a section with this tag. In the latter case, no observations will
266 : // ever be contributed.
267 : const std::optional<std::string> section_observation_key =
268 : observers::get_section_observation_key<ArraySectionIdTag>(box);
269 : ASSERT(section_observation_key != "Unused",
270 : "The identifier 'Unused' is reserved to indicate that no "
271 : "observations with this key will be contributed. Use a different "
272 : "key, or change the identifier 'Unused' to something else.");
273 : const auto& subfile_name = get<Tags::MatrixSubfileName>(box);
274 : return {
275 : observers::TypeOfObservation::Volume,
276 : observers::ObservationKey(subfile_name.value_or("Unused") +
277 : section_observation_key.value_or("Unused"))};
278 : }
279 : };
280 :
281 : } // namespace detail
282 :
283 : /// \cond
284 : template <typename BuildMatrixMetavars>
285 : struct PrepareBuildMatrix;
286 : template <typename BuildMatrixMetavars>
287 : struct AssembleFullMatrix;
288 : /// \endcond
289 :
290 : /// Dispatch global reduction to get the size of the matrix
291 : template <typename BuildMatrixMetavars>
292 1 : struct CollectTotalNumPoints {
293 : private:
294 0 : using IterationIdTag = typename BuildMatrixMetavars::iteration_id_tag;
295 0 : using OperandTag = typename BuildMatrixMetavars::operand_tag;
296 0 : using ArraySectionIdTag = typename BuildMatrixMetavars::array_section_id_tag;
297 0 : using value_type = typename BuildMatrixMetavars::value_type;
298 :
299 : public:
300 0 : using simple_tags =
301 : tmpl::list<Tags::TotalNumPoints, Tags::LocalFirstIndex, IterationIdTag,
302 : OperandTag, Tags::Matrix<value_type>>;
303 0 : using compute_tags = tmpl::list<>;
304 0 : using const_global_cache_tags =
305 : tmpl::list<logging::Tags::Verbosity<OptionTags::BuildMatrixOptionsGroup>>;
306 :
307 : template <typename DbTags, typename... InboxTags, typename Metavariables,
308 : size_t Dim, typename ActionList, typename ParallelComponent>
309 0 : static Parallel::iterable_action_return_t apply(
310 : db::DataBox<DbTags>& box,
311 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
312 : Parallel::GlobalCache<Metavariables>& cache,
313 : const ElementId<Dim>& element_id, const ActionList /*meta*/,
314 : const ParallelComponent* const /*meta*/) {
315 : // Skip everything on elements that are not part of the section
316 : if constexpr (not std::is_same_v<ArraySectionIdTag, void>) {
317 : if (not db::get<Parallel::Tags::Section<ParallelComponent,
318 : ArraySectionIdTag>>(box)
319 : .has_value()) {
320 : constexpr size_t last_action_index = tmpl::index_of<
321 : ActionList,
322 : ::Actions::Label<typename BuildMatrixMetavars::end_label>>::value;
323 : return {Parallel::AlgorithmExecution::Continue, last_action_index + 1};
324 : }
325 : }
326 : if (db::get<Tags::TotalNumPoints>(box) != 0) {
327 : // We have built the matrix already, so we can skip ahead
328 : constexpr size_t assemble_matrix_index =
329 : tmpl::index_of<ActionList,
330 : AssembleFullMatrix<BuildMatrixMetavars>>::value;
331 : return {Parallel::AlgorithmExecution::Continue, assemble_matrix_index};
332 : }
333 : db::mutate<Tags::TotalNumPoints, OperandTag>(
334 : [](const auto total_num_points, const auto operand,
335 : const size_t num_points) {
336 : // Set num points to zero until we know the total number of points
337 : *total_num_points = 0;
338 : // Size operand and fill with zeros
339 : operand->initialize(num_points, 0.);
340 : },
341 : make_not_null(&box),
342 : get<domain::Tags::Mesh<Dim>>(box).number_of_grid_points());
343 : // Collect the total number of grid points
344 : auto& section = Parallel::get_section<ParallelComponent, ArraySectionIdTag>(
345 : make_not_null(&box));
346 : Parallel::contribute_to_reduction<PrepareBuildMatrix<BuildMatrixMetavars>>(
347 : Parallel::ReductionData<
348 : Parallel::ReductionDatum<size_t, funcl::AssertEqual<>>,
349 : Parallel::ReductionDatum<std::map<ElementId<Dim>, size_t>,
350 : funcl::Merge<>>>{
351 : element_id.grid_index(),
352 : std::map<ElementId<Dim>, size_t>{
353 : std::make_pair(element_id, get<OperandTag>(box).size())}},
354 : Parallel::get_parallel_component<ParallelComponent>(cache)[element_id],
355 : Parallel::get_parallel_component<ParallelComponent>(cache),
356 : make_not_null(§ion));
357 : // Pause the algorithm for now. The reduction will be broadcast to the next
358 : // action which is responsible for restarting the algorithm.
359 : return {Parallel::AlgorithmExecution::Pause, std::nullopt};
360 : }
361 : };
362 :
363 : /// Receive the reduction and initialize the algorithm
364 : template <typename BuildMatrixMetavars>
365 1 : struct PrepareBuildMatrix {
366 : private:
367 0 : using IterationIdTag = typename BuildMatrixMetavars::iteration_id_tag;
368 0 : using OperandTag = typename BuildMatrixMetavars::operand_tag;
369 0 : using value_type = typename BuildMatrixMetavars::value_type;
370 :
371 : public:
372 : template <typename ParallelComponent, typename DbTagsList,
373 : typename Metavariables, size_t Dim>
374 0 : static void apply(
375 : db::DataBox<DbTagsList>& box, Parallel::GlobalCache<Metavariables>& cache,
376 : const ElementId<Dim>& element_id, const size_t grid_index,
377 : const std::map<ElementId<Dim>, size_t>& num_points_per_element) {
378 : if (grid_index != element_id.grid_index()) {
379 : return;
380 : }
381 : const auto [total_num_points, local_first_index] =
382 : detail::total_num_points_and_local_first_index(element_id,
383 : num_points_per_element);
384 : if (get<logging::Tags::Verbosity<OptionTags::BuildMatrixOptionsGroup>>(
385 : box) >= Verbosity::Quiet and
386 : local_first_index == 0) {
387 : Parallel::printf(
388 : "Building explicit matrix representation of size %zu x %zu.\n",
389 : total_num_points, total_num_points);
390 : }
391 : db::mutate<Tags::TotalNumPoints, Tags::LocalFirstIndex, IterationIdTag,
392 : Tags::Matrix<value_type>>(
393 : [captured_total_num_points = total_num_points,
394 : captured_local_first_index = local_first_index,
395 : local_num_points = num_points_per_element.at(element_id)](
396 : const auto stored_total_num_points,
397 : const auto stored_local_first_index, const auto iteration_id,
398 : const gsl::not_null<blaze::CompressedMatrix<value_type>*> matrix) {
399 : *stored_total_num_points = captured_total_num_points;
400 : *stored_local_first_index = captured_local_first_index;
401 : *iteration_id = 0;
402 : matrix->clear();
403 : matrix->resize(local_num_points, captured_total_num_points);
404 : },
405 : make_not_null(&box));
406 : // Proceed with algorithm
407 : Parallel::get_parallel_component<ParallelComponent>(cache)[element_id]
408 : .perform_algorithm(true);
409 : }
410 : };
411 :
412 : /// Set the operand to a unit vector with a '1' at the current grid point.
413 : /// Applying the operator to this operand gives a column of the matrix. We jump
414 : /// back to this until we have iterated over all grid points.
415 : template <typename BuildMatrixMetavars>
416 1 : struct SetUnitVector {
417 : private:
418 0 : using IterationIdTag = typename BuildMatrixMetavars::iteration_id_tag;
419 0 : using OperandTag = typename BuildMatrixMetavars::operand_tag;
420 :
421 : public:
422 : template <typename DbTags, typename... InboxTags, typename Metavariables,
423 : typename ArrayIndex, typename ActionList,
424 : typename ParallelComponent>
425 0 : static Parallel::iterable_action_return_t apply(
426 : db::DataBox<DbTags>& box,
427 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
428 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
429 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
430 : const ParallelComponent* const /*meta*/) {
431 : const std::optional<size_t> local_unit_vector_index =
432 : detail::local_unit_vector_index(get<IterationIdTag>(box),
433 : get<Tags::LocalFirstIndex>(box),
434 : get<OperandTag>(box).size());
435 : if (local_unit_vector_index.has_value()) {
436 : db::mutate<OperandTag>(
437 : [&local_unit_vector_index](const auto operand) {
438 : operand->data()[*local_unit_vector_index] = 1.;
439 : },
440 : make_not_null(&box));
441 : }
442 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
443 : }
444 : };
445 :
446 : // --- Linear operator will be applied to the unit vector here ---
447 :
448 : /// Write result out to disk, reset operand back to zero, and keep iterating
449 : template <typename BuildMatrixMetavars>
450 1 : struct StoreMatrixColumn {
451 : private:
452 0 : using IterationIdTag = typename BuildMatrixMetavars::iteration_id_tag;
453 0 : using OperandTag = typename BuildMatrixMetavars::operand_tag;
454 0 : using OperatorAppliedToOperandTag =
455 : typename BuildMatrixMetavars::operator_applied_to_operand_tag;
456 0 : using CoordsTag = typename BuildMatrixMetavars::coords_tag;
457 0 : using ArraySectionIdTag = typename BuildMatrixMetavars::array_section_id_tag;
458 0 : using value_type = typename BuildMatrixMetavars::value_type;
459 :
460 : public:
461 0 : using const_global_cache_tags = tmpl::list<Tags::MatrixSubfileName>;
462 :
463 : template <typename DbTags, typename... InboxTags, typename Metavariables,
464 : size_t Dim, typename ActionList, typename ParallelComponent>
465 0 : static Parallel::iterable_action_return_t apply(
466 : db::DataBox<DbTags>& box,
467 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
468 : Parallel::GlobalCache<Metavariables>& cache,
469 : const ElementId<Dim>& element_id, const ActionList /*meta*/,
470 : const ParallelComponent* const /*meta*/) {
471 : const size_t iteration_id = get<IterationIdTag>(box);
472 : const size_t local_first_index = get<Tags::LocalFirstIndex>(box);
473 : if (get<logging::Tags::Verbosity<OptionTags::BuildMatrixOptionsGroup>>(
474 : box) >= Verbosity::Verbose and
475 : local_first_index == 0 and (iteration_id + 1) % 100 == 0) {
476 : Parallel::printf("Column %zu / %zu done.\n", iteration_id + 1,
477 : get<Tags::TotalNumPoints>(box));
478 : }
479 : // This is the result of applying the linear operator to the unit vector. It
480 : // is a column of the operator matrix.
481 : const auto& operator_applied_to_operand =
482 : get<OperatorAppliedToOperandTag>(box);
483 : // Store in sparse matrix
484 : db::mutate<Tags::Matrix<value_type>>(
485 : [&operator_applied_to_operand, &iteration_id](
486 : const gsl::not_null<blaze::CompressedMatrix<value_type>*> matrix) {
487 : for (size_t i = 0; i < operator_applied_to_operand.size(); ++i) {
488 : const auto value = operator_applied_to_operand.data()[i];
489 : if (not equal_within_roundoff(value, 0.)) {
490 : matrix->insert(i, iteration_id, value);
491 : }
492 : }
493 : },
494 : make_not_null(&box));
495 : // Write it out to disk
496 : if (const auto& subfile_name = get<Tags::MatrixSubfileName>(box);
497 : subfile_name.has_value()) {
498 : detail::observe_matrix_column<ParallelComponent>(
499 : iteration_id, operator_applied_to_operand, element_id,
500 : get<domain::Tags::Mesh<Dim>>(box), get<CoordsTag>(box), *subfile_name,
501 : *observers::get_section_observation_key<ArraySectionIdTag>(box),
502 : cache);
503 : }
504 : // Reset operand to zero
505 : const std::optional<size_t> local_unit_vector_index =
506 : detail::local_unit_vector_index(iteration_id, local_first_index,
507 : get<OperandTag>(box).size());
508 : if (local_unit_vector_index.has_value()) {
509 : db::mutate<OperandTag>(
510 : [&local_unit_vector_index](const auto operand) {
511 : operand->data()[*local_unit_vector_index] = 0.;
512 : },
513 : make_not_null(&box));
514 : }
515 : // Keep iterating
516 : db::mutate<IterationIdTag>(
517 : [](const auto local_iteration_id) { ++(*local_iteration_id); },
518 : make_not_null(&box));
519 : if (get<IterationIdTag>(box) < get<Tags::TotalNumPoints>(box)) {
520 : constexpr size_t set_unit_vector_index =
521 : tmpl::index_of<ActionList, SetUnitVector<BuildMatrixMetavars>>::value;
522 : return {Parallel::AlgorithmExecution::Continue, set_unit_vector_index};
523 : }
524 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
525 : }
526 : };
527 :
528 : template <typename Metavariables, typename BuildMatrixMetavars>
529 0 : struct BuildMatrixSingleton {
530 0 : using chare_type = Parallel::Algorithms::Singleton;
531 0 : static constexpr bool checkpoint_data = true;
532 0 : using const_global_cache_tags = tmpl::list<>;
533 0 : using metavariables = Metavariables;
534 0 : using phase_dependent_action_list = tmpl::list<
535 : Parallel::PhaseActions<Parallel::Phase::Initialization, tmpl::list<>>>;
536 0 : using simple_tags_from_options = Parallel::get_simple_tags_from_options<
537 : Parallel::get_initialization_actions_list<phase_dependent_action_list>>;
538 :
539 0 : static void execute_next_phase(
540 : const Parallel::Phase next_phase,
541 : Parallel::CProxy_GlobalCache<Metavariables>& global_cache) {
542 : auto& local_cache = *Parallel::local_branch(global_cache);
543 : Parallel::get_parallel_component<BuildMatrixSingleton>(local_cache)
544 : .start_phase(next_phase);
545 : }
546 : };
547 :
548 : template <typename ElementArrayComponent, typename BuildMatrixMetavars>
549 : struct InvertMatrix;
550 :
551 : template <typename BuildMatrixMetavars>
552 0 : struct AssembleFullMatrix {
553 : private:
554 0 : using value_type = typename BuildMatrixMetavars::value_type;
555 0 : using ArraySectionIdTag = typename BuildMatrixMetavars::array_section_id_tag;
556 0 : using FixedSourcesTag = typename BuildMatrixMetavars::fixed_sources_tag;
557 0 : using ReductionType = std::pair<blaze::CompressedMatrix<value_type>,
558 : typename FixedSourcesTag::type>;
559 :
560 : public:
561 0 : using const_global_cache_tags = tmpl::list<Tags::EnableDirectSolve>;
562 :
563 : template <typename DbTags, typename... InboxTags, typename Metavariables,
564 : size_t Dim, typename ActionList, typename ParallelComponent>
565 0 : static Parallel::iterable_action_return_t apply(
566 : db::DataBox<DbTags>& box,
567 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
568 : Parallel::GlobalCache<Metavariables>& cache,
569 : const ElementId<Dim>& element_id, const ActionList /*meta*/,
570 : const ParallelComponent* const /*meta*/) {
571 : if (not db::get<Tags::EnableDirectSolve>(box)) {
572 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
573 : }
574 : auto& section = Parallel::get_section<ParallelComponent, ArraySectionIdTag>(
575 : make_not_null(&box));
576 : Parallel::contribute_to_reduction<
577 : InvertMatrix<ParallelComponent, BuildMatrixMetavars>>(
578 : Parallel::ReductionData<
579 : Parallel::ReductionDatum<size_t, funcl::AssertEqual<>>,
580 : Parallel::ReductionDatum<std::map<ElementId<Dim>, ReductionType>,
581 : funcl::Merge<>>>{
582 : element_id.grid_index(),
583 : std::map<ElementId<Dim>, ReductionType>{std::make_pair(
584 : element_id, ReductionType{get<Tags::Matrix<value_type>>(box),
585 : get<FixedSourcesTag>(box)})}},
586 : Parallel::get_parallel_component<ParallelComponent>(cache)[element_id],
587 : Parallel::get_parallel_component<
588 : BuildMatrixSingleton<Metavariables, BuildMatrixMetavars>>(cache),
589 : make_not_null(§ion));
590 : // Pause the algorithm for now. The reduction will be received by the next
591 : // action which is responsible for restarting the algorithm.
592 : return {Parallel::AlgorithmExecution::Pause, std::nullopt};
593 : }
594 : };
595 :
596 : template <typename BuildMatrixMetavars>
597 : struct StoreSolution;
598 :
599 : template <typename ElementArrayComponent, typename BuildMatrixMetavars>
600 0 : struct InvertMatrix {
601 : private:
602 0 : using value_type = typename BuildMatrixMetavars::value_type;
603 0 : using FixedSourcesTag = typename BuildMatrixMetavars::fixed_sources_tag;
604 0 : using ReductionType = std::pair<blaze::CompressedMatrix<value_type>,
605 : typename FixedSourcesTag::type>;
606 :
607 : public:
608 : template <typename ParallelComponent, typename DbTagsList,
609 : typename Metavariables, typename ArrayIndex, size_t Dim>
610 0 : static void apply(db::DataBox<DbTagsList>& box,
611 : Parallel::GlobalCache<Metavariables>& cache,
612 : const ArrayIndex& /*array_index*/, const size_t grid_index,
613 : const std::map<ElementId<Dim>, ReductionType>&
614 : matrix_and_sources_slices) {
615 : const size_t total_num_points =
616 : matrix_and_sources_slices.begin()->second.first.columns();
617 : blaze::DynamicMatrix<value_type> matrix(total_num_points, total_num_points);
618 : blaze::DynamicVector<value_type> source(total_num_points);
619 : // Assemble the full matrix and source
620 : size_t i = 0;
621 : for (const auto& [element_id, slice] : matrix_and_sources_slices) {
622 : const auto& [matrix_slice, source_slice] = slice;
623 : ASSERT(matrix_slice.rows() == source_slice.size(),
624 : "Matrix and source slice have different sizes.");
625 : const size_t slice_length = matrix_slice.rows();
626 : blaze::submatrix(matrix, i, 0, slice_length, matrix_slice.columns()) =
627 : matrix_slice;
628 : std::copy_n(source_slice.data(), slice_length, source.data() + i);
629 : i += slice_length;
630 : }
631 : ASSERT(i == total_num_points, "The matrix is not fully assembled. Expected "
632 : << total_num_points << " rows, got "
633 : << i);
634 : // Solve the equations Ax = b. Could use a more efficient solver here if
635 : // needed. But anything iterative could be done with the parallel linear
636 : // solver, the point here is to assemble the full matrix and invert it
637 : // directly.
638 : if (get<logging::Tags::Verbosity<OptionTags::BuildMatrixOptionsGroup>>(
639 : box) >= Verbosity::Quiet) {
640 : Parallel::printf("Inverting matrix of size %zu x %zu\n", total_num_points,
641 : total_num_points);
642 : }
643 : const blaze::DynamicVector<value_type> solution =
644 : blaze::solve(matrix, source);
645 : // Broadcast the solution to the elements
646 : Parallel::simple_action<StoreSolution<BuildMatrixMetavars>>(
647 : Parallel::get_parallel_component<ElementArrayComponent>(cache),
648 : grid_index, solution);
649 : }
650 : };
651 :
652 : template <typename BuildMatrixMetavars>
653 0 : struct StoreSolution {
654 : private:
655 0 : using value_type = typename BuildMatrixMetavars::value_type;
656 0 : using FieldsTag = typename BuildMatrixMetavars::fields_tag;
657 0 : using FixedSourcesTag = typename BuildMatrixMetavars::fixed_sources_tag;
658 :
659 : public:
660 : template <typename ParallelComponent, typename DbTagsList,
661 : typename Metavariables, size_t Dim>
662 0 : static void apply(db::DataBox<DbTagsList>& box,
663 : Parallel::GlobalCache<Metavariables>& cache,
664 : const ElementId<Dim>& element_id, const size_t grid_index,
665 : const blaze::DynamicVector<value_type>& solution) {
666 : if (grid_index != element_id.grid_index()) {
667 : return;
668 : }
669 : const size_t local_first_index = db::get<Tags::LocalFirstIndex>(box);
670 : db::mutate<
671 : FieldsTag,
672 : db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, FieldsTag>>(
673 : [&solution, &local_first_index](const auto fields,
674 : const auto operator_applied_to_fields,
675 : const auto& fixed_sources) {
676 : std::copy_n(
677 : solution.begin() + static_cast<ptrdiff_t>(local_first_index),
678 : fields->size(), fields->data());
679 : // The linear problem is solved now, so Ax = b
680 : *operator_applied_to_fields = fixed_sources;
681 : },
682 : make_not_null(&box), db::get<FixedSourcesTag>(box));
683 : // Proceed with algorithm
684 : Parallel::get_parallel_component<ParallelComponent>(cache)[element_id]
685 : .perform_algorithm(true);
686 : }
687 : };
688 :
689 : template <typename BuildMatrixMetavars>
690 0 : struct ProjectBuildMatrix : tt::ConformsTo<::amr::protocols::Projector> {
691 : private:
692 0 : using IterationIdTag = typename BuildMatrixMetavars::iteration_id_tag;
693 0 : using OperandTag = typename BuildMatrixMetavars::operand_tag;
694 0 : using value_type = typename BuildMatrixMetavars::value_type;
695 :
696 : public:
697 0 : using return_tags =
698 : tmpl::list<Tags::TotalNumPoints, Tags::Matrix<value_type>,
699 : Tags::LocalFirstIndex, IterationIdTag, OperandTag>;
700 0 : using argument_tags = tmpl::list<>;
701 :
702 : template <typename... AmrData>
703 0 : static void apply(
704 : const gsl::not_null<size_t*> total_num_points,
705 : const gsl::not_null<blaze::CompressedMatrix<value_type>*> matrix,
706 : const AmrData&... /*amr_data*/) {
707 : // Reset the built matrix when AMR changes the grid.
708 : // In case AMR is configured to keep coarse grids then the coarse-grid
709 : // elements don't run projectors during AMR, so they are not reset and just
710 : // keep the built matrix. This is good because the coarse grids are
711 : // unchanged and so the matrix is still valid (and can be used as bottom
712 : // solver in multigrid).
713 : *total_num_points = 0;
714 : matrix->clear();
715 : }
716 : };
717 :
718 : /// Dispatch global reduction to get the size of the matrix
719 : template <typename BuildMatrixMetavars>
720 1 : struct ResetBuiltMatrix {
721 : private:
722 0 : using value_type = typename BuildMatrixMetavars::value_type;
723 0 : using ArraySectionIdTag = typename BuildMatrixMetavars::array_section_id_tag;
724 :
725 : public:
726 0 : using const_global_cache_tags = tmpl::list<Tags::SkipResets>;
727 :
728 : template <typename DbTags, typename... InboxTags, typename Metavariables,
729 : size_t Dim, typename ActionList, typename ParallelComponent>
730 0 : static Parallel::iterable_action_return_t apply(
731 : db::DataBox<DbTags>& box,
732 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
733 : Parallel::GlobalCache<Metavariables>& /*cache*/,
734 : const ElementId<Dim>& /*array_index*/, const ActionList /*meta*/,
735 : const ParallelComponent* const /*meta*/) {
736 : // Skip on elements that are not part of the section
737 : if constexpr (not std::is_same_v<ArraySectionIdTag, void>) {
738 : if (not db::get<Parallel::Tags::Section<ParallelComponent,
739 : ArraySectionIdTag>>(box)
740 : .has_value()) {
741 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
742 : }
743 : }
744 : if (db::get<Tags::SkipResets>(box)) {
745 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
746 : }
747 : db::mutate<Tags::TotalNumPoints, Tags::Matrix<value_type>>(
748 : [](const auto total_num_points, const auto matrix) {
749 : *total_num_points = 0;
750 : matrix->clear();
751 : },
752 : make_not_null(&box));
753 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
754 : }
755 : };
756 :
757 : /*!
758 : * \brief Build the explicit matrix representation of the linear operator and
759 : * optionally invert it directly to solve the problem.
760 : *
761 : * This is useful for debugging and analysis only, not to actually solve the
762 : * elliptic problem (that should happen iteratively).
763 : *
764 : * Add the `actions` to the action list to build the matrix. The
765 : * `ApplyOperatorActions` template parameter are the actions that apply the
766 : * linear operator to the `OperandTag`. Also add the `amr_projectors` to the
767 : * list of AMR projectors and the `register_actions`
768 : *
769 : * \tparam FieldsTag The solution will be stored here (if direct solve is
770 : * enabled).
771 : * \tparam FixedSourcesTag The source `b` in the problem `Ax = b`. Only used
772 : * if direct solve is enabled.
773 : * \tparam OperandTag Will be set to unit vectors, with the '1' moving through
774 : * all points over the course of the iteration.
775 : * \tparam OperatorAppliedToOperandTag Where the `ApplyOperatorActions` store
776 : * the result of applying the linear operator to the `OperandTag`.
777 : * \tparam CoordsTag The tag of the coordinates observed alongside the matrix
778 : * for volume data visualization.
779 : * \tparam ArraySectionIdTag Can identify a subset of elements that this
780 : * algorithm should run over, e.g. in a multigrid setting.
781 : */
782 : template <typename FieldsTag, typename FixedSourcesTag, typename OperandTag,
783 : typename OperatorAppliedToOperandTag, typename CoordsTag,
784 : typename ArraySectionIdTag = void>
785 1 : struct BuildMatrix {
786 : private:
787 0 : using BuildMatrixMetavars =
788 : detail::BuildMatrixMetavars<FieldsTag, FixedSourcesTag, OperandTag,
789 : OperatorAppliedToOperandTag, CoordsTag,
790 : ArraySectionIdTag>;
791 : static_assert(
792 : std::is_same_v<FieldsTag, OperandTag>,
793 : "The operand and the fields tags must be the same. This is just so that "
794 : "at the end of the algorithm the operator is applied to the solution "
795 : "fields. This restriction can be lifted if needed by copying the fields "
796 : "into the operand and back.");
797 :
798 : public:
799 : template <typename Metavariables>
800 0 : using component_list =
801 : tmpl::list<BuildMatrixSingleton<Metavariables, BuildMatrixMetavars>>;
802 :
803 : template <typename ApplyOperatorActions>
804 0 : using actions =
805 : tmpl::list<CollectTotalNumPoints<BuildMatrixMetavars>,
806 : // PrepareBuildMatrix is called on reduction broadcast
807 : SetUnitVector<BuildMatrixMetavars>, ApplyOperatorActions,
808 : StoreMatrixColumn<BuildMatrixMetavars>,
809 : // Algorithm iterates until matrix is complete, then proceeds
810 : // below
811 : AssembleFullMatrix<BuildMatrixMetavars>,
812 : // Apply operator to the solution. Note that FieldsTag and
813 : // OperandTag must be the same for this (see assert above).
814 : // Note also that this operator application is not needed in
815 : // most cases, as the linear problem is solved exactly and
816 : // therefore Ax = b is set in 'StoreSolution'. However, this
817 : // operator application is needed if the operator changes
818 : // between nonlinear solver iterations but we skip the reset,
819 : // so the matrix represents the old operator. It's just one
820 : // more operator application, so we keep it always enabled for
821 : // now.
822 : ApplyOperatorActions,
823 : ::Actions::Label<typename BuildMatrixMetavars::end_label>>;
824 :
825 0 : using amr_projectors = tmpl::list<ProjectBuildMatrix<BuildMatrixMetavars>>;
826 :
827 : /// Add to the register phase to enable observations
828 1 : using register_actions = tmpl::list<observers::Actions::RegisterWithObservers<
829 : detail::RegisterWithVolumeObserver<ArraySectionIdTag>>>;
830 :
831 : /// Add to the action list to reset the matrix
832 1 : using reset_actions = tmpl::list<ResetBuiltMatrix<BuildMatrixMetavars>>;
833 : };
834 :
835 : } // namespace Actions
836 : } // namespace LinearSolver
|
