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