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