Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <optional>
8 : #include <string>
9 : #include <tuple>
10 : #include <type_traits>
11 : #include <utility>
12 : #include <vector>
13 :
14 : #include "DataStructures/DataBox/DataBox.hpp"
15 : #include "DataStructures/TaggedTuple.hpp"
16 : #include "Domain/Structure/ElementId.hpp"
17 : #include "Domain/Tags.hpp"
18 : #include "IO/H5/TensorData.hpp"
19 : #include "IO/Observer/GetSectionObservationKey.hpp"
20 : #include "IO/Observer/ObservationId.hpp"
21 : #include "IO/Observer/Tags.hpp"
22 : #include "IO/Observer/TypeOfObservation.hpp"
23 : #include "IO/Observer/VolumeActions.hpp"
24 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
25 : #include "Parallel/AlgorithmExecution.hpp"
26 : #include "Parallel/ArrayComponentId.hpp"
27 : #include "Parallel/GlobalCache.hpp"
28 : #include "Parallel/TypeTraits.hpp"
29 : #include "ParallelAlgorithms/LinearSolver/Schwarz/ElementCenteredSubdomainData.hpp"
30 : #include "ParallelAlgorithms/LinearSolver/Schwarz/Tags.hpp"
31 : #include "ParallelAlgorithms/LinearSolver/Tags.hpp"
32 : #include "Utilities/Algorithm.hpp"
33 : #include "Utilities/Gsl.hpp"
34 : #include "Utilities/MakeString.hpp"
35 : #include "Utilities/PrettyType.hpp"
36 : #include "Utilities/TMPL.hpp"
37 :
38 : namespace LinearSolver::Schwarz::detail {
39 :
40 : template <typename OptionsGroup, typename ArraySectionIdTag>
41 : struct RegisterWithVolumeObserver {
42 : template <typename ParallelComponent, typename DbTagsList,
43 : typename ArrayIndex>
44 : static std::pair<observers::TypeOfObservation, observers::ObservationKey>
45 : register_info(const db::DataBox<DbTagsList>& box,
46 : const ArrayIndex& /*array_index*/) {
47 : const std::optional<std::string> section_observation_key =
48 : observers::get_section_observation_key<ArraySectionIdTag>(box);
49 : const std::string subfile_path = "/" + pretty_type::name<OptionsGroup>() +
50 : section_observation_key.value_or("");
51 : return {observers::TypeOfObservation::Volume,
52 : observers::ObservationKey(subfile_path)};
53 : }
54 : };
55 :
56 : // Contribute the volume data recorded in the other actions to the observer at
57 : // the end of a step.
58 : template <typename FieldsTag, typename OptionsGroup, typename SubdomainOperator,
59 : typename ArraySectionIdTag>
60 : struct ObserveVolumeData {
61 : private:
62 : using fields_tag = FieldsTag;
63 : using residual_tag =
64 : db::add_tag_prefix<LinearSolver::Tags::Residual, fields_tag>;
65 : static constexpr size_t Dim = SubdomainOperator::volume_dim;
66 : using SubdomainData =
67 : ElementCenteredSubdomainData<Dim, typename residual_tag::tags_list>;
68 : using volume_data_tag =
69 : Tags::VolumeDataForOutput<SubdomainData, OptionsGroup>;
70 : using VolumeDataVars = typename volume_data_tag::type::ElementData;
71 :
72 : public:
73 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
74 : size_t Dim, typename ActionList, typename ParallelComponent>
75 : static Parallel::iterable_action_return_t apply(
76 : db::DataBox<DbTagsList>& box,
77 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
78 : Parallel::GlobalCache<Metavariables>& cache,
79 : const ElementId<Dim>& element_id, const ActionList /*meta*/,
80 : const ParallelComponent* const /*meta*/) {
81 : if (not db::get<LinearSolver::Tags::OutputVolumeData<OptionsGroup>>(box)) {
82 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
83 : }
84 : const auto& volume_data = db::get<volume_data_tag>(box);
85 : const auto& observation_id =
86 : db::get<LinearSolver::Tags::ObservationId<OptionsGroup>>(box);
87 : const auto& mesh = db::get<domain::Tags::Mesh<Dim>>(box);
88 : const auto& inertial_coords =
89 : db::get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(box);
90 : // Collect tensor components to observe
91 : std::vector<TensorComponent> components{};
92 : const auto record_tensor_components = [&components](
93 : const auto tensor_tag_v,
94 : const auto& tensor,
95 : const std::string& suffix = "") {
96 : using tensor_tag = std::decay_t<decltype(tensor_tag_v)>;
97 : using TensorType = std::decay_t<decltype(tensor)>;
98 : using VectorType = typename TensorType::type;
99 : using ValueType = typename VectorType::value_type;
100 : for (size_t i = 0; i < tensor.size(); ++i) {
101 : const std::string component_name =
102 : db::tag_name<tensor_tag>() + suffix + tensor.component_suffix(i);
103 : if constexpr (std::is_same_v<ValueType, std::complex<double>>) {
104 : components.emplace_back("Re(" + component_name + ")",
105 : real(tensor[i]));
106 : components.emplace_back("Im(" + component_name + ")",
107 : imag(tensor[i]));
108 : } else {
109 : components.emplace_back(component_name, tensor[i]);
110 : }
111 : }
112 : };
113 : record_tensor_components(domain::Tags::Coordinates<Dim, Frame::Inertial>{},
114 : inertial_coords);
115 : const auto& all_intruding_extents =
116 : db::get<Tags::IntrudingExtents<Dim, OptionsGroup>>(box);
117 : const VolumeDataVars zero_vars{mesh.number_of_grid_points(), 0.};
118 : tmpl::for_each<typename VolumeDataVars::tags_list>(
119 : [&volume_data, &record_tensor_components, &mesh, &all_intruding_extents,
120 : &zero_vars, &element_id](auto tag_v) {
121 : using tag = tmpl::type_from<decltype(tag_v)>;
122 : record_tensor_components(tag{}, get<tag>(volume_data.element_data),
123 : "_Center");
124 : for (const auto direction : Direction<Dim>::all_directions()) {
125 : const auto direction_predicate =
126 : [&direction](const auto& overlap_data) {
127 : return overlap_data.first.direction() == direction;
128 : };
129 : const auto num_overlaps =
130 : alg::count_if(volume_data.overlap_data, direction_predicate);
131 : if (num_overlaps == 0) {
132 : // No overlap data for this direction (e.g. external boundary),
133 : // record zero
134 : record_tensor_components(tag{}, get<tag>(zero_vars),
135 : "_Overlap" + get_output(direction));
136 : } else if (num_overlaps == 1) {
137 : // Overlap data from a single neighbor, record it
138 : const auto& overlap_data =
139 : alg::find_if(volume_data.overlap_data, direction_predicate)
140 : ->second;
141 : const auto& intruding_extents =
142 : gsl::at(all_intruding_extents, direction.dimension());
143 : record_tensor_components(
144 : tag{},
145 : get<tag>(extended_overlap_data(overlap_data, mesh.extents(),
146 : intruding_extents, direction)),
147 : "_Overlap" + get_output(direction));
148 : } else {
149 : ERROR("Multiple neighbors ("
150 : << num_overlaps << ") overlap with element " << element_id
151 : << " in direction " << direction
152 : << ", but we can record only one for volume data output.");
153 : }
154 : }
155 : });
156 :
157 : // Contribute tensor components to observer
158 : const std::optional<std::string> section_observation_key =
159 : observers::get_section_observation_key<ArraySectionIdTag>(box);
160 : const std::string subfile_path = "/" + pretty_type::name<OptionsGroup>() +
161 : section_observation_key.value_or("");
162 : observers::contribute_volume_data<
163 : not Parallel::is_nodegroup_v<ParallelComponent>>(
164 : cache, observers::ObservationId(observation_id, subfile_path),
165 : subfile_path,
166 : Parallel::make_array_component_id<ParallelComponent>(element_id),
167 : ElementVolumeData{element_id, std::move(components), mesh});
168 :
169 : // Increment observation ID
170 : db::mutate<LinearSolver::Tags::ObservationId<OptionsGroup>>(
171 : [](const auto local_observation_id) { ++(*local_observation_id); },
172 : make_not_null(&box));
173 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
174 : }
175 : };
176 :
177 : } // namespace LinearSolver::Schwarz::detail
|
