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(
160  tmpl::size<
162  0,
163  "All of the tags in TagsToProjectList must be in FaceVarsTagsList");
164  static_assert(
165  tmpl::size<tmpl::list_difference<TagsToProjectList,
166  VolumeVarsTagsList>>::value == 0,
167  "All of the tags in TagsToProjectList must be in VolumeVarsTagsList");
168  const size_t sliced_dim = direction.dimension();
169  if (volume_mesh.quadrature(sliced_dim) == Spectral::Quadrature::Gauss) {
170  const Matrix identity{};
171  auto interpolation_matrices = make_array<Dim>(std::cref(identity));
172  const std::pair<Matrix, Matrix>& matrices =
174  volume_mesh.slice_through(sliced_dim));
175  gsl::at(interpolation_matrices, sliced_dim) =
176  direction.side() == Side::Upper ? matrices.second : matrices.first;
177  tmpl::for_each<TagsToProjectList>([&face_fields, &interpolation_matrices,
178  &volume_fields,
179  &volume_mesh](auto tag_v) noexcept {
180  using tag = typename decltype(tag_v)::type;
181  auto& face_field = get<tag>(*face_fields);
182  const auto& volume_field = get<tag>(volume_fields);
183  DataVector face_view{face_field[0].data(),
184  face_field[0].size() * face_field.size()};
185  apply_matrices(make_not_null(&face_view), interpolation_matrices,
186  // NOLINTNEXTLINE(cppcoreguidelines-pro-type-const-cast)
187  DataVector{const_cast<double*>(volume_field[0].data()),
188  volume_field[0].size() * volume_field.size()},
189  volume_mesh.extents());
190  });
191  } else {
192  const size_t fixed_index = direction.side() == Side::Upper
193  ? volume_mesh.extents(sliced_dim) - 1
194  : 0;
195 
196  const size_t interface_grid_points =
197  volume_mesh.extents().slice_away(sliced_dim).product();
198  const size_t volume_grid_points = volume_mesh.number_of_grid_points();
199 
200  // Run the SliceIterator as the outer-most loop since incrementing the slice
201  // iterator is surprisingly expensive.
202  for (SliceIterator si(volume_mesh.extents(), sliced_dim, fixed_index); si;
203  ++si) {
204  tmpl::for_each<TagsToProjectList>([&face_fields, interface_grid_points,
205  &si, &volume_fields,
206  volume_grid_points](
207  auto tag_v) noexcept {
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  * \brief Projects a tensor to the face
229  *
230  * \note This function works for both Gauss and Gauss-Lobatto uniform meshes.
231  */
232 template <typename Symm, typename IndexList, size_t Dim>
234  const gsl::not_null<Tensor<DataVector, Symm, IndexList>*> face_field,
235  const Tensor<DataVector, Symm, IndexList>& volume_field,
236  const Mesh<Dim>& volume_mesh, const Direction<Dim>& direction) noexcept {
237  const size_t sliced_dim = direction.dimension();
238  if (volume_mesh.quadrature(sliced_dim) ==
239  Spectral::Quadrature::Gauss) {
240  const Matrix identity{};
241  auto interpolation_matrices = make_array<Dim>(std::cref(identity));
242  const std::pair<Matrix, Matrix>& 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 } // 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:233
Spectral.hpp
Direction< Dim >
cstddef
Assert.hpp
SliceIterator
Iterate over a (dim-1)-dimensional slice.
Definition: SliceIterator.hpp:23
identity
tnsr::Ij< DataType, Dim, Frame::NoFrame > identity(const DataType &used_for_type) noexcept
returns the Identity matrix
Definition: Identity.hpp:14
Spectral::boundary_interpolation_matrices
const std::pair< Matrix, Matrix > & boundary_interpolation_matrices(const Mesh< 1 > &mesh) noexcept
Matrices that interpolate to the lower and upper boundaries of the element.
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:49
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
brigand::list_difference
fold< Sequence2, Sequence1, lazy::remove< _state, _element > > list_difference
Obtain the elements of Sequence1 that are not in Sequence2.
Definition: TMPL.hpp:589
type_traits
TMPL.hpp
Mesh.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13