Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <algorithm>
7 : #include <cstddef>
8 : #include <optional>
9 : #include <string>
10 : #include <tuple>
11 : #include <variant>
12 : #include <vector>
13 :
14 : #include "DataStructures/DataBox/DataBox.hpp"
15 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
16 : #include "DataStructures/DataBox/Tag.hpp"
17 : #include "DataStructures/DataVector.hpp"
18 : #include "DataStructures/TaggedTuple.hpp"
19 : #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
20 : #include "DataStructures/Tensor/Tensor.hpp"
21 : #include "Domain/Domain.hpp"
22 : #include "Domain/ElementMap.hpp"
23 : #include "Domain/Structure/ElementId.hpp"
24 : #include "IO/Exporter/Exporter.hpp"
25 : #include "IO/Exporter/PointwiseInterpolator.hpp"
26 : #include "IO/Exporter/SelectObservation.hpp"
27 : #include "IO/H5/AccessType.hpp"
28 : #include "IO/H5/File.hpp"
29 : #include "IO/H5/TensorData.hpp"
30 : #include "IO/H5/VolumeData.hpp"
31 : #include "IO/Importers/ObservationSelector.hpp"
32 : #include "IO/Importers/Tags.hpp"
33 : #include "NumericalAlgorithms/Interpolation/RegularGridInterpolant.hpp"
34 : #include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp"
35 : #include "Parallel/AlgorithmExecution.hpp"
36 : #include "Parallel/ArrayCollection/IsDgElementCollection.hpp"
37 : #include "Parallel/ArrayComponentId.hpp"
38 : #include "Parallel/ArrayIndex.hpp"
39 : #include "Parallel/GlobalCache.hpp"
40 : #include "Parallel/Invoke.hpp"
41 : #include "Utilities/EqualWithinRoundoff.hpp"
42 : #include "Utilities/ErrorHandling/Assert.hpp"
43 : #include "Utilities/ErrorHandling/CaptureForError.hpp"
44 : #include "Utilities/ErrorHandling/Error.hpp"
45 : #include "Utilities/FileSystem.hpp"
46 : #include "Utilities/Gsl.hpp"
47 : #include "Utilities/Literals.hpp"
48 : #include "Utilities/Overloader.hpp"
49 : #include "Utilities/Requires.hpp"
50 : #include "Utilities/Serialization/Serialize.hpp"
51 : #include "Utilities/TMPL.hpp"
52 :
53 1 : namespace importers {
54 :
55 : /// \cond
56 : template <typename Metavariables>
57 : struct ElementDataReader;
58 : namespace Actions {
59 : template <size_t Dim, typename FieldTagsList, typename ReceiveComponent>
60 : struct ReadAllVolumeDataAndDistribute;
61 : } // namespace Actions
62 : /// \endcond
63 :
64 1 : namespace Tags {
65 : /*!
66 : * \brief Indicates an available tensor field is selected for importing, along
67 : * with the name of the dataset in the volume data file.
68 : *
69 : * Set the value to a dataset name to import the `FieldTag` from that dataset,
70 : * or to `std::nullopt` to skip importing the `FieldTag`. The dataset name
71 : * excludes tensor component suffixes like "_x" or "_xy". These suffixes will be
72 : * added automatically. A sensible value for the dataset name is often
73 : * `db::tag_name<FieldTag>()`, but the user should generally be given the
74 : * opportunity to set the dataset name in the input file.
75 : */
76 : template <typename FieldTag>
77 1 : struct Selected : db::SimpleTag {
78 0 : using type = std::optional<std::string>;
79 : };
80 : } // namespace Tags
81 :
82 : namespace detail {
83 :
84 : // Translate the importer's observation selection into the Exporter's
85 : // `ObservationVariant`
86 : inline spectre::Exporter::ObservationVariant observation_variant(
87 : const std::variant<double, ObservationSelector>& observation_value,
88 : const std::optional<double>& observation_value_epsilon) {
89 : return std::visit(
90 : Overloader{[&observation_value_epsilon](const double local_obs_value)
91 : -> spectre::Exporter::ObservationVariant {
92 : if (observation_value_epsilon.has_value()) {
93 : return spectre::Exporter::ObservationValue{
94 : local_obs_value, observation_value_epsilon.value()};
95 : }
96 : return local_obs_value;
97 : },
98 : [](const ObservationSelector local_obs_selector)
99 : -> spectre::Exporter::ObservationVariant {
100 : switch (local_obs_selector) {
101 : case ObservationSelector::First:
102 : return spectre::Exporter::ObservationStep{0};
103 : case ObservationSelector::Last:
104 : return spectre::Exporter::ObservationStep{-1};
105 : default:
106 : ERROR("Unknown importers::ObservationSelector: "
107 : << local_obs_selector);
108 : }
109 : }},
110 : observation_value);
111 : }
112 :
113 : // Read the single `tensor_name` from the `volume_file`, taking care of suffixes
114 : // like "_x" etc for its components.
115 : template <typename TensorType>
116 : void read_tensor_data(const gsl::not_null<TensorType*> tensor_data,
117 : const std::string& tensor_name,
118 : const h5::VolumeData& volume_file,
119 : const size_t observation_id) {
120 : for (size_t i = 0; i < tensor_data->size(); ++i) {
121 : const auto& tensor_component = volume_file.get_tensor_component(
122 : observation_id, tensor_name + tensor_data->component_suffix(
123 : tensor_data->get_tensor_index(i)));
124 : if (not std::holds_alternative<DataVector>(tensor_component.data)) {
125 : ERROR("The tensor component '"
126 : << tensor_component.name
127 : << "' is not a double-precision DataVector. Reading in "
128 : "single-precision volume data is not supported.");
129 : }
130 : (*tensor_data)[i] = std::get<DataVector>(tensor_component.data);
131 : }
132 : }
133 :
134 : // Read the `selected_fields` from the `volume_file`. Reads the data
135 : // for all elements in the `volume_file` at once. Invoked lazily when data
136 : // for an element in the volume file is needed.
137 : template <typename FieldTagsList>
138 : tuples::tagged_tuple_from_typelist<FieldTagsList> read_tensor_data(
139 : const h5::VolumeData& volume_file, const size_t observation_id,
140 : const tuples::tagged_tuple_from_typelist<
141 : db::wrap_tags_in<Tags::Selected, FieldTagsList>>& selected_fields) {
142 : tuples::tagged_tuple_from_typelist<FieldTagsList> all_tensor_data{};
143 : tmpl::for_each<FieldTagsList>([&all_tensor_data, &volume_file,
144 : &observation_id,
145 : &selected_fields](auto field_tag_v) {
146 : using field_tag = tmpl::type_from<decltype(field_tag_v)>;
147 : const auto& selection = get<Tags::Selected<field_tag>>(selected_fields);
148 : if (not selection.has_value()) {
149 : return;
150 : }
151 : read_tensor_data(make_not_null(&get<field_tag>(all_tensor_data)),
152 : selection.value(), volume_file, observation_id);
153 : });
154 : return all_tensor_data;
155 : }
156 :
157 : // Extract this element's data from the read-in dataset
158 : template <typename FieldTagsList>
159 : tuples::tagged_tuple_from_typelist<FieldTagsList> extract_element_data(
160 : const std::pair<size_t, size_t>& element_data_offset_and_length,
161 : const tuples::tagged_tuple_from_typelist<FieldTagsList>& all_tensor_data,
162 : const tuples::tagged_tuple_from_typelist<
163 : db::wrap_tags_in<Tags::Selected, FieldTagsList>>& selected_fields) {
164 : tuples::tagged_tuple_from_typelist<FieldTagsList> element_data{};
165 : tmpl::for_each<FieldTagsList>(
166 : [&element_data, &offset = element_data_offset_and_length.first,
167 : &num_points = element_data_offset_and_length.second, &all_tensor_data,
168 : &selected_fields](auto field_tag_v) {
169 : using field_tag = tmpl::type_from<decltype(field_tag_v)>;
170 : const auto& selection = get<Tags::Selected<field_tag>>(selected_fields);
171 : if (not selection.has_value()) {
172 : return;
173 : }
174 : auto& element_tensor_data = get<field_tag>(element_data);
175 : // Iterate independent components of the tensor
176 : for (size_t i = 0; i < element_tensor_data.size(); ++i) {
177 : const DataVector& data_tensor_component =
178 : get<field_tag>(all_tensor_data)[i];
179 : DataVector element_tensor_component{num_points};
180 : // Retrieve data from slice of the contigious dataset
181 : for (size_t j = 0; j < element_tensor_component.size(); ++j) {
182 : element_tensor_component[j] = data_tensor_component[offset + j];
183 : }
184 : element_tensor_data[i] = element_tensor_component;
185 : }
186 : });
187 : return element_data;
188 : }
189 :
190 : // Check that the inertial coordinates computed with the given domain are the
191 : // same as the ones passed to this function.
192 : // This is important to avoid hard-to-find bugs where data is loaded
193 : // to the wrong coordinates. For example, if the evolution domain deforms the
194 : // excision surfaces a bit but the initial data doesn't, then it would be wrong
195 : // to load the initial data to the evolution grid without an interpolation.
196 : template <size_t Dim>
197 : void verify_inertial_coordinates(
198 : const Domain<Dim>& domain, const double time,
199 : const domain::FunctionsOfTimeMap& functions_of_time,
200 : const ElementId<Dim>& element_id, const Mesh<Dim>& mesh,
201 : const tnsr::I<DataVector, Dim, Frame::Inertial>& inertial_coords) {
202 : const auto logical_coords = logical_coordinates(mesh);
203 : ElementMap<Dim, Frame::Inertial> element_map{
204 : element_id, domain.blocks()[element_id.block_id()]};
205 : const auto mapped_inertial_coords =
206 : element_map(logical_coords, time, functions_of_time);
207 : const double scale = blaze::max(get(magnitude(mapped_inertial_coords)));
208 : if (not equal_within_roundoff(mapped_inertial_coords, inertial_coords,
209 : std::numeric_limits<double>::epsilon() * 100.0,
210 : scale)) {
211 : DataVector diff =
212 : square(get<0>(inertial_coords) - get<0>(mapped_inertial_coords));
213 : for (size_t d = 1; d < Dim; ++d) {
214 : diff += square(inertial_coords.get(d) - mapped_inertial_coords.get(d));
215 : }
216 : diff = sqrt(diff);
217 : const double max_coord_distance = blaze::max(diff);
218 : CAPTURE_FOR_ERROR(element_id);
219 : CAPTURE_FOR_ERROR(max_coord_distance);
220 : CAPTURE_FOR_ERROR(scale);
221 : ERROR_NO_TRACE(
222 : "The source and target domain don't match. Set 'ElementsAreIdentical: "
223 : "False' to enable interpolation between the grids.");
224 : }
225 : }
226 :
227 : // Interpolate only the `selected_fields` in `source_element_data` to the
228 : // `target_mesh` (used when elements differ only by p-refinement)
229 : template <typename FieldTagsList, size_t Dim>
230 : void interpolate_selected_fields(
231 : const gsl::not_null<tuples::tagged_tuple_from_typelist<FieldTagsList>*>
232 : target_element_data,
233 : const tuples::tagged_tuple_from_typelist<FieldTagsList>&
234 : source_element_data,
235 : const Mesh<Dim>& source_mesh, const Mesh<Dim>& target_mesh,
236 : const tuples::tagged_tuple_from_typelist<
237 : db::wrap_tags_in<Tags::Selected, FieldTagsList>>& selected_fields) {
238 : const intrp::RegularGrid<Dim> interpolator{source_mesh, target_mesh};
239 : tmpl::for_each<FieldTagsList>([&source_element_data, &target_element_data,
240 : &interpolator,
241 : &selected_fields](auto field_tag_v) {
242 : using field_tag = tmpl::type_from<decltype(field_tag_v)>;
243 : const auto& selection = get<Tags::Selected<field_tag>>(selected_fields);
244 : if (not selection.has_value()) {
245 : return;
246 : }
247 : const auto& source_tensor_data = get<field_tag>(source_element_data);
248 : auto& target_tensor_data = get<field_tag>(*target_element_data);
249 : // Iterate independent components of the tensor
250 : for (size_t i = 0; i < source_tensor_data.size(); ++i) {
251 : const DataVector& source_tensor_component = source_tensor_data[i];
252 : DataVector& target_tensor_component = target_tensor_data[i];
253 : // Interpolate
254 : interpolator.interpolate(make_not_null(&target_tensor_component),
255 : source_tensor_component);
256 : }
257 : });
258 : }
259 :
260 : // Scatter the slice `[start, start + num_points)` of the flat,
261 : // component-indexed `interpolated_data` (as produced by
262 : // `spectre::Exporter::interpolate_to_points`) into a tagged tuple of tensors
263 : // for a single target element. The components are laid out in `FieldTagsList`
264 : // order, skipping unselected fields, matching the order of the
265 : // `tensor_components` passed to the interpolation.
266 : template <typename FieldTagsList>
267 : tuples::tagged_tuple_from_typelist<FieldTagsList> scatter_element_data(
268 : const std::vector<DataVector>& interpolated_data, const size_t start,
269 : const size_t num_points,
270 : const tuples::tagged_tuple_from_typelist<
271 : db::wrap_tags_in<Tags::Selected, FieldTagsList>>& selected_fields) {
272 : tuples::tagged_tuple_from_typelist<FieldTagsList> element_data{};
273 : size_t component_index = 0;
274 : tmpl::for_each<FieldTagsList>([&element_data, &interpolated_data, &start,
275 : &num_points, &selected_fields,
276 : &component_index](auto field_tag_v) {
277 : using field_tag = tmpl::type_from<decltype(field_tag_v)>;
278 : if (not get<Tags::Selected<field_tag>>(selected_fields).has_value()) {
279 : return;
280 : }
281 : auto& element_tensor_data = get<field_tag>(element_data);
282 : for (size_t i = 0; i < element_tensor_data.size(); ++i) {
283 : DataVector component{num_points};
284 : const DataVector& interpolated_component =
285 : interpolated_data[component_index];
286 : for (size_t j = 0; j < num_points; ++j) {
287 : component[j] = interpolated_component[start + j];
288 : }
289 : element_tensor_data[i] = std::move(component);
290 : ++component_index;
291 : }
292 : });
293 : return element_data;
294 : }
295 :
296 : } // namespace detail
297 :
298 0 : namespace Actions {
299 :
300 : /*!
301 : * \brief Read a volume data file and distribute the data to all registered
302 : * elements, interpolating to the target points if needed.
303 : *
304 : * \note Use this action if you want to quickly load and distribute volume data.
305 : * If you need to beyond that (such as more control over input-file options),
306 : * write a new action and dispatch to
307 : * `importers::Actions::ReadAllVolumeDataAndDistribute`.
308 : *
309 : * \details Invoke this action on the elements of an array parallel component to
310 : * dispatch reading the volume data file specified by options placed in the
311 : * `ImporterOptionsGroup`. The tensors in `FieldTagsList` will be loaded from
312 : * the file and distributed to all elements that have previously registered. Use
313 : * `importers::Actions::RegisterWithElementDataReader` to register the elements
314 : * of the array parallel component in a previous phase.
315 : *
316 : * Note that the volume data file will only be read once per node, triggered by
317 : * the first element that invokes this action. All subsequent invocations of
318 : * this action on the node will do nothing. See
319 : * `importers::Actions::ReadAllVolumeDataAndDistribute` for details.
320 : *
321 : * The data is distributed to the elements using `Parallel::receive_data`. The
322 : * elements can monitor `importers::Tags::VolumeData` in their inbox to wait for
323 : * the data and process it once it's available. We provide the action
324 : * `importers::Actions::ReceiveVolumeData` that waits for the data and moves it
325 : * directly into the DataBox. You can also implement a specialized action that
326 : * might verify and post-process the data before populating the DataBox.
327 : *
328 : * \see Dev guide on \ref dev_guide_importing
329 : */
330 : template <typename ImporterOptionsGroup, typename FieldTagsList>
331 1 : struct ReadVolumeData {
332 0 : using const_global_cache_tags =
333 : tmpl::list<Tags::ImporterOptions<ImporterOptionsGroup>>;
334 :
335 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
336 : size_t Dim, typename ActionList, typename ParallelComponent>
337 0 : static Parallel::iterable_action_return_t apply(
338 : db::DataBox<DbTagsList>& /*box*/,
339 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
340 : Parallel::GlobalCache<Metavariables>& cache,
341 : const ElementId<Dim>& /*array_index*/, const ActionList /*meta*/,
342 : const ParallelComponent* const /*meta*/) {
343 : // Not using `ckLocalBranch` here to make sure the simple action invocation
344 : // is asynchronous.
345 : auto& reader_component = Parallel::get_parallel_component<
346 : importers::ElementDataReader<Metavariables>>(cache);
347 : Parallel::simple_action<importers::Actions::ReadAllVolumeDataAndDistribute<
348 : Dim, FieldTagsList, ParallelComponent>>(
349 : reader_component,
350 : get<Tags::ImporterOptions<ImporterOptionsGroup>>(cache), 0_st);
351 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
352 : }
353 : };
354 :
355 : /*!
356 : * \brief Read a volume data file and distribute the data to all registered
357 : * elements, interpolating to the target points if needed.
358 : *
359 : * This action can be invoked on the `importers::ElementDataReader` component
360 : * once all elements have been registered with it. It opens the data file, reads
361 : * the data for each registered element and uses `Parallel::receive_data` to
362 : * distribute the data to the elements. The elements can monitor
363 : * `importers::Tags::VolumeData` in their inbox to wait for the data and process
364 : * it once it's available. You can use `importers::Actions::ReceiveVolumeData`
365 : * to wait for the data and move it directly into the DataBox, or implement a
366 : * specialized action that might verify and post-process the data.
367 : *
368 : * Note that instead of invoking this action directly on the
369 : * `importers::ElementDataReader` component you can invoke the iterable action
370 : * `importers::Actions::ReadVolumeData` on the elements of an array parallel
371 : * component for simple use cases.
372 : *
373 : * - Pass along the following arguments to the simple action invocation:
374 : * - `options`: `importers::ImporterOptions` that specify the H5 files
375 : * with volume data to load.
376 : * - `volume_data_id`: A number (or hash) that identifies this import
377 : * operation. Will also be used to identify the loaded volume data in the
378 : * inbox of the receiving elements.
379 : * - `selected_fields` (optional): See below.
380 : * - The `FieldTagsList` parameter specifies a typelist of tensor tags that
381 : * can be read from the file and provided to each element. The subset of tensors
382 : * that will actually be read and distributed can be selected at runtime with
383 : * the `selected_fields` argument that is passed to this simple action. See
384 : * importers::Tags::Selected for details. By default, all tensors in the
385 : * `FieldTagsList` are selected, and read from datasets named
386 : * `db::tag_name<Tag>() + suffix`, where the `suffix` is empty for scalars, or
387 : * `"_"` followed by the `Tensor::component_name` for each independent tensor
388 : * component.
389 : * - `Parallel::receive_data` is invoked on each registered element of the
390 : * `ReceiveComponent` to populate `importers::Tags::VolumeData` in the element's
391 : * inbox with a `tuples::tagged_tuple_from_typelist<FieldTagsList>` containing
392 : * the tensor data for that element. The `ReceiveComponent` must the the same
393 : * that was encoded into the `Parallel::ArrayComponentId` used to register the
394 : * elements. The `volume_data_id` passed to this action is used as key.
395 : *
396 : * \par Memory consumption
397 : * This action runs once on every node. Volume data files are loaded one at a
398 : * time, so memory consumption does _not_ grow with the the number of source
399 : * files. All coordinates of elements on this node and their interpolated data
400 : * is held in memory at once, so memory consumption scales with the number of
401 : * elements on this node.
402 : *
403 : * \see Dev guide on \ref dev_guide_importing
404 : */
405 : template <size_t Dim, typename FieldTagsList, typename ReceiveComponent>
406 1 : struct ReadAllVolumeDataAndDistribute {
407 : template <typename ParallelComponent, typename DataBox,
408 : typename Metavariables, typename ArrayIndex>
409 0 : static void apply(DataBox& box, Parallel::GlobalCache<Metavariables>& cache,
410 : const ArrayIndex& /*array_index*/,
411 : const ImporterOptions& options, const size_t volume_data_id,
412 : tuples::tagged_tuple_from_typelist<
413 : db::wrap_tags_in<Tags::Selected, FieldTagsList>>
414 : selected_fields = select_all_fields(FieldTagsList{})) {
415 : const bool elements_are_identical =
416 : get<OptionTags::ElementsAreIdentical>(options);
417 :
418 : // Only read and distribute the volume data once
419 : // This action will be invoked by `importers::Actions::ReadVolumeData` from
420 : // every element on the node, but only the first invocation reads the file
421 : // and distributes the data to all elements. Subsequent invocations do
422 : // nothing. The `volume_data_id` identifies whether or not we have already
423 : // read the requested data. Doing this at runtime avoids having to collect
424 : // all data files that will be read in at compile-time to initialize a flag
425 : // in the DataBox for each of them.
426 : const auto& has_read_volume_data =
427 : db::get<Tags::ElementDataAlreadyRead>(box);
428 : if (has_read_volume_data.find(volume_data_id) !=
429 : has_read_volume_data.end()) {
430 : return;
431 : }
432 : db::mutate<Tags::ElementDataAlreadyRead>(
433 : [&volume_data_id](const auto local_has_read_volume_data) {
434 : local_has_read_volume_data->insert(volume_data_id);
435 : },
436 : make_not_null(&box));
437 :
438 : // This is the subset of elements that reside on this node. They have
439 : // registered themselves before. Our job is to fill them with volume data.
440 : std::unordered_set<ElementId<Dim>> target_element_ids{};
441 : for (const auto& target_element : get<Tags::RegisteredElements<Dim>>(box)) {
442 : const auto& element_array_component_id = target_element.first;
443 : const CkArrayIndex& raw_element_index =
444 : element_array_component_id.array_index();
445 : // Check if the parallel component of the registered element matches the
446 : // callback, because it's possible that elements from other components
447 : // with the same index are also registered.
448 : // Since the way the component is encoded in `ArrayComponentId` is
449 : // private to that class, we construct one and compare.
450 : // Can't use Parallel::make_array_component_id here because we need the
451 : // original array_index type, not a CkArrayIndex.
452 : if (element_array_component_id !=
453 : Parallel::ArrayComponentId(
454 : std::add_pointer_t<ReceiveComponent>{nullptr},
455 : raw_element_index)) {
456 : continue;
457 : }
458 : const auto target_element_id =
459 : Parallel::ArrayIndex<ElementId<Dim>>(raw_element_index).get_index();
460 : target_element_ids.insert(target_element_id);
461 : }
462 : if (UNLIKELY(target_element_ids.empty())) {
463 : return;
464 : }
465 :
466 : // Resolve the file glob
467 : const std::string& file_glob = get<OptionTags::FileGlob>(options);
468 : const std::vector<std::string> file_paths = file_system::glob(file_glob);
469 : if (file_paths.empty()) {
470 : ERROR_NO_TRACE("The file glob '" << file_glob << "' matches no files.");
471 : }
472 :
473 : // Select observation to read from each file
474 : const spectre::Exporter::ObservationVariant observation =
475 : detail::observation_variant(
476 : get<OptionTags::ObservationValue>(options),
477 : get<OptionTags::ObservationValueEpsilon>(options));
478 :
479 : // When interpolation between the source and target grids is needed, reuse
480 : // spectre::Exporter::interpolate_to_points to interpolate to all target
481 : // points on this node at once, then scatter the results to the elements.
482 : if (not elements_are_identical) {
483 : // Gather all target points on this node into a single contiguous tensor,
484 : // recording the range [start, start + num_points) of each target element.
485 : std::vector<ElementId<Dim>> target_ids{};
486 : std::vector<size_t> target_starts{};
487 : std::vector<size_t> target_num_points{};
488 : std::vector<const tnsr::I<DataVector, Dim, Frame::Inertial>*>
489 : target_coords{};
490 : target_ids.reserve(target_element_ids.size());
491 : target_starts.reserve(target_element_ids.size());
492 : target_num_points.reserve(target_element_ids.size());
493 : target_coords.reserve(target_element_ids.size());
494 : size_t total_num_points = 0;
495 : for (const auto& target_element_id : target_element_ids) {
496 : const auto& target_points =
497 : get<Tags::RegisteredElements<Dim>>(box)
498 : .at(Parallel::make_array_component_id<ReceiveComponent>(
499 : target_element_id))
500 : .first;
501 : target_ids.push_back(target_element_id);
502 : target_coords.push_back(&target_points);
503 : target_starts.push_back(total_num_points);
504 : target_num_points.push_back(target_points.begin()->size());
505 : total_num_points += target_num_points.back();
506 : }
507 : tnsr::I<DataVector, Dim, Frame::Inertial> all_target_points{
508 : total_num_points};
509 : for (size_t e = 0; e < target_ids.size(); ++e) {
510 : for (size_t d = 0; d < Dim; ++d) {
511 : for (size_t i = 0; i < target_num_points[e]; ++i) {
512 : all_target_points.get(d)[target_starts[e] + i] =
513 : target_coords[e]->get(d)[i];
514 : }
515 : }
516 : }
517 :
518 : // Flat list of dataset component names for the selected fields, in
519 : // `FieldTagsList` order. The layout matches
520 : // `detail::scatter_element_data`.
521 : std::vector<std::string> tensor_components{};
522 : tmpl::for_each<FieldTagsList>([&tensor_components,
523 : &selected_fields](auto field_tag_v) {
524 : using field_tag = tmpl::type_from<decltype(field_tag_v)>;
525 : const auto& selection = get<Tags::Selected<field_tag>>(selected_fields);
526 : if (not selection.has_value()) {
527 : return;
528 : }
529 : using TensorType = typename field_tag::type;
530 : for (size_t i = 0; i < TensorType::size(); ++i) {
531 : tensor_components.push_back(selection.value() +
532 : TensorType::component_suffix(i));
533 : }
534 : });
535 :
536 : // Interpolate all target points at once. Error if any target point lies
537 : // outside the source domain. This implementation opens each file in turn
538 : // (so it doesn't hold all files in memory at once), and it is efficient
539 : // about mapping points through the blocks of the source domain.
540 : std::vector<DataVector> interpolated_data{};
541 : spectre::Exporter::interpolate_to_points(
542 : make_not_null(&interpolated_data), file_paths,
543 : "/" + get<OptionTags::Subgroup>(options), observation,
544 : tensor_components, all_target_points,
545 : get<OptionTags::ExtrapolateIntoExcisions>(options),
546 : /*error_on_missing_points=*/true,
547 : get<OptionTags::NumThreads>(options));
548 : // The target points are no longer needed; free them before distributing
549 : // the (potentially large) interpolated data to the target elements.
550 : all_target_points = tnsr::I<DataVector, Dim, Frame::Inertial>{};
551 :
552 : // Distribute the interpolated data to the target elements.
553 : for (size_t e = 0; e < target_ids.size(); ++e) {
554 : auto target_element_data = detail::scatter_element_data<FieldTagsList>(
555 : interpolated_data, target_starts[e], target_num_points[e],
556 : selected_fields);
557 : if constexpr (Parallel::is_dg_element_collection_v<ReceiveComponent>) {
558 : ERROR("Can't yet do numerical initial data with nodegroups");
559 : } else {
560 : Parallel::receive_data<Tags::VolumeData<FieldTagsList>>(
561 : Parallel::get_parallel_component<ReceiveComponent>(
562 : cache)[target_ids[e]],
563 : volume_data_id, std::move(target_element_data));
564 : }
565 : }
566 : return;
567 : } // not elements_are_identical
568 :
569 : // Now handle identical elements:
570 : // The source and target elements are the same (matching domains and
571 : // h-refinement), so data is transferred one-to-one, interpolating only
572 : // between different meshes (p-refinement).
573 : std::optional<size_t> prev_observation_id{};
574 : double observation_value = std::numeric_limits<double>::signaling_NaN();
575 : std::optional<Domain<Dim>> source_domain{};
576 : domain::FunctionsOfTimeMap source_domain_functions_of_time{};
577 : for (const std::string& file_name : file_paths) {
578 : // Open the volume data file
579 : h5::H5File<h5::AccessType::ReadOnly> h5file(file_name);
580 : constexpr size_t version_number = 0;
581 : const auto& volume_file = h5file.get<h5::VolumeData>(
582 : "/" + get<OptionTags::Subgroup>(options), version_number);
583 :
584 : // Select observation ID
585 : const size_t observation_id = std::visit(
586 : spectre::Exporter::SelectObservation{volume_file}, observation);
587 : if (prev_observation_id.has_value() and
588 : prev_observation_id.value() != observation_id) {
589 : ERROR("Inconsistent selection of observation ID in file "
590 : << file_name
591 : << ". Make sure all files select the same observation ID.");
592 : }
593 : prev_observation_id = observation_id;
594 : observation_value = volume_file.get_observation_value(observation_id);
595 :
596 : // Memory buffer for the tensor data stored in this file. The data is
597 : // loaded lazily when it is needed, so we can skip loading files that
598 : // contain none of the elements on this node.
599 : std::optional<tuples::tagged_tuple_from_typelist<FieldTagsList>>
600 : all_tensor_data{};
601 :
602 : // Retrieve the information needed to reconstruct which element the data
603 : // belongs to
604 : const auto source_grid_names = volume_file.get_grid_names(observation_id);
605 : const auto source_extents = volume_file.get_extents(observation_id);
606 : const auto source_bases = volume_file.get_bases(observation_id);
607 : const auto source_quadratures =
608 : volume_file.get_quadratures(observation_id);
609 : // Reconstruct domain from volume data file
610 : const std::optional<std::vector<char>> serialized_domain =
611 : volume_file.get_domain();
612 : if (serialized_domain.has_value()) {
613 : if (source_domain.has_value()) {
614 : #ifdef SPECTRE_DEBUG
615 : // Check that the domain is the same in all files (only in debug
616 : // mode)
617 : const auto deserialized_domain =
618 : deserialize<Domain<Dim>>(serialized_domain->data());
619 : if (*source_domain != deserialized_domain) {
620 : ERROR_NO_TRACE(
621 : "The domain in all volume files must be the same. Domain in "
622 : "file '"
623 : << file_name << volume_file.subfile_path()
624 : << "' differs from a previously read file.");
625 : }
626 : #endif
627 : } else {
628 : source_domain = deserialize<Domain<Dim>>(serialized_domain->data());
629 : }
630 : } else {
631 : Parallel::printf(
632 : "WARNING: No serialized domain found in file. "
633 : "Verification that elements in the source and target domain "
634 : "match will be skipped.\n");
635 : }
636 : // Reconstruct functions of time from volume data file
637 : if (source_domain_functions_of_time.empty() and
638 : source_domain.has_value() and
639 : alg::any_of(source_domain->blocks(), [](const auto& block) {
640 : return block.is_time_dependent();
641 : })) {
642 : const std::optional<std::vector<char>> serialized_functions_of_time =
643 : volume_file.get_functions_of_time(observation_id);
644 : if (not serialized_functions_of_time.has_value()) {
645 : ERROR_NO_TRACE("No domain functions of time found in file '"
646 : << file_name << volume_file.subfile_path()
647 : << "'. The functions of time are needed to verify the "
648 : "inertial coordinates with time-dependent maps.");
649 : }
650 : source_domain_functions_of_time =
651 : deserialize<domain::FunctionsOfTimeMap>(
652 : serialized_functions_of_time->data());
653 : }
654 :
655 : // Transfer the data to the target elements contained in this file. We
656 : // erase target elements when they are complete, so subsequent files only
657 : // search for the remaining elements and we can stop early.
658 : std::unordered_set<ElementId<Dim>> completed_target_elements{};
659 : for (const auto& target_element_id : target_element_ids) {
660 : const auto& [target_points, target_mesh] =
661 : get<Tags::RegisteredElements<Dim>>(box).at(
662 : Parallel::make_array_component_id<ReceiveComponent>(
663 : target_element_id));
664 : const auto target_grid_name = get_output(target_element_id);
665 : // Process this element only if it's in the file
666 : if (std::find(source_grid_names.begin(), source_grid_names.end(),
667 : target_grid_name) == source_grid_names.end()) {
668 : continue;
669 : }
670 :
671 : // Lazily load the tensor data from the file
672 : if (not all_tensor_data.has_value()) {
673 : all_tensor_data = detail::read_tensor_data<FieldTagsList>(
674 : volume_file, observation_id, selected_fields);
675 : }
676 :
677 : const auto source_mesh = h5::mesh_for_grid<Dim>(
678 : target_grid_name, source_grid_names, source_extents, source_bases,
679 : source_quadratures);
680 : const auto element_data_offset_and_length =
681 : h5::offset_and_length_for_grid(target_grid_name, source_grid_names,
682 : source_extents);
683 : auto source_element_data = detail::extract_element_data<FieldTagsList>(
684 : element_data_offset_and_length, *all_tensor_data, selected_fields);
685 :
686 : // Verify that the source and target elements really are the same
687 : if (source_domain.has_value()) {
688 : detail::verify_inertial_coordinates(*source_domain, observation_value,
689 : source_domain_functions_of_time,
690 : target_element_id, target_mesh,
691 : target_points);
692 : }
693 :
694 : // Transfer the data one-to-one, interpolating only if the meshes differ
695 : // by p-refinement
696 : tuples::tagged_tuple_from_typelist<FieldTagsList> target_element_data{};
697 : if (source_mesh == target_mesh) {
698 : target_element_data = std::move(source_element_data);
699 : } else {
700 : detail::interpolate_selected_fields<FieldTagsList>(
701 : make_not_null(&target_element_data), source_element_data,
702 : source_mesh, target_mesh, selected_fields);
703 : }
704 : if constexpr (Parallel::is_dg_element_collection_v<ReceiveComponent>) {
705 : ERROR("Can't yet do numerical initial data with nodegroups");
706 : } else {
707 : Parallel::receive_data<Tags::VolumeData<FieldTagsList>>(
708 : Parallel::get_parallel_component<ReceiveComponent>(
709 : cache)[target_element_id],
710 : volume_data_id, std::move(target_element_data));
711 : }
712 : completed_target_elements.insert(target_element_id);
713 : } // loop over registered elements
714 : for (const auto& completed_element_id : completed_target_elements) {
715 : target_element_ids.erase(completed_element_id);
716 : }
717 : // Stop early when all target elements are complete
718 : if (target_element_ids.empty()) {
719 : break;
720 : }
721 : } // loop over volume files
722 :
723 : // Have we completed all target elements? If we haven't, the source and
724 : // target domains probably don't match.
725 : if (not target_element_ids.empty()) {
726 : ERROR_NO_TRACE("The following "
727 : << target_element_ids.size()
728 : << " element(s) were not found in the source volume "
729 : "data files:\n"
730 : << target_element_ids
731 : << "\nMake sure the source and target domains match "
732 : "when 'ElementsAreIdentical' is enabled, or set "
733 : "it to 'False' to interpolate between the grids.");
734 : }
735 : }
736 :
737 : private:
738 : template <typename... LocalFieldTags>
739 : static tuples::TaggedTuple<Tags::Selected<LocalFieldTags>...>
740 0 : select_all_fields(tmpl::list<LocalFieldTags...> /*meta*/) {
741 : return {db::tag_name<LocalFieldTags>()...};
742 : }
743 : };
744 :
745 : } // namespace Actions
746 : } // namespace importers
|
