Line data Source code
1 1 : // Distributed under the MIT License. 2 : // See LICENSE.txt for details. 3 : 4 : /// \file 5 : /// Declares functions for interpolating data in volume files to target points. 6 : /// This file is intended to be included in external programs, so it 7 : /// intentionally has no dependencies on any other headers. 8 : 9 : #pragma once 10 : 11 : #include <array> 12 : #include <cstddef> 13 : #include <optional> 14 : #include <string> 15 : #include <variant> 16 : #include <vector> 17 : 18 : /// Functions that are intended to be used by external programs, e.g. to 19 : /// interpolate data in volume files to target points. 20 0 : namespace spectre::Exporter { 21 : 22 : /// Identifies an observation by its ID in the volume data file. 23 1 : struct ObservationId { 24 0 : ObservationId() = default; 25 0 : explicit ObservationId(size_t local_value) : value(local_value) {} 26 0 : size_t value; 27 : }; 28 : 29 : /// Identifies an observation by its index in the ordered list of observations. 30 : /// Negative indices are counted from the end of the list. 31 1 : struct ObservationStep { 32 0 : ObservationStep() = default; 33 0 : explicit ObservationStep(int local_value) : value(local_value) {} 34 0 : int value; 35 : }; 36 : 37 0 : using ObservationVariant = std::variant<ObservationId, ObservationStep, double>; 38 : 39 : /*! 40 : * \brief Interpolate data in volume files to target points 41 : * 42 : * \tparam Dim Dimension of the domain 43 : * \param volume_files_or_glob The list of H5 files, or a glob pattern 44 : * \param subfile_name The name of the subfile in the H5 files containing the 45 : * volume data 46 : * \param observation Either the observation ID as a `size_t`, or the index of 47 : * the observation in the volume files to interpolate as an `int` (a value of 0 48 : * would be the first observation, and a value of -1 would be the last 49 : * observation). 50 : * \param tensor_components The tensor components to interpolate, e.g. 51 : * "Lapse", "Shift_x", "Shift_y", "Shift_z", "SpatialMetric_xx", etc. 52 : * Look into the H5 file to see what components are available. 53 : * \param target_points The points to interpolate to, in inertial coordinates. 54 : * \param extrapolate_into_excisions Enables extrapolation into excision regions 55 : * of the domain (default is `false`). This can be useful to fill the excision 56 : * region with (constraint-violating but smooth) data so it can be imported into 57 : * moving puncture codes. Specifically, we implement the strategy used in 58 : * \cite Etienne2008re adjusted for distorted excisions: we choose uniformly 59 : * spaced radial anchor points spaced as $\Delta r = 0.3 r_\mathrm{AH}$ in the 60 : * grid frame (where the excision is spherical), then map the anchor points to 61 : * the distorted frame (where we have the target point) and do a 7th order 62 : * polynomial extrapolation into the excision region. 63 : * \param error_on_missing_points If `true`, an error will be thrown if any of 64 : * the target points are outside the domain. If `false`, the result will be 65 : * filled with NaNs for points outside the domain (default is `false`). 66 : * \param num_threads The number of threads to use if OpenMP is linked in. If 67 : * not specified, OpenMP will determine the number of threads automatically. 68 : * It's also possible to set the number of threads using the environment 69 : * variable OMP_NUM_THREADS. It's an error to specify num_threads if OpenMP is 70 : * not linked in. Set num_threads to 1 to disable OpenMP. 71 : * \return std::vector<std::vector<double>> The interpolated data. The first 72 : * dimension corresponds to the selected tensor components, and the second 73 : * dimension corresponds to the target points. 74 : */ 75 : template <size_t Dim> 76 1 : std::vector<std::vector<double>> interpolate_to_points( 77 : const std::variant<std::vector<std::string>, std::string>& 78 : volume_files_or_glob, 79 : const std::string& subfile_name, const ObservationVariant& observation, 80 : const std::vector<std::string>& tensor_components, 81 : const std::array<std::vector<double>, Dim>& target_points, 82 : bool extrapolate_into_excisions = false, 83 : bool error_on_missing_points = false, 84 : std::optional<size_t> num_threads = std::nullopt); 85 : 86 : } // namespace spectre::Exporter