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             : 
      55             :   /*!
      56             :    * \brief Construct the interpolator without loading any volume data
      57             :    *
      58             :    * \param volume_files_or_glob The list of H5 files, or a glob pattern
      59             :    * \param subfile_name The name of the subfile in the H5 files containing the
      60             :    * volume data
      61             :    * \param tensor_components The tensor components to interpolate. The order of
      62             :    * the components in the vector is the order in which they will be returned.
      63             :    */
      64           1 :   SpacetimeInterpolator(
      65             :       std::variant<std::vector<std::string>, std::string> volume_files_or_glob,
      66             :       std::string subfile_name, std::vector<std::string> tensor_components);
      67             : 
      68             :   /// The maximum time bounds available in the volume data files (taking into
      69             :   /// account the number of ghost slices needed for interpolation). Only time
      70             :   /// bounds in this (inclusive) range can be requested with the
      71             :   /// `load_time_bounds()` function.
      72           1 :   std::array<double, 2> max_time_bounds() const;
      73             : 
      74             :   /*!
      75             :    * \brief Load volume data into memory such that the given time bounds are
      76             :    * covered.
      77             :    *
      78             :    * The given time bounds must be within the `max_time_bounds()` (inclusive).
      79             :    *
      80             :    * Previously loaded volume data slices outside the given time bounds will be
      81             :    * dropped and new ones loaded.
      82             :    *
      83             :    * \par Thread safety
      84             :    * This function is not generally thread safe, as it loads data from H5 files.
      85             :    * It should be called from a single thread only. It is safe to keep calling
      86             :    * interpolation routines from other threads while this function is running,
      87             :    * as long as they are requesting data within both the old and the new time
      88             :    * bounds.
      89             :    */
      90           1 :   void load_time_bounds(std::array<double, 2> time_bounds);
      91             : 
      92             :   /// The time bounds of the loaded data (inclusive). It is an error to call the
      93             :   /// interpolation routines outside these bounds.
      94           1 :   std::array<double, 2> time_bounds() const { return time_bounds_; }
      95             : 
      96             :   /// The number of loaded time slices.
      97           1 :   size_t num_loaded_slices() const { return interpolators_.size(); }
      98             : 
      99             :   /*!
     100             :    * \brief Interpolate to a single point
     101             :    *
     102             :    * This function is thread safe, so the interpolator can be constructed to
     103             :    * load the volume data once and then used by multiple threads to interpolate
     104             :    * to different points in parallel. Note that updating the `block_order` (if
     105             :    * provided) is not thread safe, so each thread should manage a separate
     106             :    * `block_order`.
     107             :    *
     108             :    * \param result the interpolated data at the target point. The vector is over
     109             :    * the number of components. Will be resized automatically.
     110             :    * \param target_point the point to interpolate to.
     111             :    * \param time the time to interpolate to. Must be within the `time_bounds()`.
     112             :    * \param block_order an optional priority order to search for the block
     113             :    * containing the target point. Will be updated when the point is found.
     114             :    * See `block_logical_coordinates_single_point` for more details.
     115             :    */
     116           1 :   void interpolate_to_point(gsl::not_null<std::vector<double>*> result,
     117             :                             const tnsr::I<double, Dim, Frame>& target_point,
     118             :                             double time,
     119             :                             std::optional<gsl::not_null<std::vector<size_t>*>>
     120             :                                 block_order = std::nullopt) const;
     121             : 
     122             :  private:
     123           0 :   std::variant<std::vector<std::string>, std::string> volume_files_or_glob_;
     124           0 :   std::string subfile_name_;
     125           0 :   std::vector<std::string> tensor_components_;
     126           0 :   constexpr static size_t time_interpolation_order_ = 3;  // must be odd
     127           0 :   constexpr static size_t num_ghost_slices_ =
     128             :       (time_interpolation_order_ + 1) / 2 - 1;
     129           0 :   domain::FunctionsOfTimeMap functions_of_time_;
     130             :   // Metadata collected from data files
     131           0 :   std::vector<size_t> all_observation_ids_;
     132           0 :   std::vector<double> all_observation_values_;
     133             :   // State
     134           0 :   std::array<double, 2> time_bounds_{
     135             :       {std::numeric_limits<double>::signaling_NaN(),
     136             :        std::numeric_limits<double>::signaling_NaN()}};
     137             :   // Loaded from data files
     138           0 :   std::vector<PointwiseInterpolator<Dim, ::Frame::Grid>> interpolators_{};
     139             :   // Mutex to protect other threads from accessing the `interpolators_` vector
     140             :   // while it is being modified
     141           0 :   mutable std::shared_mutex interpolators_mutex_;  // NOLINT(spectre-mutable)
     142             : };
     143             : 
     144             : }  // namespace spectre::Exporter