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/TaggedTuple.hpp"
13 : #include "DataStructures/Tensor/Tensor.hpp"
14 : #include "Utilities/ContainerHelpers.hpp"
15 : #include "Utilities/ErrorHandling/Assert.hpp"
16 : #include "Utilities/SetNumberOfGridPoints.hpp"
17 : #include "Utilities/TMPL.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
|
