Exact.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
7 
8 #include <array>
9 #include <cstddef>
10 #include <vector>
11 
12 #include "DataStructures/DataVector.hpp"
13 #include "DataStructures/SliceIterator.hpp"
17 #include "Domain/Structure/DirectionMap.hpp"
18 #include "Evolution/DgSubcell/SliceData.hpp"
19 #include "NumericalAlgorithms/FiniteDifference/Reconstruct.hpp"
22 #include "Utilities/Gsl.hpp"
23 #include "Utilities/MakeArray.hpp"
24 
25 namespace TestHelpers::fd::reconstruction {
26 /*!
27  * \brief Test that the reconstruction procedure can exactly reconstruct
28  * polynomials of degree equal to the degree of the basis.
29  *
30  * For example, for a minmod, MC, or other second-order TVD limiter, we must be
31  * able to reconstruct linear functions exactly. For a third-order WENO scheme,
32  * we must be able to reconstruct quadratic functions exactly. This is done
33  * by setting up the function
34  *
35  * \f{align}{
36  * u=\sum_{i=1}^p x^p + y^p + z^p
37  * \f}
38  *
39  * where \f$(x,y,z)=(\xi,\eta+4,\zeta+8)\f$ and \f$\xi,\eta,\zeta\f$ are the
40  * logical coordinates. The translation is done to catch subtle bugs that may
41  * arise from indexing in the wrong dimension.
42  */
43 template <size_t Dim, typename F, typename F1>
44 void test_reconstruction_is_exact_if_in_basis(
45  const size_t max_degree, const size_t points_per_dimension,
46  const size_t stencil_width, const F& invoke_recons,
47  const F1& invoke_reconstruct_neighbor) {
48  const size_t number_of_vars = 2; // arbitrary, 2 is "cheap but not trivial"
49 
50  const Mesh<Dim> mesh{points_per_dimension, Spectral::Basis::FiniteDifference,
51  Spectral::Quadrature::CellCentered};
52  auto logical_coords = logical_coordinates(mesh);
53  // Make the logical coordinates different in each direction
54  for (size_t i = 1; i < Dim; ++i) {
55  logical_coords.get(i) += 4.0 * i;
56  }
57 
58  // Compute polynomial on cell centers in FD cluster of points
59  const auto set_polynomial = [max_degree](
60  const gsl::not_null<DataVector*> var1_ptr,
61  const gsl::not_null<DataVector*> var2_ptr,
62  const auto& local_logical_coords) {
63  *var1_ptr = 0.0;
64  *var2_ptr = 100.0; // some constant offset to distinguish the var values
65  for (size_t degree = 1; degree <= max_degree; ++degree) {
66  for (size_t i = 0; i < Dim; ++i) {
67  *var1_ptr += pow(local_logical_coords.get(i), degree);
68  if (number_of_vars == 2) {
69  *var2_ptr += pow(local_logical_coords.get(i), degree);
70  }
71  }
72  }
73  };
74  std::vector<double> volume_vars(mesh.number_of_grid_points() * number_of_vars,
75  0.0);
76  DataVector var1(volume_vars.data(), mesh.number_of_grid_points());
77  DataVector var2(volume_vars.data() + mesh.number_of_grid_points(), // NOLINT
78  mesh.number_of_grid_points());
79  set_polynomial(&var1, &var2, logical_coords);
80 
81  // Compute the polynomial at the cell center for the neighbor data that we
82  // "received".
83  //
84  // We do this by computing the solution in our entire neighbor, then using
85  // slice_data to get the subset of points that are needed.
87  for (const auto& direction : Direction<Dim>::all_directions()) {
88  auto neighbor_logical_coords = logical_coords;
89  neighbor_logical_coords.get(direction.dimension()) +=
90  direction.sign() * 2.0;
91  std::vector<double> neighbor_vars(
92  mesh.number_of_grid_points() * number_of_vars, 0.0);
93  DataVector neighbor_var1(neighbor_vars.data(),
94  mesh.number_of_grid_points());
95  DataVector neighbor_var2(
96  neighbor_vars.data() + mesh.number_of_grid_points(), // NOLINT
97  mesh.number_of_grid_points());
98  set_polynomial(&neighbor_var1, &neighbor_var2, neighbor_logical_coords);
99 
100  DirectionMap<Dim, bool> directions_to_slice{};
101  directions_to_slice[direction.opposite()] = true;
102  const auto sliced_data = evolution::dg::subcell::detail::slice_data_impl(
103  gsl::make_span(neighbor_vars.data(), neighbor_vars.size()),
104  mesh.extents(), (stencil_width - 1) / 2 + 1, directions_to_slice);
105  REQUIRE(sliced_data.size() == 1);
106  REQUIRE(sliced_data.contains(direction.opposite()));
107  neighbor_data[direction] = sliced_data.at(direction.opposite());
108  }
109 
110  // Note: reconstructed_num_pts assumes isotropic extents
111  const size_t reconstructed_num_pts =
112  (mesh.extents(0) + 1) * mesh.slice_away(0).number_of_grid_points();
113  std::array<std::vector<double>, Dim> reconstructed_upper_side_of_face_vars =
114  make_array<Dim>(
115  std::vector<double>(reconstructed_num_pts * number_of_vars));
116  std::array<std::vector<double>, Dim> reconstructed_lower_side_of_face_vars =
117  make_array<Dim>(
118  std::vector<double>(reconstructed_num_pts * number_of_vars));
119 
120  std::array<gsl::span<double>, Dim> recons_upper_side_of_face{};
121  std::array<gsl::span<double>, Dim> recons_lower_side_of_face{};
122  for (size_t i = 0; i < Dim; ++i) {
123  auto& upper = gsl::at(reconstructed_upper_side_of_face_vars, i);
124  auto& lower = gsl::at(reconstructed_lower_side_of_face_vars, i);
125  gsl::at(recons_upper_side_of_face, i) =
126  gsl::make_span(upper.data(), upper.size());
127  gsl::at(recons_lower_side_of_face, i) =
128  gsl::make_span(lower.data(), lower.size());
129  }
130 
132  for (const auto& [direction, data] : neighbor_data) {
133  ghost_cell_vars[direction] = gsl::make_span(data.data(), data.size());
134  }
135 
136  invoke_recons(make_not_null(&recons_upper_side_of_face),
137  make_not_null(&recons_lower_side_of_face),
138  gsl::make_span(volume_vars.data(), volume_vars.size()),
139  ghost_cell_vars, mesh.extents(), number_of_vars);
140 
141  for (size_t dim = 0; dim < Dim; ++dim) {
142  CAPTURE(dim);
143  // Since the reconstruction is to the same physical points the values should
144  // be the same if the polynomial is representable in the basis of the
145  // reconstructor.
146 
147  CHECK_ITERABLE_APPROX(gsl::at(reconstructed_upper_side_of_face_vars, dim),
148  gsl::at(reconstructed_lower_side_of_face_vars, dim));
149 
150  // Compare to analytic solution on the faces.
151  const auto basis = make_array<Dim>(Spectral::Basis::FiniteDifference);
152  auto quadrature = make_array<Dim>(Spectral::Quadrature::CellCentered);
153  auto extents = make_array<Dim>(points_per_dimension);
154  gsl::at(extents, dim) = points_per_dimension + 1;
155  gsl::at(quadrature, dim) = Spectral::Quadrature::FaceCentered;
156  const Mesh<Dim> face_centered_mesh{extents, basis, quadrature};
157  auto logical_coords_face_centered = logical_coordinates(face_centered_mesh);
158  for (size_t i = 1; i < Dim; ++i) {
159  logical_coords_face_centered.get(i) =
160  logical_coords_face_centered.get(i) + 4.0 * i;
161  }
162 
163  std::vector<double> expected_volume_vars(
164  face_centered_mesh.number_of_grid_points() * number_of_vars, 0.0);
165  DataVector expected_var1(expected_volume_vars.data(),
166  face_centered_mesh.number_of_grid_points());
167  DataVector expected_var2(
168  expected_volume_vars.data() +
169  face_centered_mesh.number_of_grid_points(), // NOLINT
170  face_centered_mesh.number_of_grid_points());
171  set_polynomial(&expected_var1, &expected_var2,
172  logical_coords_face_centered);
173  CHECK_ITERABLE_APPROX(gsl::at(reconstructed_upper_side_of_face_vars, dim),
174  expected_volume_vars);
175 
176  // Test fd::reconstruction::reconstruct_neighbor
177  Index<Dim> ghost_data_extents = mesh.extents();
178  ghost_data_extents[dim] = (stencil_width + 1) / 2;
179  Index<Dim> extents_with_faces = mesh.extents();
180  ++extents_with_faces[dim];
181 
182  const Direction<Dim> upper_direction{dim, Side::Upper};
183  DataVector upper_face_var1{mesh.extents().slice_away(dim).product()};
184  DataVector upper_face_var2{mesh.extents().slice_away(dim).product()};
185  auto& upper_neighbor_data = neighbor_data.at(upper_direction);
186  DataVector upper_neighbor_var1{upper_neighbor_data.data(),
187  ghost_data_extents.product()};
188  DataVector upper_neighbor_var2{
189  // NOLINTNEXTLINE
190  upper_neighbor_data.data() + ghost_data_extents.product(),
191  ghost_data_extents.product()};
192  invoke_reconstruct_neighbor(make_not_null(&upper_face_var1), var1,
193  upper_neighbor_var1, mesh.extents(),
194  ghost_data_extents, upper_direction);
195  invoke_reconstruct_neighbor(make_not_null(&upper_face_var2), var2,
196  upper_neighbor_var2, mesh.extents(),
197  ghost_data_extents, upper_direction);
198  for (SliceIterator si(extents_with_faces, dim, extents_with_faces[dim] - 1);
199  si; ++si) {
200  CHECK(approx(upper_face_var1[si.slice_offset()]) ==
201  expected_var1[si.volume_offset()]);
202  CHECK(approx(upper_face_var2[si.slice_offset()]) ==
203  expected_var2[si.volume_offset()]);
204  }
205 
206  const Direction<Dim> lower_direction{dim, Side::Lower};
207  DataVector lower_face_var1{mesh.extents().slice_away(dim).product()};
208  DataVector lower_face_var2{mesh.extents().slice_away(dim).product()};
209  auto& lower_neighbor_data = neighbor_data.at(lower_direction);
210  DataVector lower_neighbor_var1{lower_neighbor_data.data(),
211  ghost_data_extents.product()};
212  DataVector lower_neighbor_var2{
213  // NOLINTNEXTLINE
214  lower_neighbor_data.data() + ghost_data_extents.product(),
215  ghost_data_extents.product()};
216  invoke_reconstruct_neighbor(make_not_null(&lower_face_var1), var1,
217  lower_neighbor_var1, mesh.extents(),
218  ghost_data_extents, lower_direction);
219  invoke_reconstruct_neighbor(make_not_null(&lower_face_var2), var2,
220  lower_neighbor_var2, mesh.extents(),
221  ghost_data_extents, lower_direction);
222  for (SliceIterator si(extents_with_faces, dim, 0); si; ++si) {
223  CHECK(approx(lower_face_var1[si.slice_offset()]) ==
224  expected_var1[si.volume_offset()]);
225  CHECK(approx(lower_face_var2[si.slice_offset()]) ==
226  expected_var2[si.volume_offset()]);
227  }
228  }
229 }
230 } // namespace TestHelpers::fd::reconstruction
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
pow
constexpr decltype(auto) pow(const T &t) noexcept
Compute t^N where N is an integer (positive or negative)
Definition: ConstantExpressions.hpp:160
vector
MakeArray.hpp
TestingFramework.hpp
Index
Definition: Index.hpp:31
evolution::dg::subcell::fd::mesh
Mesh< Dim > mesh(const Mesh< Dim > &dg_mesh) noexcept
Computes the cell-centered finite-difference mesh from the DG mesh, using grid points per dimension,...
CHECK_ITERABLE_APPROX
#define CHECK_ITERABLE_APPROX(a, b)
A wrapper around Catch's CHECK macro that checks approximate equality of entries in iterable containe...
Definition: TestingFramework.hpp:122
Spectral.hpp
Direction
Definition: Direction.hpp:23
cstddef
SliceIterator
Iterate over a (dim-1)-dimensional slice.
Definition: SliceIterator.hpp:24
array
LogicalCoordinates.hpp
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:46
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:48
DirectionMap
Definition: DirectionMap.hpp:15
gsl::make_span
constexpr span< ElementType > make_span(ElementType *ptr, typename span< ElementType >::index_type count)
Definition: Gsl.hpp:801
logical_coordinates
void logical_coordinates(gsl::not_null< tnsr::I< DataVector, VolumeDim, Frame::Logical > * > logical_coords, const Mesh< VolumeDim > &mesh) noexcept
Compute the logical coordinates in an Element.
Gsl.hpp
Tensor.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
std::data
T data(T... args)
Index::product
constexpr size_t product() const noexcept
The product of the indices. If Dim = 0, the product is defined as 1.
Definition: Index.hpp:72
Mesh.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecThirdOrderPiecewisePolynomial.hpp:13