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