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 <cstddef>
12 : #include <map>
13 : #include <optional>
14 : #include <utility>
15 : #include <vector>
16 :
17 : #include "DataStructures/DataBox/DataBox.hpp"
18 : #include "DataStructures/DataBox/Tag.hpp"
19 : #include "Domain/Structure/ElementId.hpp"
20 : #include "Domain/Tags.hpp"
21 : #include "IO/H5/TensorData.hpp"
22 : #include "IO/Logging/Tags.hpp"
23 : #include "IO/Logging/Verbosity.hpp"
24 : #include "IO/Observer/Actions/RegisterWithObservers.hpp"
25 : #include "IO/Observer/GetSectionObservationKey.hpp"
26 : #include "IO/Observer/ObserverComponent.hpp"
27 : #include "IO/Observer/VolumeActions.hpp"
28 : #include "Parallel/AlgorithmExecution.hpp"
29 : #include "Parallel/GetSection.hpp"
30 : #include "Parallel/Phase.hpp"
31 : #include "Parallel/Printf/Printf.hpp"
32 : #include "Parallel/Reduction.hpp"
33 : #include "Parallel/Section.hpp"
34 : #include "Parallel/Tags/Section.hpp"
35 : #include "ParallelAlgorithms/Amr/Protocols/Projector.hpp"
36 : #include "Utilities/Functional.hpp"
37 : #include "Utilities/GetOutput.hpp"
38 : #include "Utilities/ProtocolHelpers.hpp"
39 : #include "Utilities/TMPL.hpp"
40 : #include "Utilities/TaggedTuple.hpp"
41 :
42 : namespace LinearSolver {
43 0 : namespace OptionTags {
44 :
45 0 : struct BuildMatrixOptionsGroup {
46 0 : static std::string name() { return "BuildMatrix"; }
47 0 : static constexpr Options::String help = {
48 : "Options for building the explicit matrix representation of the linear "
49 : "operator. This is done by applying the linear operator to unit "
50 : "vectors and is useful for debugging and analysis only, not to actually "
51 : "solve the elliptic problem (that should happen iteratively)."};
52 : };
53 :
54 0 : struct MatrixSubfileName {
55 0 : using type = std::string;
56 0 : using group = BuildMatrixOptionsGroup;
57 0 : static constexpr Options::String help = {
58 : "Subfile name in the volume data H5 files where the matrix will be "
59 : "stored. Each observation in the subfile is a column of the matrix. The "
60 : "row index is the order of elements defined by the ElementId in the "
61 : "volume data, by the order of tensor components encoded in the name of "
62 : "the components, and by the contiguous ordering of grid points for each "
63 : "component."};
64 : };
65 :
66 : } // namespace OptionTags
67 :
68 1 : namespace Tags {
69 :
70 : /// Subfile name in the volume data H5 files where the matrix will be stored.
71 1 : struct MatrixSubfileName : db::SimpleTag {
72 0 : using type = std::string;
73 0 : using option_tags = tmpl::list<OptionTags::MatrixSubfileName>;
74 0 : static constexpr bool pass_metavariables = false;
75 0 : static type create_from_options(const type& value) { return value; }
76 : };
77 :
78 : /// Size of the matrix: number of grid points times number of variables
79 1 : struct TotalNumPoints : db::SimpleTag {
80 0 : using type = size_t;
81 : };
82 :
83 : /// Index of the first point in this element
84 1 : struct LocalFirstIndex : db::SimpleTag {
85 0 : using type = size_t;
86 : };
87 :
88 : } // namespace Tags
89 :
90 1 : namespace Actions {
91 :
92 : namespace detail {
93 :
94 : /// \brief The total number of grid points (size of the matrix) and the index of
95 : /// the first grid point in this element (the offset into the matrix
96 : /// corresponding to this element).
97 : template <size_t Dim>
98 : std::pair<size_t, size_t> total_num_points_and_local_first_index(
99 : const ElementId<Dim>& element_id,
100 : const std::map<ElementId<Dim>, size_t>& num_points_per_element,
101 : size_t num_vars);
102 :
103 : /// \brief The index of the '1' of the unit vector in this element, or
104 : /// std::nullopt if the '1' is in another element.
105 : ///
106 : /// \param iteration_id enumerates all grid points
107 : /// \param local_first_index the index of the first grid point in this element
108 : /// \param local_num_points the number of grid points in this element
109 : std::optional<size_t> local_unit_vector_index(size_t iteration_id,
110 : size_t local_first_index,
111 : size_t local_num_points);
112 :
113 : /// \brief Observe matrix column as volume data.
114 : ///
115 : /// This is the format of the volume data:
116 : ///
117 : /// - Columns of the matrix are enumerated by the observation ID
118 : /// - Rows are enumerated by the order of elements defined by the ElementId
119 : /// in the volume data, by the order of tensor components encoded in the
120 : /// name of the components, and by the contiguous ordering of grid points
121 : /// for each component.
122 : /// - We also observe coordinates so the data can be plotted as volume data.
123 : ///
124 : /// The matrix can be reconstructed from this data in Python with the function
125 : /// `spectre.Elliptic.ReadH5.read_matrix`.
126 : template <typename ParallelComponent, typename OperatorAppliedToOperandTags,
127 : size_t Dim, typename CoordsFrame, typename Metavariables>
128 : void observe_matrix_column(
129 : const size_t column,
130 : const Variables<OperatorAppliedToOperandTags>& operator_applied_to_operand,
131 : const ElementId<Dim>& element_id, const Mesh<Dim>& mesh,
132 : const tnsr::I<DataVector, Dim, CoordsFrame>& coords,
133 : const std::string& subfile_name, const std::string& section_observation_key,
134 : Parallel::GlobalCache<Metavariables>& cache) {
135 : std::vector<TensorComponent> observe_components{};
136 : observe_components.reserve(
137 : Dim +
138 : Variables<
139 : OperatorAppliedToOperandTags>::number_of_independent_components);
140 : for (size_t i = 0; i < Dim; ++i) {
141 : observe_components.emplace_back(
142 : get_output(CoordsFrame{}) + "Coordinates" + coords.component_suffix(i),
143 : coords.get(i));
144 : }
145 : size_t component_i = 0;
146 : tmpl::for_each<OperatorAppliedToOperandTags>([&operator_applied_to_operand,
147 : &observe_components,
148 : &component_i](auto tag_v) {
149 : using tag = tmpl::type_from<decltype(tag_v)>;
150 : const auto& tensor = get<tag>(operator_applied_to_operand);
151 : using TensorType = std::decay_t<decltype(tensor)>;
152 : using VectorType = typename TensorType::type;
153 : using ValueType = typename VectorType::value_type;
154 : for (size_t i = 0; i < tensor.size(); ++i) {
155 : if constexpr (std::is_same_v<ValueType, std::complex<double>>) {
156 : observe_components.emplace_back(
157 : "Re(Variable_" + std::to_string(component_i) + ")",
158 : real(tensor[i]));
159 : observe_components.emplace_back(
160 : "Im(Variable_" + std::to_string(component_i) + ")",
161 : imag(tensor[i]));
162 : } else {
163 : observe_components.emplace_back(
164 : "Variable_" + std::to_string(component_i), tensor[i]);
165 : }
166 : ++component_i;
167 : }
168 : });
169 : auto& local_observer = *Parallel::local_branch(
170 : Parallel::get_parallel_component<observers::Observer<Metavariables>>(
171 : cache));
172 : Parallel::simple_action<observers::Actions::ContributeVolumeData>(
173 : local_observer,
174 : observers::ObservationId(static_cast<double>(column),
175 : subfile_name + section_observation_key),
176 : "/" + subfile_name,
177 : Parallel::make_array_component_id<ParallelComponent>(element_id),
178 : ElementVolumeData{element_id, std::move(observe_components), mesh});
179 : }
180 :
181 : /// \brief Register with the volume observer
182 : template <typename ArraySectionIdTag>
183 : struct RegisterWithVolumeObserver {
184 : template <typename ParallelComponent, typename DbTagsList,
185 : typename ArrayIndex>
186 : static std::pair<observers::TypeOfObservation, observers::ObservationKey>
187 : register_info(const db::DataBox<DbTagsList>& box,
188 : const ArrayIndex& /*array_index*/) {
189 : // Get the observation key, or "Unused" if the element does not belong
190 : // to a section with this tag. In the latter case, no observations will
191 : // ever be contributed.
192 : const std::optional<std::string> section_observation_key =
193 : observers::get_section_observation_key<ArraySectionIdTag>(box);
194 : ASSERT(section_observation_key != "Unused",
195 : "The identifier 'Unused' is reserved to indicate that no "
196 : "observations with this key will be contributed. Use a different "
197 : "key, or change the identifier 'Unused' to something else.");
198 : return {
199 : observers::TypeOfObservation::Volume,
200 : observers::ObservationKey(get<Tags::MatrixSubfileName>(box) +
201 : section_observation_key.value_or("Unused"))};
202 : }
203 : };
204 :
205 : } // namespace detail
206 :
207 : /// \cond
208 : template <typename IterationIdTag, typename OperandTag,
209 : typename OperatorAppliedToOperandTag, typename CoordsTag,
210 : typename ArraySectionIdTag>
211 : struct PrepareBuildMatrix;
212 : template <typename IterationIdTag, typename OperandTag,
213 : typename OperatorAppliedToOperandTag, typename CoordsTag,
214 : typename ArraySectionIdTag>
215 : struct StoreMatrixColumn;
216 : /// \endcond
217 :
218 : /// Dispatch global reduction to get the size of the matrix
219 : template <typename IterationIdTag, typename OperandTag,
220 : typename OperatorAppliedToOperandTag, typename CoordsTag,
221 : typename ArraySectionIdTag>
222 1 : struct CollectTotalNumPoints {
223 0 : using simple_tags = tmpl::list<Tags::TotalNumPoints, Tags::LocalFirstIndex,
224 : IterationIdTag, OperandTag>;
225 0 : using compute_tags = tmpl::list<>;
226 0 : using const_global_cache_tags =
227 : tmpl::list<logging::Tags::Verbosity<OptionTags::BuildMatrixOptionsGroup>>;
228 :
229 : template <typename DbTags, typename... InboxTags, typename Metavariables,
230 : size_t Dim, typename ActionList, typename ParallelComponent>
231 0 : static Parallel::iterable_action_return_t apply(
232 : db::DataBox<DbTags>& box,
233 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
234 : Parallel::GlobalCache<Metavariables>& cache,
235 : const ElementId<Dim>& array_index, const ActionList /*meta*/,
236 : const ParallelComponent* const /*meta*/) {
237 : // Skip everything on elements that are not part of the section
238 : if constexpr (not std::is_same_v<ArraySectionIdTag, void>) {
239 : if (not db::get<Parallel::Tags::Section<ParallelComponent,
240 : ArraySectionIdTag>>(box)
241 : .has_value()) {
242 : constexpr size_t last_action_index = tmpl::index_of<
243 : ActionList, StoreMatrixColumn<IterationIdTag, OperandTag,
244 : OperatorAppliedToOperandTag,
245 : CoordsTag, ArraySectionIdTag>>::value;
246 : return {Parallel::AlgorithmExecution::Continue, last_action_index + 1};
247 : }
248 : }
249 : db::mutate<Tags::TotalNumPoints, OperandTag>(
250 : [](const auto total_num_points, const auto operand,
251 : const size_t num_points) {
252 : // Set num points to zero until we know the total number of points
253 : *total_num_points = 0;
254 : // Size operand and fill with zeros
255 : operand->initialize(num_points, 0.);
256 : },
257 : make_not_null(&box),
258 : get<domain::Tags::Mesh<Dim>>(box).number_of_grid_points());
259 : // Collect the total number of grid points
260 : auto& section = Parallel::get_section<ParallelComponent, ArraySectionIdTag>(
261 : make_not_null(&box));
262 : Parallel::contribute_to_reduction<PrepareBuildMatrix<
263 : IterationIdTag, OperandTag, OperatorAppliedToOperandTag, CoordsTag,
264 : ArraySectionIdTag>>(
265 : Parallel::ReductionData<Parallel::ReductionDatum<
266 : std::map<ElementId<Dim>, size_t>, funcl::Merge<>>>{
267 : std::map<ElementId<Dim>, size_t>{
268 : std::make_pair(array_index, get<OperandTag>(box).size())}},
269 : Parallel::get_parallel_component<ParallelComponent>(cache)[array_index],
270 : Parallel::get_parallel_component<ParallelComponent>(cache),
271 : make_not_null(§ion));
272 : // Pause the algorithm for now. The reduction will be broadcast to the next
273 : // action which is responsible for restarting the algorithm.
274 : return {Parallel::AlgorithmExecution::Pause, std::nullopt};
275 : }
276 : };
277 :
278 : /// Receive the reduction and initialize the algorithm
279 : template <typename IterationIdTag, typename OperandTag,
280 : typename OperatorAppliedToOperandTag, typename CoordsTag,
281 : typename ArraySectionIdTag>
282 1 : struct PrepareBuildMatrix {
283 : public:
284 : template <typename ParallelComponent, typename DbTagsList,
285 : typename Metavariables, size_t Dim>
286 0 : static void apply(
287 : db::DataBox<DbTagsList>& box, Parallel::GlobalCache<Metavariables>& cache,
288 : const ElementId<Dim>& element_id,
289 : const std::map<ElementId<Dim>, size_t>& num_points_per_element) {
290 : const auto [total_num_points, local_first_index] =
291 : detail::total_num_points_and_local_first_index(
292 : element_id, num_points_per_element,
293 : OperandTag::type::number_of_independent_components);
294 : if (get<logging::Tags::Verbosity<OptionTags::BuildMatrixOptionsGroup>>(
295 : box) >= Verbosity::Quiet and
296 : local_first_index == 0) {
297 : Parallel::printf(
298 : "Building explicit matrix representation of size %zu x %zu.\n",
299 : total_num_points, total_num_points);
300 : }
301 : db::mutate<Tags::TotalNumPoints, Tags::LocalFirstIndex, IterationIdTag>(
302 : [captured_total_num_points = total_num_points,
303 : captured_local_first_index = local_first_index](
304 : const auto stored_total_num_points,
305 : const auto stored_local_first_index, const auto iteration_id) {
306 : *stored_total_num_points = captured_total_num_points;
307 : *stored_local_first_index = captured_local_first_index;
308 : *iteration_id = 0;
309 : },
310 : make_not_null(&box));
311 : // Proceed with algorithm
312 : Parallel::get_parallel_component<ParallelComponent>(cache)[element_id]
313 : .perform_algorithm(true);
314 : }
315 : };
316 :
317 : /// Set the operand to a unit vector with a '1' at the current grid point.
318 : /// Applying the operator to this operand gives a column of the matrix. We jump
319 : /// back to this until we have iterated over all grid points.
320 : template <typename IterationIdTag, typename OperandTag,
321 : typename OperatorAppliedToOperandTag, typename CoordsTag,
322 : typename ArraySectionIdTag>
323 1 : struct SetUnitVector {
324 : template <typename DbTags, typename... InboxTags, typename Metavariables,
325 : typename ArrayIndex, typename ActionList,
326 : typename ParallelComponent>
327 0 : static Parallel::iterable_action_return_t apply(
328 : db::DataBox<DbTags>& box,
329 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
330 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
331 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
332 : const ParallelComponent* const /*meta*/) {
333 : const std::optional<size_t> local_unit_vector_index =
334 : detail::local_unit_vector_index(get<IterationIdTag>(box),
335 : get<Tags::LocalFirstIndex>(box),
336 : get<OperandTag>(box).size());
337 : if (local_unit_vector_index.has_value()) {
338 : db::mutate<OperandTag>(
339 : [&local_unit_vector_index](const auto operand) {
340 : operand->data()[*local_unit_vector_index] = 1.;
341 : },
342 : make_not_null(&box));
343 : }
344 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
345 : }
346 : };
347 :
348 : // --- Linear operator will be applied to the unit vector here ---
349 :
350 : /// Write result out to disk, reset operand back to zero, and keep iterating
351 : template <typename IterationIdTag, typename OperandTag,
352 : typename OperatorAppliedToOperandTag, typename CoordsTag,
353 : typename ArraySectionIdTag>
354 1 : struct StoreMatrixColumn {
355 0 : using const_global_cache_tags = tmpl::list<Tags::MatrixSubfileName>;
356 :
357 : template <typename DbTags, typename... InboxTags, typename Metavariables,
358 : size_t Dim, typename ActionList, typename ParallelComponent>
359 0 : static Parallel::iterable_action_return_t apply(
360 : db::DataBox<DbTags>& box,
361 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
362 : Parallel::GlobalCache<Metavariables>& cache,
363 : const ElementId<Dim>& element_id, const ActionList /*meta*/,
364 : const ParallelComponent* const /*meta*/) {
365 : const size_t iteration_id = get<IterationIdTag>(box);
366 : const size_t local_first_index = get<Tags::LocalFirstIndex>(box);
367 : if (get<logging::Tags::Verbosity<OptionTags::BuildMatrixOptionsGroup>>(
368 : box) >= Verbosity::Verbose and
369 : local_first_index == 0) {
370 : Parallel::printf("Column %zu / %zu done.\n", iteration_id + 1,
371 : get<Tags::TotalNumPoints>(box));
372 : }
373 : // This is the result of applying the linear operator to the unit vector. It
374 : // is a column of the operator matrix.
375 : const auto& operator_applied_to_operand =
376 : get<OperatorAppliedToOperandTag>(box);
377 : // Write it out to disk
378 : detail::observe_matrix_column<ParallelComponent>(
379 : iteration_id, operator_applied_to_operand, element_id,
380 : get<domain::Tags::Mesh<Dim>>(box), get<CoordsTag>(box),
381 : get<Tags::MatrixSubfileName>(box),
382 : *observers::get_section_observation_key<ArraySectionIdTag>(box), cache);
383 : // Reset operand to zero
384 : const std::optional<size_t> local_unit_vector_index =
385 : detail::local_unit_vector_index(iteration_id, local_first_index,
386 : get<OperandTag>(box).size());
387 : if (local_unit_vector_index.has_value()) {
388 : db::mutate<OperandTag>(
389 : [&local_unit_vector_index](const auto operand) {
390 : operand->data()[*local_unit_vector_index] = 0.;
391 : },
392 : make_not_null(&box));
393 : }
394 : // Keep iterating
395 : db::mutate<IterationIdTag>(
396 : [](const auto local_iteration_id) { ++(*local_iteration_id); },
397 : make_not_null(&box));
398 : if (get<IterationIdTag>(box) < get<Tags::TotalNumPoints>(box)) {
399 : constexpr size_t set_unit_vector_index = tmpl::index_of<
400 : ActionList,
401 : SetUnitVector<IterationIdTag, OperandTag, OperatorAppliedToOperandTag,
402 : CoordsTag, ArraySectionIdTag>>::value;
403 : return {Parallel::AlgorithmExecution::Continue, set_unit_vector_index};
404 : }
405 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
406 : }
407 : };
408 :
409 : template <typename IterationIdTag, typename OperandTag,
410 : typename OperatorAppliedToOperandTag, typename CoordsTag,
411 : typename ArraySectionIdTag>
412 0 : struct ProjectBuildMatrix : tt::ConformsTo<::amr::protocols::Projector> {
413 0 : using return_tags = tmpl::list<Tags::TotalNumPoints, Tags::LocalFirstIndex,
414 : IterationIdTag, OperandTag>;
415 0 : using argument_tags = tmpl::list<>;
416 :
417 : template <typename... AmrData>
418 0 : static void apply(const gsl::not_null<size_t*> /*unused*/,
419 : const AmrData&... /*amr_data*/) {
420 : // Nothing to do. Everything gets initialized at the start of the algorithm.
421 : }
422 : };
423 :
424 : /*!
425 : * \brief Build the explicit matrix representation of the linear operator.
426 : *
427 : * This is useful for debugging and analysis only, not to actually solve the
428 : * elliptic problem (that should happen iteratively).
429 : *
430 : * Add the `actions` to the action list to build the matrix. The
431 : * `ApplyOperatorActions` template parameter are the actions that apply the
432 : * linear operator to the `OperandTag`. Also add the `amr_projectors` to the
433 : * list of AMR projectors and the `register_actions`
434 : *
435 : * \tparam IterationIdTag Used to keep track of the iteration over all matrix
436 : * columns. Should be the same that's used to identify iterations in the
437 : * `ApplyOperatorActions`.
438 : * \tparam OperandTag Will be set to unit vectors, with the '1' moving through
439 : * all points over the course of the iteration.
440 : * \tparam OperatorAppliedToOperandTag Where the `ApplyOperatorActions` store
441 : * the result of applying the linear operator to the `OperandTag`.
442 : * \tparam CoordsTag The tag of the coordinates observed alongside the matrix
443 : * for volume data visualization.
444 : * \tparam ArraySectionIdTag Can identify a subset of elements that this
445 : * algorithm should run over, e.g. in a multigrid setting.
446 : */
447 : template <typename IterationIdTag, typename OperandTag,
448 : typename OperatorAppliedToOperandTag, typename CoordsTag,
449 : typename ArraySectionIdTag = void>
450 1 : struct BuildMatrix {
451 : template <typename ApplyOperatorActions>
452 0 : using actions = tmpl::list<
453 : CollectTotalNumPoints<IterationIdTag, OperandTag,
454 : OperatorAppliedToOperandTag, CoordsTag,
455 : ArraySectionIdTag>,
456 : // PrepareBuildMatrix is called on reduction broadcast
457 : SetUnitVector<IterationIdTag, OperandTag, OperatorAppliedToOperandTag,
458 : CoordsTag, ArraySectionIdTag>,
459 : ApplyOperatorActions,
460 : StoreMatrixColumn<IterationIdTag, OperandTag, OperatorAppliedToOperandTag,
461 : CoordsTag, ArraySectionIdTag>>;
462 :
463 0 : using amr_projectors =
464 : tmpl::list<ProjectBuildMatrix<IterationIdTag, OperandTag,
465 : OperatorAppliedToOperandTag, CoordsTag,
466 : ArraySectionIdTag>>;
467 :
468 : /// Add to the register phase to enable observations
469 1 : using register_actions = tmpl::list<observers::Actions::RegisterWithObservers<
470 : detail::RegisterWithVolumeObserver<ArraySectionIdTag>>>;
471 : };
472 :
473 : } // namespace Actions
474 : } // namespace LinearSolver
|
