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"
15 #include "ErrorHandling/Assert.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 
85  const Mesh<Dim> uniform_gauss_mesh(volume_mesh.extents(0),
86  volume_mesh.basis(0),
87  Spectral::Quadrature::Gauss);
88  if (volume_mesh == uniform_gauss_mesh) {
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(direction.dimension()));
93  gsl::at(interpolation_matrices, direction.dimension()) =
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  DataVector face_view{
104  first_face_field[0].data(),
105  first_face_field[0].size() * number_of_independent_components};
106 
107  apply_matrices(make_not_null(&face_view), interpolation_matrices,
108  // NOLINTNEXTLINE(cppcoreguidelines-pro-type-const-cast)
109  DataVector{const_cast<double*>(first_volume_field[0].data()),
110  first_volume_field[0].size() *
111  number_of_independent_components},
112  volume_mesh.extents());
113  } else {
114  ASSERT(Mesh<Dim>(volume_mesh.extents(0), volume_mesh.basis(0),
115  Spectral::Quadrature::GaussLobatto) == volume_mesh,
116  "The current implementation assumes the mesh be either a uniform "
117  "Gauss or Gauss-Lobatto mesh, but got "
118  << volume_mesh << ".");
119  const size_t sliced_dim = direction.dimension();
120  const size_t fixed_index = direction.side() == Side::Upper
121  ? volume_mesh.extents(sliced_dim) - 1
122  : 0;
123 
124  const size_t interface_grid_points =
125  volume_mesh.extents().slice_away(sliced_dim).product();
126  const size_t volume_grid_points = volume_mesh.number_of_grid_points();
127 
128  const double* vars_data = volume_fields.data();
129  // Since the face fields are a superset of the volume tags we need to find
130  // the first volume tag on the face and get the pointer for that.
131  double* interface_vars_data = get<first_volume_tag>(*face_fields)[0].data();
132 
133  // The reason we can't use data_on_slice is because the volume and face tags
134  // are not the same, but data_on_slice assumes they are. In general, this
135  // function should replace data_on_slice in the long term since in
136  // additional to supporting different volume and face tags, it also supports
137  // Gauss and Gauss-Lobatto points.
138  //
139  // Run the SliceIterator as the outer-most loop since incrementing the slice
140  // iterator is surprisingly expensive.
141  for (SliceIterator si(volume_mesh.extents(), sliced_dim, fixed_index); si;
142  ++si) {
143  for (size_t i = 0; i < number_of_independent_components; ++i) {
144  // clang-tidy: do not use pointer arithmetic
145  interface_vars_data[si.slice_offset() + // NOLINT
146  i * interface_grid_points] = // NOLINT
147  vars_data[si.volume_offset() + i * volume_grid_points]; // NOLINT
148  }
149  }
150  }
151 }
152 
153 /*!
154  * \brief Projects a subset of the tensors in the `volume_fields` onto the face
155  *
156  * The tensors to project are listed in the `TagsToProjectList`.
157  *
158  * \note This function works for both Gauss and Gauss-Lobatto uniform meshes.
159  */
160 template <typename TagsToProjectList, typename VolumeVarsTagsList,
161  typename FaceVarsTagsList, size_t Dim>
163  const gsl::not_null<Variables<FaceVarsTagsList>*> face_fields,
164  const Variables<VolumeVarsTagsList>& volume_fields,
165  const Mesh<Dim>& volume_mesh, const Direction<Dim>& direction) noexcept {
166  static_assert(tmpl::size<VolumeVarsTagsList>::value != 0,
167  "Must have non-zero number of volume fields");
168  static_assert(tmpl::size<FaceVarsTagsList>::value >=
169  tmpl::size<VolumeVarsTagsList>::value,
170  "There must not be more volume tags than there are face tags.");
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 Mesh<Dim> uniform_gauss_mesh(volume_mesh.extents(0),
181  volume_mesh.basis(0),
182  Spectral::Quadrature::Gauss);
183  if (volume_mesh.quadrature() == uniform_gauss_mesh.quadrature()) {
184  const Matrix identity{};
185  auto interpolation_matrices = make_array<Dim>(std::cref(identity));
186  const std::pair<Matrix, Matrix>& matrices =
187  Spectral::boundary_interpolation_matrices(
188  volume_mesh.slice_through(direction.dimension()));
189  gsl::at(interpolation_matrices, direction.dimension()) =
190  direction.side() == Side::Upper ? matrices.second : matrices.first;
191  tmpl::for_each<TagsToProjectList>([&face_fields, &interpolation_matrices,
192  &volume_fields,
193  &volume_mesh](auto tag_v) noexcept {
194  using tag = typename decltype(tag_v)::type;
195  auto& face_field = get<tag>(*face_fields);
196  const auto& volume_field = get<tag>(volume_fields);
197  DataVector face_view{face_field[0].data(),
198  face_field[0].size() * face_field.size()};
199  apply_matrices(make_not_null(&face_view), interpolation_matrices,
200  // NOLINTNEXTLINE(cppcoreguidelines-pro-type-const-cast)
201  DataVector{const_cast<double*>(volume_field[0].data()),
202  volume_field[0].size() * volume_field.size()},
203  volume_mesh.extents());
204  });
205  } else {
206  ASSERT(Mesh<Dim>(volume_mesh.extents(0), volume_mesh.basis(0),
207  Spectral::Quadrature::GaussLobatto) == volume_mesh,
208  "The current implementation assumes the mesh be either a uniform "
209  "Gauss or Gauss-Lobatto mesh, but got "
210  << volume_mesh);
211 
212  const size_t sliced_dim = direction.dimension();
213  const size_t fixed_index = direction.side() == Side::Upper
214  ? volume_mesh.extents(sliced_dim) - 1
215  : 0;
216 
217  const size_t interface_grid_points =
218  volume_mesh.extents().slice_away(sliced_dim).product();
219  const size_t volume_grid_points = volume_mesh.number_of_grid_points();
220 
221  // Run the SliceIterator as the outer-most loop since incrementing the slice
222  // iterator is surprisingly expensive.
223  for (SliceIterator si(volume_mesh.extents(), sliced_dim, fixed_index); si;
224  ++si) {
225  tmpl::for_each<TagsToProjectList>([&face_fields, interface_grid_points,
226  &si, &volume_fields,
227  volume_grid_points](
228  auto tag_v) noexcept {
229  using tag = typename decltype(tag_v)::type;
230 
231  const double* vars_data = get<tag>(volume_fields)[0].data();
232  double* interface_vars_data = get<tag>(*face_fields)[0].data();
233  static constexpr size_t number_of_independent_components_in_tensor =
234  std::decay_t<decltype(get<tag>(volume_fields))>::size();
235 
236  for (size_t i = 0; i < number_of_independent_components_in_tensor;
237  ++i) {
238  // clang-tidy: do not use pointer arithmetic
239  interface_vars_data[si.slice_offset() + // NOLINT
240  i * interface_grid_points] = // NOLINT
241  vars_data[si.volume_offset() + i * volume_grid_points]; // NOLINT
242  }
243  });
244  }
245  }
246 }
247 } // 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
Spectral.hpp
Direction
Definition: Direction.hpp:23
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:42
std::decay_t
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:51
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:162
Variables.hpp
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:47
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
Mesh::quadrature
const std::array< Spectral::Quadrature, Dim > & quadrature() const noexcept
The quadrature chosen in each dimension of the grid.
Definition: Mesh.hpp:158
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: Gsl.hpp:183