Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <optional>
8 : #include <shared_mutex>
9 : #include <string>
10 : #include <utility>
11 : #include <variant>
12 : #include <vector>
13 :
14 : #include "DataStructures/Tensor/TypeAliases.hpp"
15 : #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
16 : #include "IO/Exporter/PointwiseInterpolator.hpp"
17 : #include "Utilities/Gsl.hpp"
18 :
19 : namespace spectre::Exporter {
20 :
21 : /*!
22 : * \brief Interpolates data in volume files in both space and time
23 : *
24 : * This class reads the volume data at multiple time steps into memory and can
25 : * then interpolate it to any number of target points at any time within the
26 : * time bounds of the loaded data.
27 : *
28 : * \details The interpolation is first done in space and then in time. For the
29 : * interpolation in space, the PointwiseInterpolator class is used, which does
30 : * spectral interpolation within the elements. For the interpolation in time, a
31 : * polynomial interpolation with fixed order is used (currently third order,
32 : * meaning four time slices are used).
33 : *
34 : * \par Coordinate frames
35 : * The target points are given in the `Frame` specified as template parameter
36 : * (defaults to `Frame::Inertial`). However, the interpolation is done in the
37 : * grid frame, which is time-independent and therefore best suitable for
38 : * interpolation in time. Details of this idea can be found in
39 : * \cite Bohn:2016afc .
40 : *
41 : * \par Thread safety
42 : * Loading volume data from H5 files on multiple threads at once is not
43 : * generally thread safe, so `load_time_bounds` may only be called on a single
44 : * thread. However, once the data is loaded, the `interpolate_to_point` function
45 : * is thread safe. This means the volume data can be loaded once in a single
46 : * thread and then used by multiple threads to interpolate to different points
47 : * in parallel. Also, `load_time_bounds` can be called again to load new data
48 : * while other threads are still accessing the old data, as long as they only
49 : * request data within both the old and the new time bounds.
50 : */
51 : template <size_t Dim, typename Frame = ::Frame::Inertial>
52 1 : struct SpacetimeInterpolator {
53 0 : SpacetimeInterpolator() = default;
54 0 : SpacetimeInterpolator(const SpacetimeInterpolator&) = delete;
55 0 : SpacetimeInterpolator& operator=(const SpacetimeInterpolator&) = delete;
56 0 : SpacetimeInterpolator(SpacetimeInterpolator&&);
57 0 : SpacetimeInterpolator& operator=(SpacetimeInterpolator&&);
58 0 : ~SpacetimeInterpolator() = default;
59 :
60 : /*!
61 : * \brief Construct the interpolator without loading any volume data
62 : *
63 : * \param volume_files_or_glob The list of H5 files, or a glob pattern
64 : * \param subfile_name The name of the subfile in the H5 files containing the
65 : * volume data
66 : * \param tensor_components The tensor components to interpolate. The order of
67 : * the components in the vector is the order in which they will be returned.
68 : */
69 1 : SpacetimeInterpolator(
70 : std::variant<std::vector<std::string>, std::string> volume_files_or_glob,
71 : std::string subfile_name, std::vector<std::string> tensor_components);
72 :
73 : /// The maximum time bounds available in the volume data files (taking into
74 : /// account the number of ghost slices needed for interpolation). Only time
75 : /// bounds in this (inclusive) range can be requested with the
76 : /// `load_time_bounds()` function.
77 1 : std::array<double, 2> max_time_bounds() const;
78 :
79 : /*!
80 : * \brief Load volume data into memory such that the given time bounds are
81 : * covered.
82 : *
83 : * The given time bounds must be within the `max_time_bounds()` (inclusive).
84 : *
85 : * Previously loaded volume data slices outside the given time bounds will be
86 : * dropped and new ones loaded.
87 : *
88 : * \par Thread safety
89 : * This function is not generally thread safe, as it loads data from H5 files.
90 : * It should be called from a single thread only. It is safe to keep calling
91 : * interpolation routines from other threads while this function is running,
92 : * as long as they are requesting data within both the old and the new time
93 : * bounds.
94 : */
95 1 : void load_time_bounds(std::array<double, 2> time_bounds);
96 :
97 : /// The time bounds of the loaded data (inclusive). It is an error to call the
98 : /// interpolation routines outside these bounds.
99 1 : std::array<double, 2> time_bounds() const { return time_bounds_; }
100 :
101 : /// The number of loaded time slices.
102 1 : size_t num_loaded_slices() const { return interpolators_.size(); }
103 :
104 : /*!
105 : * \brief Interpolate to a single point
106 : *
107 : * This function is thread safe, so the interpolator can be constructed to
108 : * load the volume data once and then used by multiple threads to interpolate
109 : * to different points in parallel. Note that updating the `block_order` (if
110 : * provided) is not thread safe, so each thread should manage a separate
111 : * `block_order`.
112 : *
113 : * \param result the interpolated data at the target point. The vector is over
114 : * the number of components. Will be resized automatically.
115 : * \param target_point the point to interpolate to.
116 : * \param time the time to interpolate to. Must be within the `time_bounds()`.
117 : * \param block_order an optional priority order to search for the block
118 : * containing the target point. Will be updated when the point is found.
119 : * See `block_logical_coordinates_single_point` for more details.
120 : */
121 1 : void interpolate_to_point(gsl::not_null<std::vector<double>*> result,
122 : const tnsr::I<double, Dim, Frame>& target_point,
123 : double time,
124 : std::optional<gsl::not_null<std::vector<size_t>*>>
125 : block_order = std::nullopt) const;
126 :
127 : private:
128 0 : std::variant<std::vector<std::string>, std::string> volume_files_or_glob_;
129 0 : std::string subfile_name_;
130 0 : std::vector<std::string> tensor_components_;
131 0 : constexpr static size_t time_interpolation_order_ = 3; // must be odd
132 0 : constexpr static size_t num_ghost_slices_ =
133 : (time_interpolation_order_ + 1) / 2 - 1;
134 0 : domain::FunctionsOfTimeMap functions_of_time_;
135 : // Metadata collected from data files
136 0 : std::vector<size_t> all_observation_ids_;
137 0 : std::vector<double> all_observation_values_;
138 : // State
139 0 : std::array<double, 2> time_bounds_{
140 : {std::numeric_limits<double>::signaling_NaN(),
141 : std::numeric_limits<double>::signaling_NaN()}};
142 : // Loaded from data files
143 0 : std::vector<PointwiseInterpolator<Dim, ::Frame::Grid>> interpolators_{};
144 : // Mutex to protect other threads from accessing the `interpolators_` vector
145 : // while it is being modified
146 0 : mutable std::shared_mutex interpolators_mutex_; // NOLINT(spectre-mutable)
147 : };
148 :
149 : } // namespace spectre::Exporter
|
