Python.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 <random>
11 #include <string>
12 #include <vector>
13 
14 #include "DataStructures/DataVector.hpp"
15 #include "DataStructures/Index.hpp"
16 #include "DataStructures/SliceIterator.hpp"
18 #include "Domain/Structure/DirectionMap.hpp"
20 #include "Framework/Pypp.hpp"
23 #include "Utilities/Gsl.hpp"
24 #include "Utilities/MakeArray.hpp"
25 
26 namespace TestHelpers::fd::reconstruction {
27 /*!
28  * \brief Compare the output of reconstruction done once in python and once in
29  * C++.
30  *
31  * - `extents` The extents are the volume extents of the subcell grid in a
32  * DG-subcell scheme. In general, this should be the size of the reconstruction
33  * stencil plus one or larger to provide a good test.
34  * - `stencil_width` the width of the reconstruction stencil. For example, for
35  * minmod or MC this would be 3, for 5-point 5th-order WENO this would
36  * be 5, and for three three-point stencil CWENO this would be 5 as well.
37  * -`python_test_file` the file in which the python function to compare the
38  * output result is defined in.
39  * - `python_function` the python function whose output will be compared to
40  * - `recons_function` an invokable (usually a lambda) that has the following
41  * signature:
42  * \code{.cpp}
43  * (const gsl::not_null<std::array<gsl::span<double>, Dim>*>
44  * reconstructed_upper_side_of_face_vars,
45  * const gsl::not_null<std::array<gsl::span<double>, Dim>*>
46  * reconstructed_lower_side_of_face_vars,
47  * const gsl::span<const double>& volume_vars,
48  * const DirectionMap<Dim, gsl::span<const double>>& ghost_cell_vars,
49  * const Index<Dim>& volume_extents, const size_t number_of_variables)
50  * \endcode
51  */
52 template <size_t Dim, typename ReconsFunction, typename F1>
53 void test_with_python(const Index<Dim>& extents, const size_t stencil_width,
54  const std::string& python_test_file,
55  const std::string& python_function,
56  const ReconsFunction& recons_function,
57  const F1& invoke_reconstruct_neighbor) {
58  CAPTURE(extents);
59  CAPTURE(Dim);
60  MAKE_GENERATOR(gen);
61  std::uniform_real_distribution<> dist(-1.0, 1.0);
62 
63  const size_t ghost_zone_size = (stencil_width - 1) / 2;
64  const size_t ghost_cells_from_neighbor = ghost_zone_size + 1;
65  Index<Dim> extents_with_ghosts = extents;
66  for (size_t i = 0; i < Dim; ++i) {
67  // 2 because we extend in both upper and lower direction.
68  extents_with_ghosts[i] += 2 * ghost_cells_from_neighbor;
69  }
70  const size_t number_of_vars = 1;
71  std::vector<double> vars_with_ghosts(extents_with_ghosts.product() *
72  number_of_vars);
73  fill_with_random_values(make_not_null(&vars_with_ghosts), make_not_null(&gen),
74  make_not_null(&dist));
75  // Copy out interior and neighbor data since that's what we will be passing
76  // to the C++ implementation.
77  const auto in_volume = [ghost_cells_from_neighbor, extents](
78  const auto& idx, const size_t dim) {
79  return idx >= ghost_cells_from_neighbor and
80  idx < extents[dim] + ghost_cells_from_neighbor;
81  };
82  const auto in_upper_neighbor = [ghost_cells_from_neighbor, extents](
83  const auto& idx, const size_t dim) {
84  return idx >= extents[dim] + ghost_cells_from_neighbor;
85  };
86  const auto in_lower_neighbor = [ghost_cells_from_neighbor](const auto& idx) {
87  return idx < ghost_cells_from_neighbor;
88  };
89 
91  for (const auto& direction : Direction<Dim>::all_directions()) {
92  ghost_cells[direction] = std::vector<double>{};
93  }
94  std::vector<double> volume_vars(extents.product() * number_of_vars);
95  for (size_t k = 0; k < (Dim > 2 ? extents_with_ghosts[2] : 1); ++k) {
96  for (size_t j = 0; j < (Dim > 1 ? extents_with_ghosts[1] : 1); ++j) {
97  for (size_t i = 0; i < extents_with_ghosts[0]; ++i) {
98  if constexpr (Dim > 1) {
99  // Check if we are in a Voronoi neighbor (ignore then)
100  std::vector<bool> has_neighbor(Dim);
101  has_neighbor.at(0) = in_lower_neighbor(i) or in_upper_neighbor(i, 0);
102  has_neighbor.at(1) = in_lower_neighbor(j) or in_upper_neighbor(j, 1);
103  if constexpr (Dim > 2) {
104  has_neighbor.at(2) =
105  in_lower_neighbor(k) or in_upper_neighbor(k, 2);
106  }
107  if (std::count(has_neighbor.begin(), has_neighbor.end(), true) >= 2) {
108  continue;
109  }
110  }
111 
112  std::array volume_with_ghost_indices{i, j, k};
113  Index<Dim> volume_with_ghosts_index{};
114  for (size_t d = 0; d < Dim; ++d) {
115  volume_with_ghosts_index[d] = gsl::at(volume_with_ghost_indices, d);
116  }
117 
118  // If in volume, copy over...
119  if (in_volume(i, 0) and
120  (Dim == 1 or (in_volume(j, 1) and (Dim == 2 or in_volume(k, 2))))) {
121  Index<Dim> volume_index{};
122  for (size_t d = 0; d < Dim; ++d) {
123  volume_index[d] = gsl::at(volume_with_ghost_indices, d) -
124  (extents_with_ghosts[d] - extents[d]) / 2;
125  }
126 
127  volume_vars[collapsed_index(volume_index, extents)] =
128  vars_with_ghosts[collapsed_index(volume_with_ghosts_index,
129  extents_with_ghosts)];
130  } else {
131  // We are dealing with something in the neighbor data. We need to
132  // identify which neighbor and then set.
133  Side side = Side::Lower;
134  size_t dimension = 0;
135  if (in_upper_neighbor(i, 0)) {
136  side = Side::Upper;
137  }
138  if (Dim > 1 and in_lower_neighbor(j)) {
139  dimension = 1;
140  }
141  if (Dim > 1 and in_upper_neighbor(j, 1)) {
142  side = Side::Upper;
143  dimension = 1;
144  }
145  if (Dim > 2 and in_lower_neighbor(k)) {
146  dimension = 2;
147  }
148  if (Dim > 2 and in_upper_neighbor(k, 2)) {
149  side = Side::Upper;
150  dimension = 2;
151  }
152 
153  const Direction<Dim> direction{dimension, side};
154  ghost_cells.at(direction).push_back(vars_with_ghosts[collapsed_index(
155  volume_with_ghosts_index, extents_with_ghosts)]);
156  }
157  }
158  }
159  }
160 
161  const size_t reconstructed_num_pts =
162  (extents[0] + 1) * extents.slice_away(0).product();
163  std::array<std::vector<double>, Dim> reconstructed_upper_side_of_face_vars =
164  make_array<Dim>(
165  std::vector<double>(reconstructed_num_pts * number_of_vars));
166  std::array<std::vector<double>, Dim> reconstructed_lower_side_of_face_vars =
167  make_array<Dim>(
168  std::vector<double>(reconstructed_num_pts * number_of_vars));
169 
170  std::array<gsl::span<double>, Dim> recons_upper_side_of_face{};
171  std::array<gsl::span<double>, Dim> recons_lower_side_of_face{};
172  for (size_t i = 0; i < Dim; ++i) {
173  auto& upper = gsl::at(reconstructed_upper_side_of_face_vars, i);
174  auto& lower = gsl::at(reconstructed_lower_side_of_face_vars, i);
175  gsl::at(recons_upper_side_of_face, i) =
176  gsl::make_span(upper.data(), upper.size());
177  gsl::at(recons_lower_side_of_face, i) =
178  gsl::make_span(lower.data(), lower.size());
179  }
180 
182  for (const auto& [direction, data] : ghost_cells) {
183  ghost_cell_vars[direction] = gsl::make_span(data.data(), data.size());
184  }
185 
186  recons_function(make_not_null(&recons_upper_side_of_face),
187  make_not_null(&recons_lower_side_of_face),
188  gsl::make_span(volume_vars), ghost_cell_vars, extents,
189  number_of_vars);
190 
191  const std::vector<size_t> extents_with_ghosts_vector(
192  extents_with_ghosts.begin(), extents_with_ghosts.end());
193  const std::vector<size_t> ghost_zones(Dim, ghost_zone_size);
194  for (size_t var_index = 0; var_index < number_of_vars; ++var_index) {
195  CAPTURE(var_index);
196  const std::vector<double> var(
197  vars_with_ghosts.begin() +
198  static_cast<std::ptrdiff_t>(var_index *
199  extents_with_ghosts.product()),
200  vars_with_ghosts.begin() +
201  static_cast<std::ptrdiff_t>((var_index + 1) *
202  extents_with_ghosts.product()));
203  const auto result =
204  pypp::call<std::vector<std::vector<std::vector<double>>>>(
205  python_test_file, python_function, var, extents_with_ghosts_vector,
206  Dim);
207 
208  const std::vector<std::vector<double>>& python_recons_on_lower = result[0];
209  const std::vector<std::vector<double>>& python_recons_on_upper = result[1];
210  for (size_t d = 0; d < Dim; ++d) {
211  CAPTURE(d);
212  const std::vector<double> recons_upper_side_this_var(
213  gsl::at(recons_upper_side_of_face, d).begin() +
214  var_index * reconstructed_num_pts,
215  gsl::at(recons_upper_side_of_face, d).begin() +
216  (var_index + 1) * reconstructed_num_pts);
217  CHECK_ITERABLE_APPROX(recons_upper_side_this_var,
218  python_recons_on_upper[d]);
219 
220  const std::vector<double> recons_lower_side_this_var(
221  gsl::at(recons_lower_side_of_face, d).begin() +
222  var_index * reconstructed_num_pts,
223  gsl::at(recons_lower_side_of_face, d).begin() +
224  (var_index + 1) * reconstructed_num_pts);
225  CHECK_ITERABLE_APPROX(recons_lower_side_this_var,
226  python_recons_on_lower[d]);
227 
228  // Test fd::reconstruction::reconstruct_neighbor
229  Index<Dim> ghost_data_extents = extents;
230  ghost_data_extents[d] = (stencil_width + 1) / 2;
231  Index<Dim> extents_with_faces = extents;
232  ++extents_with_faces[d];
233  const DataVector dv_var{
234  const_cast<double*>(volume_vars.data() +
235  var_index * extents.product()),
236  extents.product()};
237 
238  const Direction<Dim> upper_direction{d, Side::Upper};
239  DataVector upper_face_var{extents.slice_away(d).product()};
240  auto& upper_neighbor_data = ghost_cells.at(upper_direction);
241  DataVector upper_neighbor_var{
242  // NOLINTNEXTLINE
243  upper_neighbor_data.data() + var_index * ghost_data_extents.product(),
244  ghost_data_extents.product()};
245  invoke_reconstruct_neighbor(make_not_null(&upper_face_var), dv_var,
246  upper_neighbor_var, extents,
247  ghost_data_extents, upper_direction);
248  for (SliceIterator si(extents_with_faces, d, extents_with_faces[d] - 1);
249  si; ++si) {
250  CHECK(approx(upper_face_var[si.slice_offset()]) ==
251  recons_upper_side_this_var[si.volume_offset()]);
252  }
253 
254  const Direction<Dim> lower_direction{d, Side::Lower};
255  DataVector lower_face_var{extents.slice_away(d).product()};
256  auto& lower_neighbor_data = ghost_cells.at(lower_direction);
257  DataVector lower_neighbor_var{
258  // NOLINTNEXTLINE
259  lower_neighbor_data.data() + var_index * ghost_data_extents.product(),
260  ghost_data_extents.product()};
261  invoke_reconstruct_neighbor(make_not_null(&lower_face_var), dv_var,
262  lower_neighbor_var, extents,
263  ghost_data_extents, lower_direction);
264  for (SliceIterator si(extents_with_faces, d, 0); si; ++si) {
265  CHECK(approx(lower_face_var[si.slice_offset()]) ==
266  recons_lower_side_this_var[si.volume_offset()]);
267  }
268  }
269  }
270 }
271 } // 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
std::string
vector
MakeArray.hpp
TestingFramework.hpp
Side.hpp
random
fill_with_random_values
void fill_with_random_values(const gsl::not_null< T * > data, const gsl::not_null< UniformRandomBitGenerator * > generator, const gsl::not_null< RandomNumberDistribution * > distribution) noexcept
Fill an existing data structure with random values.
Definition: MakeWithRandomValues.hpp:165
MakeWithRandomValues.hpp
Index
Definition: Index.hpp:31
collapsed_index
size_t collapsed_index(const Index< N > &index, const Index< N > &extents) noexcept
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
Direction
Definition: Direction.hpp:23
std::uniform_real_distribution
TestHelpers.hpp
Pypp.hpp
cstddef
SliceIterator
Iterate over a (dim-1)-dimensional slice.
Definition: SliceIterator.hpp:24
array
Index.hpp
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:46
MAKE_GENERATOR
#define MAKE_GENERATOR(...)
MAKE_GENERATOR(NAME [, SEED]) declares a variable of name NAME containing a generator of type std::mt...
Definition: TestHelpers.hpp:433
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
Side
Side
Definition: Side.hpp:17
Gsl.hpp
std::begin
T begin(T... args)
std::ptrdiff_t
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
Index::slice_away
Index< Dim - 1 > slice_away(const size_t d) const noexcept
Return a smaller Index with the d-th element removed.
Definition: Index.hpp:87
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
string