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