Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <Exporter.hpp>  // The SpEC Exporter
       7             : #include <cstddef>
       8             : #include <memory>
       9             : #include <utility>
      10             : #include <vector>
      11             : 
      12             : #include "DataStructures/Tensor/Tensor.hpp"
      13             : #include "Utilities/ContainerHelpers.hpp"
      14             : #include "Utilities/ErrorHandling/Assert.hpp"
      15             : #include "Utilities/SetNumberOfGridPoints.hpp"
      16             : #include "Utilities/TMPL.hpp"
      17             : #include "Utilities/TaggedTuple.hpp"
      18             : 
      19             : namespace io {
      20             : 
      21             : /*!
      22             :  * \brief Interpolate numerical SpEC initial data to arbitrary points
      23             :  *
      24             :  * \tparam Tags List of tags to load. The tags must correspond exactly to the
      25             :  * list of variables that the `spec_exporter` was configured with. This function
      26             :  * does not support interpolating only a subset of the variables. This is a
      27             :  * limitation of the `spec::Exporter`.
      28             :  * \tparam DataType `double` or `DataVector`.
      29             :  * \tparam CoordFrame The frame of the coordinates `x`. These coordinates are
      30             :  * always assumed to be in SpEC's "grid" frame.
      31             :  *
      32             :  * \param spec_exporter Has the numerical data already loaded. It must be
      33             :  * configured with the same number of variables that were passed as `Tags`, and
      34             :  * in the same order.
      35             :  * \param x Interpolate to these coordinates. They are assumed to be in SpEC's
      36             :  * "grid" frame.
      37             :  * \param which_interpolator Index of the interpolator. See `spec::Exporter`
      38             :  * documentation for details.
      39             :  */
      40             : template <typename Tags, typename DataType, typename CoordFrame>
      41           1 : tuples::tagged_tuple_from_typelist<Tags> interpolate_from_spec(
      42             :     const gsl::not_null<spec::Exporter*> spec_exporter,
      43             :     const tnsr::I<DataType, 3, CoordFrame>& x,
      44             :     const size_t which_interpolator) {
      45             :   // The `spec::Exporter` doesn't currently expose its `vars_to_interpolate`.
      46             :   // Once it does, we can assert the number of variables here.
      47             :   // const auto& dataset_names = spec_exporter.vars_to_interpolate();
      48             :   // ASSERT(tmpl::size<Tags>::value == dataset_names.size(),
      49             :   //        "Mismatch between number of tags and dataset names. The SpEC
      50             :   //        exporter " "was configured with " +
      51             :   //            std::to_string(dataset_names.size()) +
      52             :   //            " variables to interpolate, but requested " +
      53             :   //            std::to_string(tmpl::size<Tags>::value) + " tags.");
      54             :   // Transform coordinates into SpEC's expected format. SpEC expects grid
      55             :   // coordinates. We assume that the grid coordinates coincide with inertial
      56             :   // coordinates for the initial data.
      57             :   const size_t num_points = get_size(get<0>(x));
      58             :   std::vector<std::vector<double>> spec_grid_coords(num_points);
      59             :   for (size_t i = 0; i < num_points; ++i) {
      60             :     spec_grid_coords[i] = std::vector<double>{get_element(get<0>(x), i),
      61             :                                               get_element(get<1>(x), i),
      62             :                                               get_element(get<2>(x), i)};
      63             :   }
      64             :   // Allocate memory for the interpolation and point into it
      65             :   tuples::tagged_tuple_from_typelist<Tags> interpolation_buffer{};
      66             :   std::vector<std::vector<double*>> buffer_pointers(tmpl::size<Tags>::value);
      67             :   size_t var_i = 0;
      68             :   tmpl::for_each<Tags>([&interpolation_buffer, &buffer_pointers, &num_points,
      69             :                         &var_i](const auto tag_v) {
      70             :     using tag = tmpl::type_from<std::decay_t<decltype(tag_v)>>;
      71             :     auto& tensor = get<tag>(interpolation_buffer);
      72             :     set_number_of_grid_points(make_not_null(&tensor), num_points);
      73             :     const size_t num_components = tensor.size();
      74             :     buffer_pointers[var_i].resize(num_components);
      75             :     // The SpEC exporter supports tensors up to symmetric rank 2, which are
      76             :     // ordered xx, yx, zx, yy, zy, zz. Because this is also the order in
      77             :     // which we store components in the Tensor class, we don't have to do
      78             :     // anything special here.
      79             :     // WARNING: If the Tensor storage order changes for some reason, this
      80             :     // code needs to be updated.
      81             :     for (size_t component_i = 0; component_i < num_components; ++component_i) {
      82             : #ifdef SPECTRE_DEBUG
      83             :       const auto component_name =
      84             :           tensor.component_name(tensor.get_tensor_index(component_i));
      85             :       if constexpr (std::decay_t<decltype(tensor)>::rank() == 1) {
      86             :         const std::array<std::string, 3> spec_component_order{{"x", "y", "z"}};
      87             :         ASSERT(component_name == gsl::at(spec_component_order, component_i),
      88             :                "Unexpected Tensor component order. Expected "
      89             :                    << spec_component_order);
      90             :       } else if constexpr (std::decay_t<decltype(tensor)>::rank() == 2 and
      91             :                            std::decay_t<decltype(tensor)>::size() == 6) {
      92             :         const std::array<std::string, 6> spec_component_order{
      93             :             {"xx", "yx", "zx", "yy", "zy", "zz"}};
      94             :         ASSERT(component_name == gsl::at(spec_component_order, component_i),
      95             :                "Unexpected Tensor component order. Expected "
      96             :                    << spec_component_order);
      97             :       } else if constexpr (std::decay_t<decltype(tensor)>::rank() > 0) {
      98             :         ASSERT(
      99             :             false,
     100             :             "Unsupported Tensor type for SpEC import. Only scalars, "
     101             :             "vectors, and symmetric rank-2 tensors are currently supported.");
     102             :       }
     103             : #endif
     104             :       auto& component = tensor[component_i];
     105             :       auto& component_pointer = buffer_pointers[var_i][component_i];
     106             :       if constexpr (std::is_same_v<DataType, double>) {
     107             :         component_pointer = &component;
     108             :       } else {
     109             :         component_pointer = component.data();
     110             :       }
     111             :     }
     112             :     ++var_i;
     113             :   });
     114             :   // Interpolate!
     115             :   spec_exporter->interpolate(buffer_pointers, spec_grid_coords,
     116             :                              which_interpolator);
     117             :   return interpolation_buffer;
     118             : }
     119             : 
     120             : }  // namespace io