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 <type_traits>
8 :
9 : #include "DataStructures/ApplyMatrices.hpp"
10 : #include "DataStructures/DataVector.hpp"
11 : #include "DataStructures/Matrix.hpp"
12 : #include "DataStructures/SliceIterator.hpp"
13 : #include "DataStructures/Variables.hpp"
14 : #include "Domain/Structure/Direction.hpp"
15 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
16 : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
17 : #include "NumericalAlgorithms/Spectral/Spectral.hpp"
18 : #include "Utilities/ErrorHandling/Assert.hpp"
19 : #include "Utilities/Gsl.hpp"
20 : #include "Utilities/MakeArray.hpp"
21 : #include "Utilities/TMPL.hpp"
22 :
23 : namespace dg {
24 : namespace detail {
25 : template <typename TagsList>
26 : struct NumberOfIndependentComponents;
27 :
28 : template <typename... Tags>
29 : struct NumberOfIndependentComponents<tmpl::list<Tags...>> {
30 : static constexpr size_t value = (... + Tags::type::size());
31 : };
32 : } // namespace detail
33 :
34 : /*!
35 : * \brief Projects a `Variables` of volume data to a contiguous subset of
36 : * a boundary `Variables`
37 : *
38 : * The `volume_fields` are all projected into the `face_fields` in the direction
39 : * `direction`. The tags in `VolumeVarsTagsList` must be a contiguous subset of
40 : * the tags in `FaceVarsTagsList`. That is, `FaceVarsTagsList` must be
41 : * equivalent to `tmpl::append<Before, VolumeVarsTagsList, After>` where
42 : * `Before` and `After` are `tmpl::list`s of arbitrary size. This is because the
43 : * projection is applied on all of the tensor components of the volume variables
44 : * and is written into contiguous memory on the boundary.
45 : *
46 : * In general, this function will be used for projecting all the evolved
47 : * variables or all the volume fluxes to the faces. The function
48 : * `dg::project_tensors_to_boundary()` should be used for projecting
49 : * individual tensors to the face.
50 : *
51 : * \note This function works for both Gauss and Gauss-Lobatto uniform meshes.
52 : */
53 : template <typename VolumeVarsTagsList, typename FaceVarsTagsList, size_t Dim>
54 1 : void project_contiguous_data_to_boundary(
55 : const gsl::not_null<Variables<FaceVarsTagsList>*> face_fields,
56 : const Variables<VolumeVarsTagsList>& volume_fields,
57 : const Mesh<Dim>& volume_mesh, const Direction<Dim>& direction) {
58 : static_assert(tmpl::size<VolumeVarsTagsList>::value != 0,
59 : "Must have non-zero number of volume fields");
60 : static_assert(tmpl::size<FaceVarsTagsList>::value >=
61 : tmpl::size<VolumeVarsTagsList>::value,
62 : "There must not be more volume tags than there are face tags.");
63 : static_assert(
64 : tmpl::list_contains_v<FaceVarsTagsList, tmpl::front<VolumeVarsTagsList>>,
65 : "The first tag of VolumeVarsTagsList is not in the face tags. The "
66 : "VolumeVarsTagsList must be a subset of the FaceVarsTagsList");
67 : static_assert(
68 : tmpl::list_contains_v<FaceVarsTagsList, tmpl::back<VolumeVarsTagsList>>,
69 : "The last tag of VolumeVarsTagsList is not in the face tags. The "
70 : "VolumeVarsTagsList must be a subset of the FaceVarsTagsList");
71 : using face_vars_excluding_extras_at_end = tmpl::front<
72 : tmpl::split_at<FaceVarsTagsList,
73 : tmpl::next<tmpl::index_of<
74 : FaceVarsTagsList, tmpl::back<VolumeVarsTagsList>>>>>;
75 : using front_face_vars_split = tmpl::split_at<
76 : face_vars_excluding_extras_at_end,
77 : tmpl::index_of<FaceVarsTagsList, tmpl::front<VolumeVarsTagsList>>>;
78 : using volume_vars_face_subset_list = tmpl::back<front_face_vars_split>;
79 : static_assert(
80 : std::is_same_v<volume_vars_face_subset_list, VolumeVarsTagsList>,
81 : "The VolumeVarsTagsList must be a subset of the FaceVarsTagsList.");
82 : constexpr const size_t number_of_independent_components =
83 : Variables<VolumeVarsTagsList>::number_of_independent_components;
84 : using first_volume_tag = tmpl::front<VolumeVarsTagsList>;
85 : const size_t sliced_dim = direction.dimension();
86 : if (volume_mesh.quadrature(sliced_dim) == Spectral::Quadrature::Gauss) {
87 : const Matrix identity{};
88 : auto interpolation_matrices = make_array<Dim>(std::cref(identity));
89 : const auto& matrix = Spectral::boundary_interpolation_matrices(
90 : volume_mesh.slice_through(sliced_dim));
91 : gsl::at(interpolation_matrices, sliced_dim) =
92 : direction.side() == Side::Upper ? matrix.second : matrix.first;
93 :
94 : auto& first_face_field = get<first_volume_tag>(*face_fields);
95 : auto& first_volume_field = get<first_volume_tag>(volume_fields);
96 :
97 : // The size is the number of tensor components we are projecting times the
98 : // number of grid points on the face. Note that this is _not_ equal to the
99 : // size of face_fields->size() since face_fields is a superset of the
100 : // volume variables.
101 : DataVector face_view{
102 : first_face_field[0].data(),
103 : first_face_field[0].size() * number_of_independent_components};
104 :
105 : apply_matrices(make_not_null(&face_view), interpolation_matrices,
106 : // NOLINTNEXTLINE(cppcoreguidelines-pro-type-const-cast)
107 : DataVector{const_cast<double*>(first_volume_field[0].data()),
108 : first_volume_field[0].size() *
109 : number_of_independent_components},
110 : volume_mesh.extents());
111 : } else {
112 : const size_t fixed_index = direction.side() == Side::Upper
113 : ? volume_mesh.extents(sliced_dim) - 1
114 : : 0;
115 :
116 : const size_t interface_grid_points =
117 : volume_mesh.extents().slice_away(sliced_dim).product();
118 : const size_t volume_grid_points = volume_mesh.number_of_grid_points();
119 :
120 : const double* vars_data = volume_fields.data();
121 : // Since the face fields are a superset of the volume tags we need to find
122 : // the first volume tag on the face and get the pointer for that.
123 : double* interface_vars_data = get<first_volume_tag>(*face_fields)[0].data();
124 :
125 : // The reason we can't use data_on_slice is because the volume and face tags
126 : // are not the same, but data_on_slice assumes they are. In general, this
127 : // function should replace data_on_slice in the long term since in
128 : // additional to supporting different volume and face tags, it also supports
129 : // Gauss and Gauss-Lobatto points.
130 : //
131 : // Run the SliceIterator as the outer-most loop since incrementing the slice
132 : // iterator is surprisingly expensive.
133 : for (SliceIterator si(volume_mesh.extents(), sliced_dim, fixed_index); si;
134 : ++si) {
135 : for (size_t i = 0; i < number_of_independent_components; ++i) {
136 : // clang-tidy: do not use pointer arithmetic
137 : interface_vars_data[si.slice_offset() + // NOLINT
138 : i * interface_grid_points] = // NOLINT
139 : vars_data[si.volume_offset() + i * volume_grid_points]; // NOLINT
140 : }
141 : }
142 : }
143 : }
144 :
145 : /*!
146 : * \brief Projects a subset of the tensors in the `volume_fields` onto the face
147 : *
148 : * The tensors to project are listed in the `TagsToProjectList`.
149 : *
150 : * \note This function works for both Gauss and Gauss-Lobatto uniform meshes.
151 : */
152 : template <typename TagsToProjectList, typename VolumeVarsTagsList,
153 : typename FaceVarsTagsList, size_t Dim>
154 1 : void project_tensors_to_boundary(
155 : const gsl::not_null<Variables<FaceVarsTagsList>*> face_fields,
156 : const Variables<VolumeVarsTagsList>& volume_fields,
157 : const Mesh<Dim>& volume_mesh, const Direction<Dim>& direction) {
158 : static_assert(tmpl::size<VolumeVarsTagsList>::value != 0,
159 : "Must have non-zero number of volume fields");
160 : static_assert(
161 : tmpl::size<
162 : tmpl::list_difference<TagsToProjectList, FaceVarsTagsList>>::value ==
163 : 0,
164 : "All of the tags in TagsToProjectList must be in FaceVarsTagsList");
165 : static_assert(
166 : tmpl::size<tmpl::list_difference<TagsToProjectList,
167 : VolumeVarsTagsList>>::value == 0,
168 : "All of the tags in TagsToProjectList must be in VolumeVarsTagsList");
169 : const size_t sliced_dim = direction.dimension();
170 : if (volume_mesh.quadrature(sliced_dim) == Spectral::Quadrature::Gauss) {
171 : const Matrix identity{};
172 : auto interpolation_matrices = make_array<Dim>(std::cref(identity));
173 : const std::pair<Matrix, Matrix>& matrices =
174 : Spectral::boundary_interpolation_matrices(
175 : volume_mesh.slice_through(sliced_dim));
176 : gsl::at(interpolation_matrices, sliced_dim) =
177 : direction.side() == Side::Upper ? matrices.second : matrices.first;
178 : tmpl::for_each<TagsToProjectList>([&face_fields, &interpolation_matrices,
179 : &volume_fields,
180 : &volume_mesh](auto tag_v) {
181 : using tag = typename decltype(tag_v)::type;
182 : auto& face_field = get<tag>(*face_fields);
183 : const auto& volume_field = get<tag>(volume_fields);
184 : DataVector face_view{face_field[0].data(),
185 : face_field[0].size() * face_field.size()};
186 : apply_matrices(make_not_null(&face_view), interpolation_matrices,
187 : // NOLINTNEXTLINE(cppcoreguidelines-pro-type-const-cast)
188 : DataVector{const_cast<double*>(volume_field[0].data()),
189 : volume_field[0].size() * volume_field.size()},
190 : volume_mesh.extents());
191 : });
192 : } else {
193 : const size_t fixed_index = direction.side() == Side::Upper
194 : ? volume_mesh.extents(sliced_dim) - 1
195 : : 0;
196 :
197 : const size_t interface_grid_points =
198 : volume_mesh.extents().slice_away(sliced_dim).product();
199 : const size_t volume_grid_points = volume_mesh.number_of_grid_points();
200 :
201 : // Run the SliceIterator as the outer-most loop since incrementing the slice
202 : // iterator is surprisingly expensive.
203 : for (SliceIterator si(volume_mesh.extents(), sliced_dim, fixed_index); si;
204 : ++si) {
205 : tmpl::for_each<TagsToProjectList>([&face_fields, interface_grid_points,
206 : &si, &volume_fields,
207 : volume_grid_points](auto tag_v) {
208 : using tag = typename decltype(tag_v)::type;
209 :
210 : const double* vars_data = get<tag>(volume_fields)[0].data();
211 : double* interface_vars_data = get<tag>(*face_fields)[0].data();
212 : static constexpr size_t number_of_independent_components_in_tensor =
213 : std::decay_t<decltype(get<tag>(volume_fields))>::size();
214 :
215 : for (size_t i = 0; i < number_of_independent_components_in_tensor;
216 : ++i) {
217 : // clang-tidy: do not use pointer arithmetic
218 : interface_vars_data[si.slice_offset() + // NOLINT
219 : i * interface_grid_points] = // NOLINT
220 : vars_data[si.volume_offset() + i * volume_grid_points]; // NOLINT
221 : }
222 : });
223 : }
224 : }
225 : }
226 :
227 : /// @{
228 : /*!
229 : * \brief Projects a tensor to the face
230 : *
231 : * \note This function works for both Gauss and Gauss-Lobatto uniform meshes.
232 : */
233 : template <typename Symm, typename IndexList, size_t Dim>
234 1 : void project_tensor_to_boundary(
235 : const gsl::not_null<Tensor<DataVector, Symm, IndexList>*> face_field,
236 : const Tensor<DataVector, Symm, IndexList>& volume_field,
237 : const Mesh<Dim>& volume_mesh, const Direction<Dim>& direction) {
238 : const size_t sliced_dim = direction.dimension();
239 : if (volume_mesh.quadrature(sliced_dim) == Spectral::Quadrature::Gauss) {
240 : const Matrix identity{};
241 : auto interpolation_matrices = make_array<Dim>(std::cref(identity));
242 : const std::pair<Matrix, Matrix>& matrices =
243 : Spectral::boundary_interpolation_matrices(
244 : volume_mesh.slice_through(sliced_dim));
245 : gsl::at(interpolation_matrices, sliced_dim) =
246 : direction.side() == Side::Upper ? matrices.second : matrices.first;
247 : for (size_t tensor_storage_index = 0;
248 : tensor_storage_index < volume_field.size(); ++tensor_storage_index) {
249 : apply_matrices(make_not_null(&(*face_field)[tensor_storage_index]),
250 : interpolation_matrices, volume_field[tensor_storage_index],
251 : volume_mesh.extents());
252 : }
253 : } else {
254 : const size_t fixed_index = direction.side() == Side::Upper
255 : ? volume_mesh.extents(sliced_dim) - 1
256 : : 0;
257 :
258 : // Run the SliceIterator as the outer-most loop since incrementing the slice
259 : // iterator is surprisingly expensive.
260 : for (SliceIterator si(volume_mesh.extents(), sliced_dim, fixed_index); si;
261 : ++si) {
262 : for (size_t tensor_storage_index = 0;
263 : tensor_storage_index < volume_field.size(); ++tensor_storage_index) {
264 : (*face_field)[tensor_storage_index][si.slice_offset()] =
265 : volume_field[tensor_storage_index][si.volume_offset()];
266 : }
267 : }
268 : }
269 : }
270 :
271 : template <typename Symm, typename IndexList, size_t Dim>
272 1 : Tensor<DataVector, Symm, IndexList> project_tensor_to_boundary(
273 : const Tensor<DataVector, Symm, IndexList>& volume_field,
274 : const Mesh<Dim>& volume_mesh, const Direction<Dim>& direction) {
275 : Tensor<DataVector, Symm, IndexList> face_field{
276 : volume_mesh.slice_away(direction.dimension()).number_of_grid_points()};
277 : project_tensor_to_boundary(make_not_null(&face_field), volume_field,
278 : volume_mesh, direction);
279 : return face_field;
280 : }
281 : /// @}
282 :
283 : } // namespace dg
|
