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 <functional>
8 : #include <initializer_list>
9 : #include <optional>
10 : #include <pup.h>
11 : #include <string>
12 : #include <type_traits>
13 : #include <unordered_map>
14 : #include <utility>
15 : #include <vector>
16 :
17 : #include "DataStructures/DataBox/ObservationBox.hpp"
18 : #include "DataStructures/DataBox/Prefixes.hpp"
19 : #include "DataStructures/DataBox/TagName.hpp"
20 : #include "DataStructures/DataBox/ValidateSelection.hpp"
21 : #include "DataStructures/DataVector.hpp"
22 : #include "DataStructures/FloatingPointType.hpp"
23 : #include "Domain/Structure/ElementId.hpp"
24 : #include "Domain/Tags.hpp"
25 : #include "IO/H5/TensorData.hpp"
26 : #include "IO/Observer/GetSectionObservationKey.hpp"
27 : #include "IO/Observer/ObservationId.hpp"
28 : #include "IO/Observer/ObserverComponent.hpp"
29 : #include "IO/Observer/Tags.hpp"
30 : #include "IO/Observer/VolumeActions.hpp"
31 : #include "NumericalAlgorithms/Interpolation/RegularGridInterpolant.hpp"
32 : #include "Options/Auto.hpp"
33 : #include "Options/String.hpp"
34 : #include "Parallel/ArrayComponentId.hpp"
35 : #include "Parallel/ArrayIndex.hpp"
36 : #include "Parallel/GlobalCache.hpp"
37 : #include "Parallel/Invoke.hpp"
38 : #include "Parallel/Local.hpp"
39 : #include "Parallel/Printf/Printf.hpp"
40 : #include "ParallelAlgorithms/Events/Tags.hpp"
41 : #include "ParallelAlgorithms/EventsAndTriggers/Event.hpp"
42 : #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
43 : #include "Utilities/Algorithm.hpp"
44 : #include "Utilities/ErrorHandling/Assert.hpp"
45 : #include "Utilities/Literals.hpp"
46 : #include "Utilities/MakeString.hpp"
47 : #include "Utilities/Numeric.hpp"
48 : #include "Utilities/OptionalHelpers.hpp"
49 : #include "Utilities/Serialization/CharmPupable.hpp"
50 : #include "Utilities/Serialization/PupStlCpp17.hpp"
51 : #include "Utilities/StdHelpers.hpp"
52 : #include "Utilities/TMPL.hpp"
53 : #include "Utilities/TypeTraits/IsA.hpp"
54 :
55 : /// \cond
56 : template <size_t Dim>
57 : class Mesh;
58 : namespace Frame {
59 : struct Inertial;
60 : } // namespace Frame
61 : /// \endcond
62 :
63 : namespace dg {
64 : namespace Events {
65 : /// \cond
66 : template <size_t VolumeDim, typename Tensors,
67 : typename NonTensorComputeTagsList = tmpl::list<>,
68 : typename ArraySectionIdTag = void>
69 : class ObserveFields;
70 : /// \endcond
71 :
72 : /*!
73 : * \ingroup DiscontinuousGalerkinGroup
74 : * \brief %Observe volume tensor fields.
75 : *
76 : * A class that writes volume quantities to an h5 file during the simulation.
77 : * The observed quantitites are specified in the `VariablesToObserve` option.
78 : * Any `Tensor` in the `db::DataBox` can be observed but must be listed in the
79 : * `Tensors` template parameter. Any additional compute tags that hold a
80 : * `Tensor` can also be added to the `Tensors` template parameter. Finally,
81 : * `Variables` and other non-tensor compute tags can be listed in the
82 : * `NonTensorComputeTags` to facilitate observing. Note that the
83 : * `InertialCoordinates` are always observed.
84 : *
85 : * The user may specify an `interpolation_mesh` to which the
86 : * data is interpolated.
87 : *
88 : * \note The `NonTensorComputeTags` are intended to be used for `Variables`
89 : * compute tags like `Tags::DerivCompute`
90 : *
91 : * \par Array sections
92 : * This event supports sections (see `Parallel::Section`). Set the
93 : * `ArraySectionIdTag` template parameter to split up observations into subsets
94 : * of elements. The `observers::Tags::ObservationKey<ArraySectionIdTag>` must be
95 : * available in the DataBox. It identifies the section and is used as a suffix
96 : * for the path in the output file.
97 : */
98 : template <size_t VolumeDim, typename... Tensors,
99 : typename... NonTensorComputeTags, typename ArraySectionIdTag>
100 1 : class ObserveFields<VolumeDim, tmpl::list<Tensors...>,
101 : tmpl::list<NonTensorComputeTags...>, ArraySectionIdTag>
102 : : public Event {
103 : public:
104 : /// The name of the subfile inside the HDF5 file
105 1 : struct SubfileName {
106 0 : using type = std::string;
107 0 : static constexpr Options::String help = {
108 : "The name of the subfile inside the HDF5 file without an extension and "
109 : "without a preceding '/'."};
110 : };
111 :
112 : /// \cond
113 : explicit ObserveFields(CkMigrateMessage* /*unused*/) {}
114 : using PUP::able::register_constructor;
115 : WRAPPED_PUPable_decl_template(ObserveFields); // NOLINT
116 : /// \endcond
117 :
118 0 : struct VariablesToObserve {
119 0 : static constexpr Options::String help = "Subset of variables to observe";
120 0 : using type = std::vector<std::string>;
121 0 : static size_t lower_bound_on_size() { return 1; }
122 : };
123 :
124 0 : struct InterpolateToMesh {
125 0 : using type = Options::Auto<Mesh<VolumeDim>, Options::AutoLabel::None>;
126 0 : static constexpr Options::String help =
127 : "An optional mesh to which the variables are interpolated. This mesh "
128 : "specifies any number of collocation points, basis, and quadrature on "
129 : "which the observed quantities are evaluated. If no mesh is given, the "
130 : "results will be evaluated on the mesh the simulation runs on. The "
131 : "user may add several ObserveField Events e.g. with and without an "
132 : "interpolating mesh to output the data both on the original mesh and "
133 : "on a new mesh.";
134 : };
135 :
136 : /// The floating point type/precision with which to write the data to disk.
137 : ///
138 : /// Must be specified once for all data or individually for each variable
139 : /// being observed.
140 1 : struct FloatingPointTypes {
141 0 : static constexpr Options::String help =
142 : "The floating point type/precision with which to write the data to "
143 : "disk.\n\n"
144 : "Must be specified once for all data or individually for each "
145 : "variable being observed.";
146 0 : using type = std::vector<FloatingPointType>;
147 0 : static size_t upper_bound_on_size() { return sizeof...(Tensors); }
148 0 : static size_t lower_bound_on_size() { return 1; }
149 : };
150 :
151 : /// The floating point type/precision with which to write the coordinates to
152 : /// disk.
153 1 : struct CoordinatesFloatingPointType {
154 0 : static constexpr Options::String help =
155 : "The floating point type/precision with which to write the coordinates "
156 : "to disk.";
157 0 : using type = FloatingPointType;
158 : };
159 :
160 0 : using options =
161 : tmpl::list<SubfileName, CoordinatesFloatingPointType, FloatingPointTypes,
162 : VariablesToObserve, InterpolateToMesh>;
163 :
164 0 : static constexpr Options::String help =
165 : "Observe volume tensor fields.\n"
166 : "\n"
167 : "Writes volume quantities:\n"
168 : " * InertialCoordinates\n"
169 : " * Tensors listed in the 'VariablesToObserve' option\n";
170 :
171 0 : ObserveFields() = default;
172 :
173 0 : ObserveFields(const std::string& subfile_name,
174 : FloatingPointType coordinates_floating_point_type,
175 : const std::vector<FloatingPointType>& floating_point_types,
176 : const std::vector<std::string>& variables_to_observe,
177 : std::optional<Mesh<VolumeDim>> interpolation_mesh = {},
178 : const Options::Context& context = {});
179 :
180 0 : using compute_tags_for_observation_box =
181 : tmpl::list<Tensors..., NonTensorComputeTags...>;
182 :
183 0 : using return_tags = tmpl::list<>;
184 0 : using argument_tags = tmpl::list<::Tags::ObservationBox,
185 : ::Events::Tags::ObserverMesh<VolumeDim>>;
186 :
187 : template <typename DataBoxType, typename ComputeTagsList,
188 : typename Metavariables, typename ParallelComponent>
189 0 : void operator()(const ObservationBox<DataBoxType, ComputeTagsList>& box,
190 : const Mesh<VolumeDim>& mesh,
191 : Parallel::GlobalCache<Metavariables>& cache,
192 : const ElementId<VolumeDim>& array_index,
193 : const ParallelComponent* const component,
194 : const ObservationValue& observation_value) const {
195 : // Skip observation on elements that are not part of a section
196 : const std::optional<std::string> section_observation_key =
197 : observers::get_section_observation_key<ArraySectionIdTag>(box);
198 : if (not section_observation_key.has_value()) {
199 : return;
200 : }
201 : call_operator_impl(subfile_path_ + *section_observation_key,
202 : variables_to_observe_, interpolation_mesh_, mesh, box,
203 : cache, array_index, component, observation_value);
204 : }
205 :
206 : // We factor out the work into a static member function so it can be shared
207 : // with other field observing events, like the one that deals with DG-subcell
208 : // where there are two grids. This is to avoid copy-pasting all of the code.
209 : template <typename DataBoxType, typename ComputeTagsList,
210 : typename Metavariables, typename ParallelComponent>
211 0 : static void call_operator_impl(
212 : const std::string& subfile_path,
213 : const std::unordered_map<std::string, FloatingPointType>&
214 : variables_to_observe,
215 : const std::optional<Mesh<VolumeDim>>& interpolation_mesh,
216 : const Mesh<VolumeDim>& mesh,
217 : const ObservationBox<DataBoxType, ComputeTagsList>& box,
218 : Parallel::GlobalCache<Metavariables>& cache,
219 : const ElementId<VolumeDim>& element_id,
220 : const ParallelComponent* const /*meta*/,
221 : const ObservationValue& observation_value) {
222 : // if no interpolation_mesh is provided, the interpolation is essentially
223 : // ignored by the RegularGridInterpolant except for a single copy.
224 : const intrp::RegularGrid interpolant(mesh,
225 : interpolation_mesh.value_or(mesh));
226 :
227 : // Remove tensor types, only storing individual components.
228 : std::vector<TensorComponent> components;
229 : // This is larger than we need if we are only observing some
230 : // tensors, but that's not a big deal and calculating the correct
231 : // size is nontrivial.
232 : components.reserve(alg::accumulate(
233 : std::initializer_list<size_t>{
234 : std::decay_t<decltype(value(typename Tensors::type{}))>::size()...},
235 : 0_st));
236 :
237 : const auto record_tensor_component_impl =
238 : [&components, &interpolant](const auto& tensor,
239 : const FloatingPointType floating_point_type,
240 : const std::string& tag_name) {
241 : for (size_t i = 0; i < tensor.size(); ++i) {
242 : const auto tensor_component = interpolant.interpolate(tensor[i]);
243 : if (floating_point_type == FloatingPointType::Float) {
244 : components.emplace_back(
245 : tag_name + tensor.component_suffix(i),
246 : std::vector<float>{tensor_component.begin(),
247 : tensor_component.end()});
248 : } else {
249 : components.emplace_back(tag_name + tensor.component_suffix(i),
250 : tensor_component);
251 : }
252 : }
253 : };
254 : const auto record_tensor_components =
255 : [&box, &record_tensor_component_impl,
256 : &variables_to_observe](const auto tensor_tag_v) {
257 : using tensor_tag = tmpl::type_from<decltype(tensor_tag_v)>;
258 : const std::string tag_name = db::tag_name<tensor_tag>();
259 : if (const auto var_to_observe = variables_to_observe.find(tag_name);
260 : var_to_observe != variables_to_observe.end()) {
261 : const auto& tensor = get<tensor_tag>(box);
262 : if (not has_value(tensor)) {
263 : // This will only print a warning the first time it's called on a
264 : // node.
265 : [[maybe_unused]] static bool t =
266 : ObserveFields::print_warning_about_optional<tensor_tag>();
267 : return;
268 : }
269 : const auto floating_point_type = var_to_observe->second;
270 : record_tensor_component_impl(value(tensor), floating_point_type,
271 : tag_name);
272 : }
273 : };
274 : EXPAND_PACK_LEFT_TO_RIGHT(record_tensor_components(tmpl::type_<Tensors>{}));
275 :
276 : const Parallel::ArrayComponentId array_component_id{
277 : std::add_pointer_t<ParallelComponent>{nullptr},
278 : Parallel::ArrayIndex<ElementId<VolumeDim>>{element_id}};
279 : ElementVolumeData element_volume_data{element_id, std::move(components),
280 : interpolation_mesh.value_or(mesh)};
281 : observers::ObservationId observation_id{observation_value.value,
282 : subfile_path + ".vol"};
283 :
284 : auto& local_observer = *Parallel::local_branch(
285 : Parallel::get_parallel_component<
286 : tmpl::conditional_t<Parallel::is_nodegroup_v<ParallelComponent>,
287 : observers::ObserverWriter<Metavariables>,
288 : observers::Observer<Metavariables>>>(cache));
289 :
290 : if constexpr (Parallel::is_nodegroup_v<ParallelComponent>) {
291 : // Send data to reduction observer writer (nodegroup)
292 : std::unordered_map<Parallel::ArrayComponentId,
293 : std::vector<ElementVolumeData>>
294 : data_to_send{};
295 : data_to_send[array_component_id] =
296 : std::vector{std::move(element_volume_data)};
297 : Parallel::threaded_action<
298 : observers::ThreadedActions::ContributeVolumeDataToWriter>(
299 : local_observer, std::move(observation_id), array_component_id,
300 : subfile_path, std::move(data_to_send));
301 : } else {
302 : // Send data to volume observer
303 : Parallel::simple_action<observers::Actions::ContributeVolumeData>(
304 : local_observer, std::move(observation_id), subfile_path,
305 : array_component_id, std::move(element_volume_data));
306 : }
307 : }
308 :
309 0 : using observation_registration_tags = tmpl::list<::Tags::DataBox>;
310 :
311 : template <typename DbTagsList>
312 : std::optional<
313 : std::pair<observers::TypeOfObservation, observers::ObservationKey>>
314 0 : get_observation_type_and_key_for_registration(
315 : const db::DataBox<DbTagsList>& box) const {
316 : const std::optional<std::string> section_observation_key =
317 : observers::get_section_observation_key<ArraySectionIdTag>(box);
318 : if (not section_observation_key.has_value()) {
319 : return std::nullopt;
320 : }
321 : return {{observers::TypeOfObservation::Volume,
322 : observers::ObservationKey{
323 : subfile_path_ + section_observation_key.value() + ".vol"}}};
324 : }
325 :
326 0 : using is_ready_argument_tags = tmpl::list<>;
327 :
328 : template <typename Metavariables, typename ArrayIndex, typename Component>
329 0 : bool is_ready(Parallel::GlobalCache<Metavariables>& /*cache*/,
330 : const ArrayIndex& /*array_index*/,
331 : const Component* const /*meta*/) const {
332 : return true;
333 : }
334 :
335 1 : bool needs_evolved_variables() const override { return true; }
336 :
337 : // NOLINTNEXTLINE(google-runtime-references)
338 0 : void pup(PUP::er& p) override {
339 : Event::pup(p);
340 : p | subfile_path_;
341 : p | variables_to_observe_;
342 : p | interpolation_mesh_;
343 : }
344 :
345 : private:
346 : template <typename Tag>
347 0 : static bool print_warning_about_optional() {
348 : Parallel::printf(
349 : "Warning: ObserveFields is trying to dump the tag %s "
350 : "but it is stored as a std::optional and has not been "
351 : "evaluated. This most commonly occurs when you are "
352 : "trying to either observe an analytic solution or errors when "
353 : "no analytic solution is available.\n",
354 : db::tag_name<Tag>());
355 : return false;
356 : }
357 :
358 0 : std::string subfile_path_;
359 0 : std::unordered_map<std::string, FloatingPointType> variables_to_observe_{};
360 0 : std::optional<Mesh<VolumeDim>> interpolation_mesh_{};
361 : };
362 :
363 : template <size_t VolumeDim, typename... Tensors,
364 : typename... NonTensorComputeTags, typename ArraySectionIdTag>
365 : ObserveFields<VolumeDim, tmpl::list<Tensors...>,
366 : tmpl::list<NonTensorComputeTags...>, ArraySectionIdTag>::
367 : ObserveFields(const std::string& subfile_name,
368 : const FloatingPointType coordinates_floating_point_type,
369 : const std::vector<FloatingPointType>& floating_point_types,
370 : const std::vector<std::string>& variables_to_observe,
371 : std::optional<Mesh<VolumeDim>> interpolation_mesh,
372 : const Options::Context& context)
373 : : subfile_path_("/" + subfile_name),
374 : variables_to_observe_([&context, &floating_point_types,
375 : &variables_to_observe]() {
376 : if (floating_point_types.size() != 1 and
377 : floating_point_types.size() != variables_to_observe.size()) {
378 : PARSE_ERROR(context, "The number of floating point types specified ("
379 : << floating_point_types.size()
380 : << ") must be 1 or the number of variables "
381 : "specified for observing ("
382 : << variables_to_observe.size() << ")");
383 : }
384 : std::unordered_map<std::string, FloatingPointType> result{};
385 : for (size_t i = 0; i < variables_to_observe.size(); ++i) {
386 : result[variables_to_observe[i]] = floating_point_types.size() == 1
387 : ? floating_point_types[0]
388 : : floating_point_types[i];
389 : ASSERT(
390 : result.at(variables_to_observe[i]) == FloatingPointType::Float or
391 : result.at(variables_to_observe[i]) ==
392 : FloatingPointType::Double,
393 : "Floating point type for variable '"
394 : << variables_to_observe[i]
395 : << "' must be either Float or Double.");
396 : }
397 : return result;
398 : }()),
399 : interpolation_mesh_(interpolation_mesh) {
400 : ASSERT(
401 : (... or (db::tag_name<Tensors>() == "InertialCoordinates")),
402 : "There is no tag with name 'InertialCoordinates' specified "
403 : "for the observer. Please make sure you specify a tag in the 'Tensors' "
404 : "list that has the 'db::tag_name()' 'InertialCoordinates'.");
405 : db::validate_selection<tmpl::list<Tensors...>>(variables_to_observe, context);
406 : variables_to_observe_["InertialCoordinates"] =
407 : coordinates_floating_point_type;
408 : }
409 :
410 : /// \cond
411 : template <size_t VolumeDim, typename... Tensors,
412 : typename... NonTensorComputeTags, typename ArraySectionIdTag>
413 : PUP::able::PUP_ID ObserveFields<VolumeDim, tmpl::list<Tensors...>,
414 : tmpl::list<NonTensorComputeTags...>,
415 : ArraySectionIdTag>::my_PUP_ID = 0; // NOLINT
416 : /// \endcond
417 : } // namespace Events
418 : } // namespace dg
|
