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