Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cstddef>
8 : #include <map>
9 : #include <optional>
10 : #include <string>
11 : #include <unordered_map>
12 : #include <variant>
13 : #include <vector>
14 :
15 : #include "DataStructures/Tensor/TypeAliases.hpp"
16 : #include "Domain/Domain.hpp"
17 : #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
18 : #include "Domain/Structure/ElementId.hpp"
19 : #include "Domain/Structure/ElementSearchTree.hpp"
20 : #include "IO/Logging/Verbosity.hpp"
21 : #include "NumericalAlgorithms/Interpolation/UniformCardinalBSpline.hpp"
22 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
23 : #include "Utilities/Gsl.hpp"
24 :
25 : namespace spectre::Exporter {
26 :
27 : /*!
28 : * \brief Interpolate tensor components in both space and time using modal
29 : * data. This is much more efficient and accurate than interpolating nodal
30 : * data.
31 : *
32 : * \details This class builds a time interpolant for every modal coefficient of
33 : * every tensor component on each element. It reads volume data one element at
34 : * a time, converts nodal values to modal coefficients, and constructs a
35 : * cardinal cubic B-spline in time for each mode. At evaluation time it locates
36 : * the element containing the target point, evaluates the modal interpolants at
37 : * the requested time, and evaluates the Legendre series at the element-logical
38 : * point.
39 : *
40 : * ### How to use
41 : *
42 : * In a simulation, create a few events that write out volume data at different
43 : * time intervals using ObserveFields. The observation times need to be uniform
44 : * within each subfile, so either observe at a fixed slab interval if the slab
45 : * size is constant, or at fixed times using a dense trigger. The idea is that
46 : * there is a very coarse grid written very frequently and a finer grid written
47 : * less frequently. This way, the lowest modes are resolved very accurately and
48 : * the higher modes are approximated with a larger error. A useful choice in
49 : * practice is:
50 : *
51 : * - Observe the lowest 3-4 modes very frequently, e.g. every 0.5M or even more
52 : * often.
53 : * - Observe the lowest 6-7 modes 10 times less frequently.
54 : * - Observe all modes 10 times less frequently again.
55 : *
56 : * Important: use the `ProjectToMesh` option of ObserveFields to ensure that
57 : * the data is truncated cleanly in modal space. Do not use `InterpolateToMesh`
58 : * as this creates a catastrophic aliasing error. Also, do not use the
59 : * `CleanFunctionsOfTime` event, as this removes the required global history of
60 : * the functions of time.
61 : *
62 : * Then construct the ModalSpacetimeInterpolator from these subfiles in
63 : * priority order, from the coarse mesh observed most frequently to the fine
64 : * mesh that defines the final spatial resolution. The first subfile containing
65 : * a mode owns it, and later subfiles only supply the modes that are missing on
66 : * the earlier meshes. The subfile meshes must therefore be nested, i.e. have
67 : * component-wise nondecreasing extents, and all subfiles must contain the same
68 : * elements and the same domain. The subfiles may have different time steps and
69 : * nonaligned observation times; the available time interval is their
70 : * intersection.
71 : *
72 : * ### Error tolerance and compression
73 : *
74 : * The time series of every owned mode is compressed into an
75 : * `intrp::UniformCardinalBSpline` on the coarsest uniform time grid that
76 : * reproduces the samples to within an absolute error tolerance, and modes
77 : * whose amplitude never exceeds the tolerance are dropped entirely since they
78 : * are not important for reconstructing the field (this can be up to 80% of the
79 : * modes in practice). This way, the interpolator only keeps the data that is
80 : * necessary.
81 : *
82 : * By default, the tolerance for each element and tensor component is estimated
83 : * from the lowest mode in the first subfile: it is the interpolation error of
84 : * that mode when using all of its observed data, so all other modes are
85 : * interpolated with a conservatively similar error. Note that this infers the
86 : * tolerance from the temporal variability of the data rather than from an
87 : * accuracy requirement, so specify the tolerance explicitly when you know the
88 : * accuracy you need. The estimate requires at least 9 observations in the
89 : * first subfile; an explicit tolerance does not. Either way this is a per-mode
90 : * compression heuristic. It does not bound the total reconstructed field error
91 : * and does not detect temporal undersampling in the original volume data.
92 : *
93 : * ### Limitations
94 : *
95 : * Only the Legendre basis is supported for now. Changes to the element
96 : * topology or per-element mesh within one time series (e.g. AMR), restarts
97 : * that produce non-uniform observation times, and element migration between
98 : * files within one time series are not supported. Files from different
99 : * simulation segments must be joined first, with overlapping observations
100 : * removed. Time-dependent domains use the global functions of time stored
101 : * with the final subfile, and all functions required by the domain must
102 : * cover the full available time interval.
103 : *
104 : * Since the time grid of each mode is uniform, it is a good idea to use a
105 : * separate interpolator for parts of the evolution that are dynamically
106 : * different (e.g. junk radiation vs. inspiral vs. merger vs. ringdown). A
107 : * wrapper that patches together several interpolators may be added in the
108 : * future.
109 : *
110 : * \note Due to a bug in the boost cardinal cubic B-spline implementation that
111 : * was fixed only recently, see
112 : * https://github.com/boostorg/math/commit/4809e714d4806c07da3a3def0c4550daa0529b8d,
113 : * the interpolator can only be built with boost version 1.81 or later. See
114 : * `intrp::UniformCardinalBSpline` for details.
115 : */
116 : template <size_t Dim, typename Frame = ::Frame::Inertial>
117 1 : class ModalSpacetimeInterpolator {
118 : public:
119 : /*!
120 : * \brief Construct from one or more volume files and nested volume
121 : * subfiles.
122 : *
123 : * \param volume_files_or_glob A list of volume H5 files, or a glob string
124 : * that resolves to volume files. The files may distribute the elements
125 : * of each observation across nodes.
126 : * \param subfiles_in_priority_order Volume subfile names ordered from the
127 : * preferred source of low modes to the source of the final spatial
128 : * mesh. The meshes must have component-wise nondecreasing extents.
129 : * \param tensor_components Tensor component names to interpolate. The
130 : * returned values have this order.
131 : * \param start_time Optional lower bound used to select observations from
132 : * every subfile.
133 : * \param end_time Optional upper bound used to select observations from
134 : * every subfile.
135 : * \param error_tolerance Optional absolute tolerance used to drop and
136 : * compress modes. By default the tolerance is estimated from the
137 : * lowest mode in the first subfile, which then needs at least 9
138 : * observations.
139 : * \param verbosity Controls diagnostic output during construction.
140 : */
141 1 : ModalSpacetimeInterpolator(
142 : const std::variant<std::vector<std::string>, std::string>&
143 : volume_files_or_glob,
144 : std::vector<std::string> subfiles_in_priority_order,
145 : std::vector<std::string> tensor_components,
146 : std::optional<double> start_time = std::nullopt,
147 : std::optional<double> end_time = std::nullopt,
148 : std::optional<double> error_tolerance = std::nullopt,
149 : Verbosity verbosity = Verbosity::Quiet);
150 :
151 0 : ModalSpacetimeInterpolator(ModalSpacetimeInterpolator&&) = default;
152 0 : ModalSpacetimeInterpolator& operator=(ModalSpacetimeInterpolator&&) = default;
153 0 : ModalSpacetimeInterpolator(const ModalSpacetimeInterpolator&) = delete;
154 0 : ModalSpacetimeInterpolator& operator=(const ModalSpacetimeInterpolator&) =
155 : delete;
156 0 : ~ModalSpacetimeInterpolator() = default;
157 :
158 : /*!
159 : * \brief Interpolate the tensor components at a spacetime point.
160 : *
161 : * \param result Output buffer, resized to the number of tensor components.
162 : * \param target_point Coordinates of the target point in `Frame`.
163 : * \param time Time at which to evaluate the interpolant.
164 : * \param block_order Optional block-ordering hint that speeds up repeated
165 : * block searches.
166 : */
167 1 : void interpolate_to_point(gsl::not_null<std::vector<double>*> result,
168 : const tnsr::I<double, Dim, Frame>& target_point,
169 : double time,
170 : std::optional<gsl::not_null<std::vector<size_t>*>>
171 : block_order = std::nullopt) const;
172 :
173 : /// Tensor components in the order returned by `interpolate_to_point()`
174 1 : const std::vector<std::string>& tensor_components() const {
175 : return tensor_components_;
176 : }
177 :
178 : /// Inclusive time interval shared by all input subfiles
179 1 : const std::array<double, 2>& time_bounds() const { return time_bounds_; }
180 :
181 : private:
182 0 : using ModeInterpolant = std::optional<intrp::UniformCardinalBSpline>;
183 :
184 0 : struct ElementData {
185 0 : Mesh<Dim> mesh{};
186 : // Indexed by [component][mode on the final mesh]
187 0 : std::vector<std::vector<ModeInterpolant>> interpolants;
188 : };
189 :
190 0 : std::vector<std::string> tensor_components_;
191 0 : std::array<double, 2> time_bounds_{};
192 0 : Domain<Dim> domain_{};
193 0 : domain::FunctionsOfTimeMap functions_of_time_;
194 0 : std::map<size_t, domain::ElementSearchTree<Dim>> element_search_trees_{};
195 0 : std::unordered_map<ElementId<Dim>, ElementData> element_data_{};
196 : };
197 :
198 : } // namespace spectre::Exporter
|
