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