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