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 <optional>
       9             : #include <string>
      10             : #include <unordered_map>
      11             : #include <variant>
      12             : #include <vector>
      13             : 
      14             : #include "DataStructures/DataBox/TagName.hpp"
      15             : #include "DataStructures/Tensor/Tensor.hpp"
      16             : #include "Domain/BlockLogicalCoordinates.hpp"
      17             : #include "Domain/Domain.hpp"
      18             : #include "Domain/Structure/ElementSearchTree.hpp"
      19             : #include "IO/Exporter/Exporter.hpp"
      20             : #include "Utilities/TaggedTuple.hpp"
      21             : 
      22             : namespace spectre::Exporter {
      23             : 
      24             : /// Collect tensor component names from Tags list
      25             : template <typename Tags>
      26           1 : auto get_tensor_components() {
      27             :   std::vector<std::string> tensor_components{};
      28             :   tmpl::for_each<Tags>([&tensor_components](auto tag_v) {
      29             :     using tensor_tag = tmpl::type_from<decltype(tag_v)>;
      30             :     using TensorType = typename tensor_tag::type;
      31             :     const std::string tag_name = db::tag_name<tensor_tag>();
      32             :     for (size_t i = 0; i < TensorType::size(); ++i) {
      33             :       const std::string component_name =
      34             :           tag_name + TensorType::component_suffix(i);
      35             :       tensor_components.push_back(component_name);
      36             :     }
      37             :   });
      38             :   return tensor_components;
      39             : }
      40             : 
      41             : /// Convert tensor components to a tagged_tuple
      42             : template <typename Tags, typename DataType>
      43           1 : auto make_tagged_tuple(std::vector<DataType> interpolated_data) {
      44             :   tuples::tagged_tuple_from_typelist<Tags> result{};
      45             :   size_t component_index = 0;
      46             :   tmpl::for_each<Tags>(
      47             :       [&component_index, &interpolated_data, &result](auto tag_v) {
      48             :         using tensor_tag = tmpl::type_from<decltype(tag_v)>;
      49             :         using TensorType = typename tensor_tag::type;
      50             :         for (size_t i = 0; i < TensorType::size(); ++i) {
      51             :           get<tensor_tag>(result)[i] =
      52             :               std::move(interpolated_data[component_index]);
      53             :           ++component_index;
      54             :         }
      55             :       });
      56             :   return result;
      57             : }
      58             : 
      59             : /// @{
      60             : /*!
      61             :  * \brief Interpolates data in volume files to target points
      62             :  *
      63             :  * These are overloads of the `interpolate_to_points` function that work with
      64             :  * Tensor types and tags, rather than the raw C++ types that are used in
      65             :  * Exporter.hpp so it can be used by external programs.
      66             :  *
      67             :  * The `Tags` template parameter is a typelist of tags that should be read from
      68             :  * the volume files. The dataset names to read are constructed from the tag
      69             :  * names. Here is an example of how to use this function:
      70             :  *
      71             :  * \snippet Test_Exporter.cpp interpolate_tensors_to_points_example
      72             :  */
      73             : template <typename ResultDataType, size_t Dim, typename Frame>
      74           1 : void interpolate_to_points(
      75             :     gsl::not_null<std::vector<ResultDataType>*> result,
      76             :     const std::variant<std::vector<std::string>, std::string>&
      77             :         volume_files_or_glob,
      78             :     const std::string& subfile_name, const ObservationVariant& observation,
      79             :     const std::vector<std::string>& tensor_components,
      80             :     const tnsr::I<DataVector, Dim, Frame>& target_points,
      81             :     bool extrapolate_into_excisions = false,
      82             :     bool error_on_missing_points = false,
      83             :     std::optional<size_t> num_threads = std::nullopt);
      84             : 
      85             : template <size_t Dim, typename Frame>
      86           1 : std::vector<DataVector> interpolate_to_points(
      87             :     const std::variant<std::vector<std::string>, std::string>&
      88             :         volume_files_or_glob,
      89             :     const std::string& subfile_name, const ObservationVariant& observation,
      90             :     const std::vector<std::string>& tensor_components,
      91             :     const tnsr::I<DataVector, Dim, Frame>& target_points,
      92             :     bool extrapolate_into_excisions = false,
      93             :     bool error_on_missing_points = false,
      94             :     std::optional<size_t> num_threads = std::nullopt) {
      95             :   std::vector<DataVector> interpolated_data{};
      96             :   interpolate_to_points(make_not_null(&interpolated_data), volume_files_or_glob,
      97             :                         subfile_name, observation, tensor_components,
      98             :                         target_points, extrapolate_into_excisions,
      99             :                         error_on_missing_points, num_threads);
     100             :   return interpolated_data;
     101             : }
     102             : 
     103             : template <typename Tags, size_t Dim, typename Frame>
     104           1 : tuples::tagged_tuple_from_typelist<Tags> interpolate_to_points(
     105             :     const std::variant<std::vector<std::string>, std::string>&
     106             :         volume_files_or_glob,
     107             :     const std::string& subfile_name, const ObservationVariant& observation,
     108             :     const tnsr::I<DataVector, Dim, Frame>& target_points,
     109             :     bool extrapolate_into_excisions = false,
     110             :     const bool error_on_missing_points = false,
     111             :     std::optional<size_t> num_threads = std::nullopt) {
     112             :   return make_tagged_tuple<Tags>(interpolate_to_points(
     113             :       volume_files_or_glob, subfile_name, observation,
     114             :       get_tensor_components<Tags>(), target_points, extrapolate_into_excisions,
     115             :       error_on_missing_points, num_threads));
     116             : }
     117             : /// @}
     118             : 
     119             : /*!
     120             :  * \brief Interpolates data in volume files to target points by reading the
     121             :  * volume data into memory
     122             :  *
     123             :  * This class reads the volume data at the requested time into memory and can
     124             :  * then interpolate it to any number of target points.
     125             :  *
     126             :  * \par Thread safety
     127             :  * Constructing the `PointwiseInterpolator` on multiple threads at once is not
     128             :  * thread safe (unless HDF5 was built with thread-safety support) because the
     129             :  * constructor opens the H5 files and reads in the data. However, once the data
     130             :  * is loaded, the `interpolate_to_points` and `interpolate_to_point` functions
     131             :  * are thread safe. This means the volume data can be loaded once in a single
     132             :  * thread and then used by multiple threads to interpolate to different points
     133             :  * in parallel.
     134             :  */
     135             : template <size_t Dim, typename Frame = ::Frame::Inertial>
     136           1 : struct PointwiseInterpolator {
     137           0 :   PointwiseInterpolator() = default;
     138           0 :   PointwiseInterpolator(const std::variant<std::vector<std::string>,
     139             :                                            std::string>& volume_files_or_glob,
     140             :                         const std::string& subfile_name,
     141             :                         const ObservationVariant& observation,
     142             :                         const std::vector<std::string>& tensor_components);
     143             : 
     144             :   /*!
     145             :    * \brief Interpolate to many points
     146             :    *
     147             :    * \param result the interpolated data at the target points. The outer vector
     148             :    * is the number of components and the inner vector is the number of target
     149             :    * points. Will be resized automatically.
     150             :    * \param target_points the points to interpolate to
     151             :    * \param extrapolate_into_excisions whether to extrapolate into excised
     152             :    * regions
     153             :    * \param error_on_missing_points whether to throw an error if any of the
     154             :    * target points are outside the domain
     155             :    * \param num_threads the number of OpenMP threads to use to parallelize the
     156             :    * interpolation over the target points. If not provided, the default number
     157             :    * of threads will be used.
     158             :    */
     159           1 :   void interpolate_to_points(
     160             :       gsl::not_null<std::vector<DataVector>*> result,
     161             :       const tnsr::I<DataVector, Dim, Frame>& target_points,
     162             :       bool extrapolate_into_excisions = false,
     163             :       bool error_on_missing_points = false,
     164             :       std::optional<size_t> num_threads = std::nullopt) const;
     165             : 
     166             :   /*!
     167             :    * \brief Interpolate to a single point
     168             :    *
     169             :    * This function is thread safe, so the interpolator can be constructed to
     170             :    * load the volume data once and then used by multiple threads to interpolate
     171             :    * to different points in parallel. Note that updating the `block_order` (if
     172             :    * provided) is not thread safe, so each thread should manage a separate
     173             :    * `block_order`.
     174             :    *
     175             :    * \param result the interpolated data at the target point. The vector is over
     176             :    * the number of components. Will be resized automatically.
     177             :    * \param target_point the point to interpolate to
     178             :    * \param block_order an optional priority order to search for the block
     179             :    * containing the target point. Will be updated when the point is found.
     180             :    * See `block_logical_coordinates_single_point` for more details.
     181             :    */
     182           1 :   void interpolate_to_point(gsl::not_null<std::vector<double>*> result,
     183             :                             const tnsr::I<double, Dim, Frame>& target_point,
     184             :                             std::optional<gsl::not_null<std::vector<size_t>*>>
     185             :                                 block_order = std::nullopt) const;
     186             : 
     187             :   /*!
     188             :    * \brief Interpolate to a single point in block-logical coordinates
     189             :    *
     190             :    * \param result the interpolated data at the target point. The vector is over
     191             :    * the number of components. Will be resized automatically.
     192             :    * \param target_point the point to interpolate to in block-logical
     193             :    * coordinates.
     194             :    */
     195           1 :   void interpolate_to_point(
     196             :       gsl::not_null<std::vector<double>*> result,
     197             :       const IdPair<domain::BlockId,
     198             :                    tnsr::I<double, Dim, ::Frame::BlockLogical>>& target_point)
     199             :       const;
     200             : 
     201           0 :   size_t obs_id() const { return obs_id_; }
     202           0 :   double time() const { return time_; }
     203           0 :   const Domain<Dim>& domain() const { return domain_; }
     204           0 :   const domain::FunctionsOfTimeMap& functions_of_time() const {
     205             :     return functions_of_time_;
     206             :   }
     207             : 
     208             :  private:
     209             :   // Loaded from data files
     210           0 :   size_t obs_id_ = std::numeric_limits<size_t>::max();
     211           0 :   double time_ = std::numeric_limits<double>::signaling_NaN();
     212           0 :   Domain<Dim> domain_;
     213           0 :   domain::FunctionsOfTimeMap functions_of_time_;
     214             :   // Outer vector is the source data file
     215           0 :   std::vector<std::vector<ElementId<Dim>>> element_ids_;
     216             :   std::vector<std::map<size_t, domain::ElementSearchTree<Dim>>>
     217           0 :       element_search_trees_;
     218             :   std::vector<
     219             :       std::unordered_map<ElementId<Dim>, std::tuple<Mesh<Dim>, size_t, size_t>>>
     220           0 :       meshes_;
     221             :   std::vector<std::vector<std::variant<DataVector, std::vector<float>>>>
     222           0 :       tensor_data_;
     223             : };
     224             : 
     225             : }  // namespace spectre::Exporter