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