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/BoundaryInterpolationMatrices.hpp"
16 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
17 : #include "NumericalAlgorithms/Spectral/Quadrature.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 : using VectorType = typename Variables<FaceVarsTagsList>::vector_type;
86 : using ValueType = typename VectorType::value_type;
87 : const size_t sliced_dim = direction.dimension();
88 : if (volume_mesh.quadrature(sliced_dim) == Spectral::Quadrature::Gauss) {
89 : const Matrix identity{};
90 : auto interpolation_matrices = make_array<Dim>(std::cref(identity));
91 : const auto& matrix = Spectral::boundary_interpolation_matrices(
92 : volume_mesh.slice_through(sliced_dim));
93 : gsl::at(interpolation_matrices, sliced_dim) =
94 : direction.side() == Side::Upper ? matrix.second : matrix.first;
95 :
96 : auto& first_face_field = get<first_volume_tag>(*face_fields);
97 : auto& first_volume_field = get<first_volume_tag>(volume_fields);
98 :
99 : // The size is the number of tensor components we are projecting times the
100 : // number of grid points on the face. Note that this is _not_ equal to the
101 : // size of face_fields->size() since face_fields is a superset of the
102 : // volume variables.
103 : VectorType face_view{
104 : first_face_field[0].data(),
105 : first_face_field[0].size() * number_of_independent_components};
106 :
107 : apply_matrices(
108 : make_not_null(&face_view), interpolation_matrices,
109 : VectorType{
110 : // NOLINTNEXTLINE(cppcoreguidelines-pro-type-const-cast)
111 : const_cast<ValueType*>(first_volume_field[0].data()),
112 : first_volume_field[0].size() * number_of_independent_components},
113 : volume_mesh.extents());
114 : } else {
115 : ASSERT(volume_mesh.quadrature(sliced_dim) ==
116 : Spectral::Quadrature::GaussLobatto or
117 : (volume_mesh.quadrature(sliced_dim) ==
118 : Spectral::Quadrature::GaussRadauUpper and
119 : direction.side() == Side::Upper),
120 : "Got quadrature without boundary collocation point at "
121 : << direction << " with " << volume_mesh.quadrature(sliced_dim));
122 : const size_t fixed_index = direction.side() == Side::Upper
123 : ? volume_mesh.extents(sliced_dim) - 1
124 : : 0;
125 :
126 : const size_t interface_grid_points =
127 : volume_mesh.extents().slice_away(sliced_dim).product();
128 : const size_t volume_grid_points = volume_mesh.number_of_grid_points();
129 :
130 : const ValueType* vars_data = volume_fields.data();
131 : // Since the face fields are a superset of the volume tags we need to find
132 : // the first volume tag on the face and get the pointer for that.
133 : ValueType* interface_vars_data =
134 : get<first_volume_tag>(*face_fields)[0].data();
135 :
136 : // The reason we can't use data_on_slice is because the volume and face tags
137 : // are not the same, but data_on_slice assumes they are. In general, this
138 : // function should replace data_on_slice in the long term since in
139 : // additional to supporting different volume and face tags, it also supports
140 : // Gauss and Gauss-Lobatto points.
141 : //
142 : // Run the SliceIterator as the outer-most loop since incrementing the slice
143 : // iterator is surprisingly expensive.
144 : for (SliceIterator si(volume_mesh.extents(), sliced_dim, fixed_index); si;
145 : ++si) {
146 : for (size_t i = 0; i < number_of_independent_components; ++i) {
147 : // clang-tidy: do not use pointer arithmetic
148 : interface_vars_data[si.slice_offset() + // NOLINT
149 : i * interface_grid_points] = // NOLINT
150 : vars_data[si.volume_offset() + i * volume_grid_points]; // NOLINT
151 : }
152 : }
153 : }
154 : }
155 :
156 : /*!
157 : * \brief Projects a subset of the tensors in the `volume_fields` onto the face
158 : *
159 : * The tensors to project are listed in the `TagsToProjectList`.
160 : *
161 : * \note This function works for both Gauss and Gauss-Lobatto uniform meshes.
162 : */
163 : template <typename TagsToProjectList, typename VolumeVarsTagsList,
164 : typename FaceVarsTagsList, size_t Dim>
165 1 : void project_tensors_to_boundary(
166 : const gsl::not_null<Variables<FaceVarsTagsList>*> face_fields,
167 : const Variables<VolumeVarsTagsList>& volume_fields,
168 : const Mesh<Dim>& volume_mesh, const Direction<Dim>& direction) {
169 : static_assert(tmpl::size<VolumeVarsTagsList>::value != 0,
170 : "Must have non-zero number of volume fields");
171 : static_assert(
172 : tmpl::size<
173 : tmpl::list_difference<TagsToProjectList, FaceVarsTagsList>>::value ==
174 : 0,
175 : "All of the tags in TagsToProjectList must be in FaceVarsTagsList");
176 : static_assert(
177 : tmpl::size<tmpl::list_difference<TagsToProjectList,
178 : VolumeVarsTagsList>>::value == 0,
179 : "All of the tags in TagsToProjectList must be in VolumeVarsTagsList");
180 : const size_t sliced_dim = direction.dimension();
181 : if (volume_mesh.quadrature(sliced_dim) == Spectral::Quadrature::Gauss) {
182 : const Matrix identity{};
183 : auto interpolation_matrices = make_array<Dim>(std::cref(identity));
184 : const std::pair<Matrix, Matrix>& matrices =
185 : Spectral::boundary_interpolation_matrices(
186 : volume_mesh.slice_through(sliced_dim));
187 : gsl::at(interpolation_matrices, sliced_dim) =
188 : direction.side() == Side::Upper ? matrices.second : matrices.first;
189 : tmpl::for_each<TagsToProjectList>([&face_fields, &interpolation_matrices,
190 : &volume_fields,
191 : &volume_mesh](auto tag_v) {
192 : using tag = typename decltype(tag_v)::type;
193 : auto& face_field = get<tag>(*face_fields);
194 : const auto& volume_field = get<tag>(volume_fields);
195 : DataVector face_view{face_field[0].data(),
196 : face_field[0].size() * face_field.size()};
197 : apply_matrices(make_not_null(&face_view), interpolation_matrices,
198 : // NOLINTNEXTLINE(cppcoreguidelines-pro-type-const-cast)
199 : DataVector{const_cast<double*>(volume_field[0].data()),
200 : volume_field[0].size() * volume_field.size()},
201 : volume_mesh.extents());
202 : });
203 : } else {
204 : ASSERT(volume_mesh.quadrature(sliced_dim) ==
205 : Spectral::Quadrature::GaussLobatto or
206 : (volume_mesh.quadrature(sliced_dim) ==
207 : Spectral::Quadrature::GaussRadauUpper and
208 : direction.side() == Side::Upper),
209 : "Got quadrature without boundary collocation point at "
210 : << direction << " with " << volume_mesh.quadrature(sliced_dim));
211 : const size_t fixed_index = direction.side() == Side::Upper
212 : ? volume_mesh.extents(sliced_dim) - 1
213 : : 0;
214 :
215 : const size_t interface_grid_points =
216 : volume_mesh.extents().slice_away(sliced_dim).product();
217 : const size_t volume_grid_points = volume_mesh.number_of_grid_points();
218 :
219 : // Run the SliceIterator as the outer-most loop since incrementing the slice
220 : // iterator is surprisingly expensive.
221 : for (SliceIterator si(volume_mesh.extents(), sliced_dim, fixed_index); si;
222 : ++si) {
223 : tmpl::for_each<TagsToProjectList>([&face_fields, interface_grid_points,
224 : &si, &volume_fields,
225 : volume_grid_points](auto tag_v) {
226 : using tag = typename decltype(tag_v)::type;
227 :
228 : const double* vars_data = get<tag>(volume_fields)[0].data();
229 : double* interface_vars_data = get<tag>(*face_fields)[0].data();
230 : static constexpr size_t number_of_independent_components_in_tensor =
231 : std::decay_t<decltype(get<tag>(volume_fields))>::size();
232 :
233 : for (size_t i = 0; i < number_of_independent_components_in_tensor;
234 : ++i) {
235 : // clang-tidy: do not use pointer arithmetic
236 : interface_vars_data[si.slice_offset() + // NOLINT
237 : i * interface_grid_points] = // NOLINT
238 : vars_data[si.volume_offset() + i * volume_grid_points]; // NOLINT
239 : }
240 : });
241 : }
242 : }
243 : }
244 :
245 : /// @{
246 : /*!
247 : * \brief Projects a tensor to the face
248 : *
249 : * \note This function works for both Gauss and Gauss-Lobatto uniform meshes.
250 : */
251 : template <typename Symm, typename IndexList, size_t Dim>
252 1 : void project_tensor_to_boundary(
253 : const gsl::not_null<Tensor<DataVector, Symm, IndexList>*> face_field,
254 : const Tensor<DataVector, Symm, IndexList>& volume_field,
255 : const Mesh<Dim>& volume_mesh, const Direction<Dim>& direction) {
256 : const size_t sliced_dim = direction.dimension();
257 : if (volume_mesh.quadrature(sliced_dim) == Spectral::Quadrature::Gauss) {
258 : const Matrix identity{};
259 : auto interpolation_matrices = make_array<Dim>(std::cref(identity));
260 : const std::pair<Matrix, Matrix>& matrices =
261 : Spectral::boundary_interpolation_matrices(
262 : volume_mesh.slice_through(sliced_dim));
263 : gsl::at(interpolation_matrices, sliced_dim) =
264 : direction.side() == Side::Upper ? matrices.second : matrices.first;
265 : for (size_t tensor_storage_index = 0;
266 : tensor_storage_index < volume_field.size(); ++tensor_storage_index) {
267 : apply_matrices(make_not_null(&(*face_field)[tensor_storage_index]),
268 : interpolation_matrices, volume_field[tensor_storage_index],
269 : volume_mesh.extents());
270 : }
271 : } else {
272 : ASSERT(volume_mesh.quadrature(sliced_dim) ==
273 : Spectral::Quadrature::GaussLobatto or
274 : (volume_mesh.quadrature(sliced_dim) ==
275 : Spectral::Quadrature::GaussRadauUpper and
276 : direction.side() == Side::Upper),
277 : "Got quadrature without boundary collocation point at "
278 : << direction << " with " << volume_mesh.quadrature(sliced_dim));
279 : const size_t fixed_index = direction.side() == Side::Upper
280 : ? volume_mesh.extents(sliced_dim) - 1
281 : : 0;
282 :
283 : // Run the SliceIterator as the outer-most loop since incrementing the slice
284 : // iterator is surprisingly expensive.
285 : for (SliceIterator si(volume_mesh.extents(), sliced_dim, fixed_index); si;
286 : ++si) {
287 : for (size_t tensor_storage_index = 0;
288 : tensor_storage_index < volume_field.size(); ++tensor_storage_index) {
289 : (*face_field)[tensor_storage_index][si.slice_offset()] =
290 : volume_field[tensor_storage_index][si.volume_offset()];
291 : }
292 : }
293 : }
294 : }
295 :
296 : template <typename Symm, typename IndexList, size_t Dim>
297 1 : Tensor<DataVector, Symm, IndexList> project_tensor_to_boundary(
298 : const Tensor<DataVector, Symm, IndexList>& volume_field,
299 : const Mesh<Dim>& volume_mesh, const Direction<Dim>& direction) {
300 : Tensor<DataVector, Symm, IndexList> face_field{
301 : volume_mesh.slice_away(direction.dimension()).number_of_grid_points()};
302 : project_tensor_to_boundary(make_not_null(&face_field), volume_field,
303 : volume_mesh, direction);
304 : return face_field;
305 : }
306 : /// @}
307 :
308 : } // namespace dg
|
