BoundaryCorrections.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
7 
8 #include <cstddef>
9 #include <optional>
10 #include <random>
11 #include <string>
12 #include <type_traits>
13 
14 #include "DataStructures/DataBox/PrefixHelpers.hpp"
16 #include "DataStructures/DataBox/TagName.hpp"
17 #include "DataStructures/DataVector.hpp"
19 #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
22 #include "Framework/Pypp.hpp"
25 #include "Helpers/Evolution/DiscontinuousGalerkin/NormalVectors.hpp"
26 #include "Helpers/Evolution/DiscontinuousGalerkin/Range.hpp"
27 #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
29 #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
30 #include "Utilities/Gsl.hpp"
31 #include "Utilities/TMPL.hpp"
32 #include "Utilities/TaggedTuple.hpp"
33 #include "Utilities/TypeTraits/CreateHasTypeAlias.hpp"
34 
35 namespace TestHelpers::evolution::dg {
36 /// Indicate if the boundary correction should be zero when the solution is
37 /// smooth (should pretty much always be `Yes`)
38 enum class ZeroOnSmoothSolution { Yes, No };
39 
40 namespace detail {
41 template <bool HasPrimitiveVars = false>
42 struct get_correction_primitive_vars_impl {
43  template <typename BoundaryCorrection>
44  using f = tmpl::list<>;
45 };
46 
47 template <>
48 struct get_correction_primitive_vars_impl<true> {
49  template <typename BoundaryCorrection>
50  using f = typename BoundaryCorrection::dg_package_data_primitive_tags;
51 };
52 
53 template <bool HasPrimitiveVars, typename BoundaryCorrection>
54 using get_correction_primitive_vars =
55  typename get_correction_primitive_vars_impl<HasPrimitiveVars>::template f<
56  BoundaryCorrection>;
57 
58 template <bool HasPrimitiveVars = false>
59 struct get_system_primitive_vars_impl {
60  template <typename System>
61  using f = tmpl::list<>;
62 };
63 
64 template <>
65 struct get_system_primitive_vars_impl<true> {
66  template <typename System>
67  using f = typename System::primitive_variables_tag::tags_list;
68 };
69 
70 template <bool HasPrimitiveVars, typename System>
71 using get_system_primitive_vars = typename get_system_primitive_vars_impl<
72  HasPrimitiveVars>::template f<System>;
73 
74 template <typename BoundaryCorrection, typename... PackageTags,
75  typename... FaceTags, typename... VolumeTags,
76  typename... FaceTagsToForward, size_t Dim>
77 double call_dg_package_data(
78  const gsl::not_null<Variables<tmpl::list<PackageTags...>>*> package_data,
79  const BoundaryCorrection& correction,
80  const Variables<tmpl::list<FaceTags...>>& face_variables,
81  const tuples::TaggedTuple<VolumeTags...>& volume_data,
82  const tnsr::i<DataVector, Dim, Frame::Inertial>& unit_normal_covector,
83  const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
84  mesh_velocity,
85  tmpl::list<FaceTagsToForward...> /*meta*/) {
86  std::optional<Scalar<DataVector>> normal_dot_mesh_velocity{};
87  if (mesh_velocity.has_value()) {
88  normal_dot_mesh_velocity =
89  dot_product(*mesh_velocity, unit_normal_covector);
90  }
91  const double max_speed = correction.dg_package_data(
92  make_not_null(&get<PackageTags>(*package_data))...,
93  get<FaceTagsToForward>(face_variables)..., unit_normal_covector,
94  mesh_velocity, normal_dot_mesh_velocity, get<VolumeTags>(volume_data)...);
95  return max_speed;
96 }
97 
98 template <typename BoundaryCorrection, typename... PackageTags,
99  typename... FaceTags, typename... VolumeTags,
100  typename... FaceTagsToForward, size_t Dim>
101 double call_dg_package_data(
102  const gsl::not_null<Variables<tmpl::list<PackageTags...>>*> package_data,
103  const BoundaryCorrection& correction,
104  const Variables<tmpl::list<FaceTags...>>& face_variables,
105  const tuples::TaggedTuple<VolumeTags...>& volume_data,
106  const tnsr::i<DataVector, Dim, Frame::Inertial>& unit_normal_covector,
107  const tnsr::I<DataVector, Dim, Frame::Inertial>& unit_normal_vector,
108  const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
109  mesh_velocity,
110  tmpl::list<FaceTagsToForward...> /*meta*/) {
111  std::optional<Scalar<DataVector>> normal_dot_mesh_velocity{};
112  if (mesh_velocity.has_value()) {
113  normal_dot_mesh_velocity =
114  dot_product(*mesh_velocity, unit_normal_covector);
115  }
116  const double max_speed = correction.dg_package_data(
117  make_not_null(&get<PackageTags>(*package_data))...,
118  get<FaceTagsToForward>(face_variables)..., unit_normal_covector,
119  unit_normal_vector, mesh_velocity, normal_dot_mesh_velocity,
120  get<VolumeTags>(volume_data)...);
121  return max_speed;
122 }
123 
124 template <typename BoundaryCorrection, typename... BoundaryCorrectionTags,
125  typename... PackageTags>
126 void call_dg_boundary_terms(
127  const gsl::not_null<Variables<tmpl::list<BoundaryCorrectionTags...>>*>
128  boundary_corrections,
129  const BoundaryCorrection& correction,
130  const Variables<tmpl::list<PackageTags...>>& interior_package_data,
131  const Variables<tmpl::list<PackageTags...>>& exterior_package_data,
132  const ::dg::Formulation dg_formulation) {
133  correction.dg_boundary_terms(
134  make_not_null(&get<BoundaryCorrectionTags>(*boundary_corrections))...,
135  get<PackageTags>(interior_package_data)...,
136  get<PackageTags>(exterior_package_data)..., dg_formulation);
137 }
138 
139 template <typename System, typename BoundaryCorrection, size_t FaceDim,
140  typename... VolumeTags, typename... RangeTags>
141 void test_boundary_correction_conservation_impl(
142  const gsl::not_null<std::mt19937*> generator,
143  const BoundaryCorrection& correction_in, const Mesh<FaceDim>& face_mesh,
144  const tuples::TaggedTuple<VolumeTags...>& volume_data,
145  const tuples::TaggedTuple<Tags::Range<RangeTags>...>& ranges,
146  const bool use_moving_mesh, const ::dg::Formulation dg_formulation,
147  const ZeroOnSmoothSolution zero_on_smooth_solution, const double eps) {
148  CAPTURE(use_moving_mesh);
149  CAPTURE(dg_formulation);
150  CAPTURE(FaceDim);
151  constexpr bool curved_background =
152  detail::has_inverse_spatial_metric_tag_v<System>;
153  CAPTURE(curved_background);
154  Approx custom_approx = Approx::custom().epsilon(eps).scale(1.0);
155 
156  using variables_tags = typename System::variables_tag::tags_list;
157  using primitive_variables_tags = detail::get_system_primitive_vars<
158  System::has_primitive_and_conservative_vars, System>;
159  using flux_variables = typename System::flux_variables;
160  using flux_tags =
163  using temporary_tags =
164  typename System::compute_volume_time_derivative_terms::temporary_tags;
165  using dt_variables_tags = db::wrap_tags_in<::Tags::dt, variables_tags>;
166 
167  using dg_package_field_tags =
168  typename BoundaryCorrection::dg_package_field_tags;
169  using package_temporary_tags =
170  typename BoundaryCorrection::dg_package_data_temporary_tags;
171  using package_primitive_tags = detail::get_correction_primitive_vars<
172  System::has_primitive_and_conservative_vars, BoundaryCorrection>;
173 
174  const auto correction_base_ptr =
175  serialize_and_deserialize(correction_in.get_clone());
176  const auto& correction =
177  dynamic_cast<const BoundaryCorrection&>(*correction_base_ptr);
178 
179  // Check that the temporary tags needed on the boundary
180  // (package_temporary_tags) are listed as temporary tags for the volume time
181  // derivative computation (temporary_tags).
182  static_assert(
183  std::is_same_v<
185  tmpl::list<>>,
186  "There are temporary tags needed by the boundary correction that are not "
187  "computed as temporary tags by the system");
188 
189  // Check that the primitive tags needed on the boundary
190  // (package_primitive_tags) are listed as the primitive tags in
191  // the system (primitive_variables_tags).
192  static_assert(
193  std::is_same_v<tmpl::list_difference<package_primitive_tags,
194  primitive_variables_tags>,
195  tmpl::list<>>,
196  "There are primitive tags needed by the boundary correction that are not "
197  "listed in the system as being primitive variables");
198 
199  using face_tags =
200  tmpl::append<variables_tags, flux_tags, package_temporary_tags,
201  package_primitive_tags>;
202  using face_tags_with_curved_background = tmpl::conditional_t<
203  curved_background,
204  tmpl::remove_duplicates<tmpl::push_back<
205  face_tags, typename detail::inverse_spatial_metric_tag<
206  curved_background>::template f<System>>>,
207  face_tags>;
208 
209  // Sanity check: make sure each range tag is a `dg_package_data` argument
210  using invalid_range_tags =
211  tmpl::list_difference<tmpl::list<RangeTags...>,
212  face_tags_with_curved_background>;
213  static_assert(std::is_same_v<tmpl::list<>, invalid_range_tags>,
214  "Received Tags::Range for Tags that aren't arguments to "
215  "dg_package_data");
216 
217  std::uniform_real_distribution<> dist(-1.0, 1.0);
218  DataVector used_for_size{face_mesh.number_of_grid_points()};
219 
220  // Fill all fields with random values in [-1,1), then, for each tag with a
221  // specified range, overwrite with new random values in [min,max)
222  Variables<dg_package_field_tags> interior_package_data{used_for_size.size()};
223  auto interior_fields_on_face =
224  make_with_random_values<Variables<face_tags_with_curved_background>>(
225  generator, make_not_null(&dist), used_for_size);
226  tmpl::for_each<tmpl::list<RangeTags...>>([&generator,
227  &interior_fields_on_face,
228  &ranges](auto tag_v) {
229  using tag = tmpl::type_from<decltype(tag_v)>;
230  const std::array<double, 2>& range = tuples::get<Tags::Range<tag>>(ranges);
231  std::uniform_real_distribution<> local_dist(range[0], range[1]);
232  fill_with_random_values(make_not_null(&get<tag>(interior_fields_on_face)),
233  generator, make_not_null(&local_dist));
234  });
235 
236  // Same as above but now for external data
237  Variables<dg_package_field_tags> exterior_package_data{used_for_size.size()};
238  auto exterior_fields_on_face =
239  make_with_random_values<Variables<face_tags_with_curved_background>>(
240  generator, make_not_null(&dist), used_for_size);
241  tmpl::for_each<tmpl::list<RangeTags...>>([&generator,
242  &exterior_fields_on_face,
243  &ranges](auto tag_v) {
244  using tag = tmpl::type_from<decltype(tag_v)>;
245  const std::array<double, 2>& range = tuples::get<Tags::Range<tag>>(ranges);
246  std::uniform_real_distribution<> local_dist(range[0], range[1]);
247  fill_with_random_values(make_not_null(&get<tag>(exterior_fields_on_face)),
248  generator, make_not_null(&local_dist));
249  });
250 
251  // Compute the interior and exterior normal vectors so they are pointing in
252  // opposite directions.
253  auto interior_unit_normal_covector = make_with_random_values<
254  tnsr::i<DataVector, FaceDim + 1, Frame::Inertial>>(
255  generator, make_not_null(&dist), used_for_size);
256  tnsr::I<DataVector, FaceDim + 1, Frame::Inertial>
257  interior_unit_normal_vector{};
258 
259  tnsr::i<DataVector, FaceDim + 1, Frame::Inertial>
260  exterior_unit_normal_covector;
261  tnsr::I<DataVector, FaceDim + 1, Frame::Inertial>
262  exterior_unit_normal_vector{};
263  if constexpr (not curved_background) {
264  const Scalar<DataVector> interior_normal_magnitude =
265  magnitude(interior_unit_normal_covector);
266  for (auto& t : interior_unit_normal_covector) {
267  t /= get(interior_normal_magnitude);
268  }
269  exterior_unit_normal_covector = interior_unit_normal_covector;
270  for (auto& t : exterior_unit_normal_covector) {
271  t *= -1.0;
272  }
273  } else {
274  using inv_spatial_metric = typename detail::inverse_spatial_metric_tag<
275  curved_background>::template f<System>;
276  exterior_unit_normal_covector = interior_unit_normal_covector;
277  for (auto& t : exterior_unit_normal_covector) {
278  t *= -1.0;
279  }
280  detail::adjust_spatial_metric_or_inverse(
281  make_not_null(&get<inv_spatial_metric>(interior_fields_on_face)));
282 
283  // make the exterior inverse spatial metric be close to the interior one
284  // since the solution should be smooth. We can't change too much because we
285  // want the spatial metric to have diagonal terms much larger than the
286  // off-diagonal terms.
287  std::uniform_real_distribution<> inv_metric_change_dist(0.999, 1.0);
288  for (size_t i = 0; i < FaceDim + 1; ++i) {
289  for (size_t j = i; j < FaceDim + 1; ++j) {
290  get<inv_spatial_metric>(exterior_fields_on_face).get(i, j) =
291  inv_metric_change_dist(*generator) *
292  get<inv_spatial_metric>(interior_fields_on_face).get(i, j);
293  }
294  }
295  detail::normalize_vector_and_covector(
296  make_not_null(&interior_unit_normal_covector),
297  make_not_null(&interior_unit_normal_vector),
298  get<inv_spatial_metric>(interior_fields_on_face));
299  detail::normalize_vector_and_covector(
300  make_not_null(&exterior_unit_normal_covector),
301  make_not_null(&exterior_unit_normal_vector),
302  get<inv_spatial_metric>(exterior_fields_on_face));
303 
304  // if dg_package_data depends on the (not inverted) spatial metric, we also
305  // adjust the spatial metric. Note that for testing purposes, the metric and
306  // inverse metric don't need to be inverses of each other. However, we do
307  // want both to be physically reasonable, i.e., close to flat space.
308  using spatial_metric =
310  if constexpr (tmpl::list_contains_v<face_tags_with_curved_background,
311  spatial_metric>) {
312  detail::adjust_spatial_metric_or_inverse(
313  make_not_null(&get<spatial_metric>(interior_fields_on_face)));
314  // as for external inverse metric: we get the external metric by slightly
315  // tweaking the internal metric
316  for (size_t i = 0; i < FaceDim + 1; ++i) {
317  for (size_t j = i; j < FaceDim + 1; ++j) {
318  get<spatial_metric>(exterior_fields_on_face).get(i, j) =
319  inv_metric_change_dist(*generator) *
320  get<spatial_metric>(interior_fields_on_face).get(i, j);
321  }
322  }
323  }
324  }
325 
327  mesh_velocity{};
328  if (use_moving_mesh) {
329  mesh_velocity = make_with_random_values<
330  tnsr::I<DataVector, FaceDim + 1, Frame::Inertial>>(
331  generator, make_not_null(&dist), used_for_size);
332  }
333 
334  if constexpr (curved_background) {
335  call_dg_package_data(
336  make_not_null(&interior_package_data), correction,
337  interior_fields_on_face, volume_data, interior_unit_normal_covector,
338  interior_unit_normal_vector, mesh_velocity, face_tags{});
339  call_dg_package_data(
340  make_not_null(&exterior_package_data), correction,
341  exterior_fields_on_face, volume_data, exterior_unit_normal_covector,
342  exterior_unit_normal_vector, mesh_velocity, face_tags{});
343  } else {
344  call_dg_package_data(make_not_null(&interior_package_data), correction,
345  interior_fields_on_face, volume_data,
346  interior_unit_normal_covector, mesh_velocity,
347  face_tags{});
348  call_dg_package_data(make_not_null(&exterior_package_data), correction,
349  exterior_fields_on_face, volume_data,
350  exterior_unit_normal_covector, mesh_velocity,
351  face_tags{});
352  }
353 
354  Variables<dt_variables_tags> boundary_corrections{
355  face_mesh.number_of_grid_points()};
356  call_dg_boundary_terms(make_not_null(&boundary_corrections), correction,
357  interior_package_data, exterior_package_data,
358  dg_formulation);
359 
360  if (dg_formulation == ::dg::Formulation::StrongInertial) {
361  // The strong form should be (WeakForm - (n_i F^i)_{interior}).
362  // Since we also test conservation for the weak form we just need to test
363  // that the strong form satisfies the above definition.
364 
365  if constexpr (curved_background) {
366  call_dg_package_data(
367  make_not_null(&interior_package_data), correction,
368  interior_fields_on_face, volume_data, interior_unit_normal_covector,
369  interior_unit_normal_vector, mesh_velocity, face_tags{});
370  call_dg_package_data(
371  make_not_null(&exterior_package_data), correction,
372  exterior_fields_on_face, volume_data, exterior_unit_normal_covector,
373  exterior_unit_normal_vector, mesh_velocity, face_tags{});
374  } else {
375  call_dg_package_data(make_not_null(&interior_package_data), correction,
376  interior_fields_on_face, volume_data,
377  interior_unit_normal_covector, mesh_velocity,
378  face_tags{});
379  call_dg_package_data(make_not_null(&exterior_package_data), correction,
380  exterior_fields_on_face, volume_data,
381  exterior_unit_normal_covector, mesh_velocity,
382  face_tags{});
383  }
384 
385  Variables<dt_variables_tags> expected_boundary_corrections{
386  face_mesh.number_of_grid_points()};
387  call_dg_boundary_terms(make_not_null(&expected_boundary_corrections),
388  correction, interior_package_data,
389  exterior_package_data,
390  ::dg::Formulation::WeakInertial);
391 
392  tmpl::for_each<flux_variables>([&interior_package_data,
393  &expected_boundary_corrections](
394  auto flux_variable_tag_v) noexcept {
395  using flux_variable_tag = tmpl::type_from<decltype(flux_variable_tag_v)>;
396  using normal_dot_flux_tag = ::Tags::NormalDotFlux<flux_variable_tag>;
397  using dt_tag = ::Tags::dt<flux_variable_tag>;
398  const auto& normal_dot_flux =
399  get<normal_dot_flux_tag>(interior_package_data);
400  auto& expected_boundary_correction =
401  get<dt_tag>(expected_boundary_corrections);
402  for (size_t tensor_index = 0;
403  tensor_index < expected_boundary_correction.size(); ++tensor_index) {
404  expected_boundary_correction[tensor_index] -=
405  normal_dot_flux[tensor_index];
406  }
407  });
408  {
409  INFO("Check weak and strong boundary terms match.");
411  boundary_corrections, expected_boundary_corrections, custom_approx);
412  }
413 
414  if (zero_on_smooth_solution == ZeroOnSmoothSolution::Yes) {
415  INFO(
416  "Testing that if the solution is the same on both sides the "
417  "StrongInertial correction is identically zero.");
418  Variables<dg_package_field_tags> interior_package_data_opposite_signs{
419  used_for_size.size()};
420  if constexpr (curved_background) {
421  // If the solution is the same on both sides then the interior and
422  // exterior normal vectors only differ by a sign.
423  for (size_t i = 0; i < FaceDim + 1; ++i) {
424  exterior_unit_normal_vector.get(i) =
425  -interior_unit_normal_vector.get(i);
426  exterior_unit_normal_covector.get(i) =
427  -interior_unit_normal_covector.get(i);
428  }
429  call_dg_package_data(
430  make_not_null(&interior_package_data_opposite_signs), correction,
431  interior_fields_on_face, volume_data, exterior_unit_normal_covector,
432  exterior_unit_normal_vector, mesh_velocity, face_tags{});
433  } else {
434  call_dg_package_data(
435  make_not_null(&interior_package_data_opposite_signs), correction,
436  interior_fields_on_face, volume_data, exterior_unit_normal_covector,
437  mesh_velocity, face_tags{});
438  }
439  Variables<dt_variables_tags> zero_boundary_correction{
440  face_mesh.number_of_grid_points()};
441  call_dg_boundary_terms(make_not_null(&zero_boundary_correction),
442  correction, interior_package_data,
443  interior_package_data_opposite_signs,
444  ::dg::Formulation::StrongInertial);
445  Variables<dt_variables_tags> expected_zero_boundary_correction{
446  face_mesh.number_of_grid_points(), 0.0};
447  tmpl::for_each<dt_variables_tags>([&custom_approx,
448  &expected_zero_boundary_correction,
449  &zero_boundary_correction](
450  auto dt_variables_tag_v) noexcept {
451  using dt_variables_tag = tmpl::type_from<decltype(dt_variables_tag_v)>;
452  const std::string tag_name = db::tag_name<dt_variables_tag>();
453  CAPTURE(tag_name);
455  get<dt_variables_tag>(zero_boundary_correction),
456  get<dt_variables_tag>(expected_zero_boundary_correction),
457  custom_approx);
458  });
459  }
460  } else if (dg_formulation == ::dg::Formulation::WeakInertial) {
461  INFO(
462  "Checking that swapping the two sides results in an overall minus "
463  "sign.");
464  Variables<dt_variables_tags> reverse_side_boundary_corrections{
465  face_mesh.number_of_grid_points()};
466  call_dg_boundary_terms(make_not_null(&reverse_side_boundary_corrections),
467  correction, exterior_package_data,
468  interior_package_data, dg_formulation);
469  // Check that the flux leaving one element equals the flux entering its
470  // neighbor, i.e., F*(interior, exterior) == -F*(exterior, interior)
471  reverse_side_boundary_corrections *= -1.0;
472  tmpl::for_each<flux_variables>([&boundary_corrections, &custom_approx,
473  &reverse_side_boundary_corrections](
474  auto flux_variable_tag_v) {
475  using flux_variable_tag = tmpl::type_from<decltype(flux_variable_tag_v)>;
476  const std::string tag_name = db::tag_name<flux_variable_tag>();
477  CAPTURE(tag_name);
479  get<::Tags::dt<flux_variable_tag>>(reverse_side_boundary_corrections),
480  get<::Tags::dt<flux_variable_tag>>(boundary_corrections),
481  custom_approx);
482  });
483  } else {
484  ERROR("DG formulation must be StrongInertial or WeakInertial, not "
485  << dg_formulation);
486  }
487 }
488 } // namespace detail
489 
490 /*!
491  * \ingroup TestingFrameworkGroup
492  * \brief Checks that the boundary correction is conservative and that for
493  * smooth solutions the strong-form correction is zero.
494  *
495  * By default, each input tensor for `dg_package_data` is randomly generated
496  * from the interval `[-1,1)` (except the metric, which for systems with a
497  * metric is generated to be close to flat space). The argument `ranges` is a
498  * `TaggedTuple` of `TestHelpers::evolution::dg::Tags::Range<tag>` that
499  * enables the caller to pick a custom interval for generating the input
500  * `tag`. This is useful, for example, for tensors that need positive values.
501  * Each `tag` in `ranges` must be an argument of `dg_package_data`.
502  */
503 template <typename System, typename BoundaryCorrection, size_t FaceDim,
504  typename... VolumeTags, typename... RangeTags>
506  const gsl::not_null<std::mt19937*> generator,
507  const BoundaryCorrection& correction, const Mesh<FaceDim>& face_mesh,
508  const tuples::TaggedTuple<VolumeTags...>& volume_data,
510  const ZeroOnSmoothSolution zero_on_smooth_solution =
511  ZeroOnSmoothSolution::Yes,
512  const double eps = 1.0e-12) {
513  for (const auto use_moving_mesh : {true, false}) {
514  for (const auto& dg_formulation :
515  {::dg::Formulation::StrongInertial, ::dg::Formulation::WeakInertial}) {
516  detail::test_boundary_correction_conservation_impl<System>(
517  generator, correction, face_mesh, volume_data, ranges,
518  use_moving_mesh, dg_formulation, zero_on_smooth_solution, eps);
519  }
520  }
521 }
522 
523 namespace detail {
524 template <typename System, typename ConversionClassList, typename VariablesTags,
525  typename BoundaryCorrection, size_t FaceDim, typename... FaceTags,
526  typename... VolumeTags, typename... RangeTags,
527  typename... DgPackageDataTags>
528 void test_with_python(
529  const gsl::not_null<std::mt19937*> generator,
530  const std::string& python_module,
531  const std::array<
532  std::string,
533  tmpl::size<typename BoundaryCorrection::dg_package_field_tags>::value>&
534  python_dg_package_data_functions,
535  const std::array<std::string, tmpl::size<VariablesTags>::value>&
536  python_dg_boundary_terms_functions,
537  const BoundaryCorrection& correction, const Mesh<FaceDim>& face_mesh,
538  const tuples::TaggedTuple<VolumeTags...>& volume_data,
539  const tuples::TaggedTuple<Tags::Range<RangeTags>...>& ranges,
540  const bool use_moving_mesh, const ::dg::Formulation dg_formulation,
541  const double epsilon, tmpl::list<FaceTags...> /*meta*/,
542  tmpl::list<DgPackageDataTags...> /*meta*/) {
543  CAPTURE(face_mesh);
544  CAPTURE(dg_formulation);
545  CAPTURE(use_moving_mesh);
546  REQUIRE(face_mesh.number_of_grid_points() >= 1);
547  constexpr bool curved_background =
548  detail::has_inverse_spatial_metric_tag_v<System>;
549  using dg_package_field_tags =
550  typename BoundaryCorrection::dg_package_field_tags;
551 
552  using face_tags = tmpl::list<FaceTags...>;
553  using face_tags_with_curved_background = tmpl::conditional_t<
554  curved_background,
555  tmpl::remove_duplicates<tmpl::push_back<
556  face_tags, typename TestHelpers::evolution::dg::detail::
557  inverse_spatial_metric_tag<
558  curved_background>::template f<System>>>,
559  face_tags>;
560 
561  // Sanity check: make sure each range tag is a `dg_package_data` argument
562  using invalid_range_tags =
563  tmpl::list_difference<tmpl::list<RangeTags...>,
564  face_tags_with_curved_background>;
565  static_assert(std::is_same_v<tmpl::list<>, invalid_range_tags>,
566  "Received Tags::Range for Tags that aren't arguments to "
567  "dg_package_data");
568 
569  std::uniform_real_distribution<> dist(-1.0, 1.0);
570  DataVector used_for_size{face_mesh.number_of_grid_points()};
571 
572  // Fill all fields with random values in [-1,1), then, for each tag with a
573  // specified range, overwrite with new random values in [min,max)
574  auto fields_on_face =
575  make_with_random_values<Variables<face_tags_with_curved_background>>(
576  generator, make_not_null(&dist), used_for_size);
577  tmpl::for_each<tmpl::list<RangeTags...>>([&generator, &fields_on_face,
578  &ranges](auto tag_v) {
579  using tag = tmpl::type_from<decltype(tag_v)>;
580  const std::array<double, 2>& range = tuples::get<Tags::Range<tag>>(ranges);
581  std::uniform_real_distribution<> local_dist(range[0], range[1]);
582  fill_with_random_values(make_not_null(&get<tag>(fields_on_face)), generator,
583  make_not_null(&local_dist));
584  });
585 
586  auto unit_normal_covector = make_with_random_values<
587  tnsr::i<DataVector, FaceDim + 1, Frame::Inertial>>(
588  generator, make_not_null(&dist), used_for_size);
589  tnsr::I<DataVector, FaceDim + 1, Frame::Inertial> unit_normal_vector{};
590 
591  if constexpr (not curved_background) {
592  const Scalar<DataVector> normal_magnitude = magnitude(unit_normal_covector);
593  for (auto& t : unit_normal_covector) {
594  t /= get(normal_magnitude);
595  }
596  } else {
597  using inv_spatial_metric = typename detail::inverse_spatial_metric_tag<
598  curved_background>::template f<System>;
599  detail::adjust_spatial_metric_or_inverse(
600  make_not_null(&get<inv_spatial_metric>(fields_on_face)));
601  detail::normalize_vector_and_covector(
602  make_not_null(&unit_normal_covector),
603  make_not_null(&unit_normal_vector),
604  get<inv_spatial_metric>(fields_on_face));
605 
606  using spatial_metric =
608  if constexpr (tmpl::list_contains_v<face_tags_with_curved_background,
609  spatial_metric>) {
610  detail::adjust_spatial_metric_or_inverse(
611  make_not_null(&get<spatial_metric>(fields_on_face)));
612  }
613  }
614 
616  mesh_velocity{};
617  if (use_moving_mesh) {
618  mesh_velocity = make_with_random_values<
619  tnsr::I<DataVector, FaceDim + 1, Frame::Inertial>>(
620  generator, make_not_null(&dist), used_for_size);
621  }
622  std::optional<Scalar<DataVector>> normal_dot_mesh_velocity{};
623  if (mesh_velocity.has_value()) {
624  normal_dot_mesh_velocity =
625  dot_product(*mesh_velocity, unit_normal_covector);
626  }
627 
628  // Call C++ implementation of dg_package_data
629  Variables<dg_package_field_tags> package_data{used_for_size.size()};
630  if constexpr (curved_background) {
631  call_dg_package_data(make_not_null(&package_data), correction,
632  fields_on_face, volume_data, unit_normal_covector,
633  unit_normal_vector, mesh_velocity,
634  tmpl::list<FaceTags...>{});
635  } else {
636  call_dg_package_data(make_not_null(&package_data), correction,
637  fields_on_face, volume_data, unit_normal_covector,
638  mesh_velocity, tmpl::list<FaceTags...>{});
639  }
640 
641  // Call python implementation of dg_package_data
642  size_t function_name_index = 0;
643  tmpl::for_each<dg_package_field_tags>(
644  [epsilon, &fields_on_face, &function_name_index, &mesh_velocity,
645  &normal_dot_mesh_velocity, &package_data,
646  &python_dg_package_data_functions, &python_module, &unit_normal_covector,
647  &unit_normal_vector, &volume_data](auto package_data_tag_v) {
648  // avoid compiler warnings if there isn't any volume data
649  (void)volume_data;
650  INFO("Testing package data");
651  using package_data_tag = tmpl::type_from<decltype(package_data_tag_v)>;
652  using ResultType = typename package_data_tag::type;
653  const std::string tag_name = db::tag_name<package_data_tag>();
654  CAPTURE(tag_name);
655  try {
656  CAPTURE(python_module);
657  CAPTURE(
658  gsl::at(python_dg_package_data_functions, function_name_index));
659  CAPTURE(pretty_type::short_name<ResultType>());
660  if constexpr (curved_background) {
661  const auto python_result =
662  pypp::call<ResultType, ConversionClassList>(
663  python_module,
664  gsl::at(python_dg_package_data_functions,
665  function_name_index),
666  get<FaceTags>(fields_on_face)..., unit_normal_covector,
667  unit_normal_vector, mesh_velocity, normal_dot_mesh_velocity,
668  get<VolumeTags>(volume_data)...);
670  get<package_data_tag>(package_data), python_result,
671  Approx::custom().epsilon(epsilon).scale(1.0));
672  } else {
673  (void)unit_normal_vector;
674  const auto python_result =
675  pypp::call<ResultType, ConversionClassList>(
676  python_module,
677  gsl::at(python_dg_package_data_functions,
678  function_name_index),
679  get<FaceTags>(fields_on_face)..., unit_normal_covector,
680  mesh_velocity, normal_dot_mesh_velocity,
681  get<VolumeTags>(volume_data)...);
683  get<package_data_tag>(package_data), python_result,
684  Approx::custom().epsilon(epsilon).scale(1.0));
685  }
686  } catch (const std::exception& e) {
687  INFO("On line " << __LINE__ << " Python call to "
688  << gsl::at(python_dg_package_data_functions,
689  function_name_index)
690  << " in module " << python_module
691  << " failed: " << e.what());
692  REQUIRE(false);
693  }
694  ++function_name_index;
695  });
696 
697  // Now we need to check the dg_boundary_terms function.
698  //
699  // In order for the package data to have a self-consistent state (not all
700  // things being packaged are independent), we compute the packaged data.
701  // We reuse the previously computed packaged data as the interior packaged
702  // data.
703  const auto& interior_package_data = package_data;
704 
705  auto exterior_fields_on_face =
706  make_with_random_values<Variables<face_tags_with_curved_background>>(
707  generator, make_not_null(&dist), used_for_size);
708  tmpl::for_each<tmpl::list<RangeTags...>>([&generator,
709  &exterior_fields_on_face,
710  &ranges](auto tag_v) {
711  using tag = tmpl::type_from<decltype(tag_v)>;
712  const std::array<double, 2>& range = tuples::get<Tags::Range<tag>>(ranges);
713  std::uniform_real_distribution<> local_dist(range[0], range[1]);
714  fill_with_random_values(make_not_null(&get<tag>(exterior_fields_on_face)),
715  generator, make_not_null(&local_dist));
716  });
717 
718  tnsr::i<DataVector, FaceDim + 1, Frame::Inertial>
719  exterior_unit_normal_covector;
720  tnsr::I<DataVector, FaceDim + 1, Frame::Inertial>
721  exterior_unit_normal_vector{};
722  if constexpr (not curved_background) {
723  exterior_unit_normal_covector = unit_normal_covector;
724  for (auto& t : exterior_unit_normal_covector) {
725  t *= -1.0;
726  }
727  } else {
728  using inv_spatial_metric = typename detail::inverse_spatial_metric_tag<
729  curved_background>::template f<System>;
730  exterior_unit_normal_covector = unit_normal_covector;
731  for (auto& t : exterior_unit_normal_covector) {
732  t *= -1.0;
733  }
734 
735  // make the exterior inverse spatial metric be close to the interior one
736  // since the solution should be smooth. We can't change too much because we
737  // want the spatial metric to have diagonal terms much larger than the
738  // off-diagonal terms.
739  std::uniform_real_distribution<> inv_metric_change_dist(0.999, 1.0);
740  for (size_t i = 0; i < FaceDim + 1; ++i) {
741  for (size_t j = i; j < FaceDim + 1; ++j) {
742  get<inv_spatial_metric>(exterior_fields_on_face).get(i, j) =
743  inv_metric_change_dist(*generator) *
744  get<inv_spatial_metric>(fields_on_face).get(i, j);
745  }
746  }
747  detail::normalize_vector_and_covector(
748  make_not_null(&exterior_unit_normal_covector),
749  make_not_null(&exterior_unit_normal_vector),
750  get<inv_spatial_metric>(exterior_fields_on_face));
751 
752  // if dg_package_data depends on the (not inverted) spatial metric, we also
753  // adjust the spatial metric. Note that for testing purposes, the metric and
754  // inverse metric don't need to be inverses of each other. However, we do
755  // want both to be physically reasonable, i.e., close to flat space.
756  using spatial_metric =
758  if constexpr (tmpl::list_contains_v<face_tags_with_curved_background,
759  spatial_metric>) {
760  // as for external inverse metric: we get the external metric by slightly
761  // tweaking the internal metric
762  for (size_t i = 0; i < FaceDim + 1; ++i) {
763  for (size_t j = i; j < FaceDim + 1; ++j) {
764  get<spatial_metric>(exterior_fields_on_face).get(i, j) =
765  inv_metric_change_dist(*generator) *
766  get<spatial_metric>(fields_on_face).get(i, j);
767  }
768  }
769  }
770  }
771 
772  Variables<dg_package_field_tags> exterior_package_data{used_for_size.size()};
773  if constexpr (curved_background) {
774  call_dg_package_data(
775  make_not_null(&exterior_package_data), correction,
776  exterior_fields_on_face, volume_data, exterior_unit_normal_covector,
777  exterior_unit_normal_vector, mesh_velocity, face_tags{});
778  } else {
779  call_dg_package_data(make_not_null(&exterior_package_data), correction,
780  exterior_fields_on_face, volume_data,
781  exterior_unit_normal_covector, mesh_velocity,
782  face_tags{});
783  }
784 
785  // We don't need to prefix the VariablesTags with anything because we are not
786  // interacting with any code that cares about what the tags are, just that the
787  // types matched the evolved variables.
788  Variables<VariablesTags> boundary_corrections{
789  face_mesh.number_of_grid_points()};
790 
791  // Call C++ implementation of dg_boundary_terms
792  call_dg_boundary_terms(make_not_null(&boundary_corrections), correction,
793  interior_package_data, exterior_package_data,
794  dg_formulation);
795 
796  // Call python implementation of dg_boundary_terms
797  function_name_index = 0;
798  tmpl::for_each<VariablesTags>([&boundary_corrections, &dg_formulation,
799  epsilon, &exterior_package_data,
800  &function_name_index, &interior_package_data,
801  &python_dg_boundary_terms_functions,
802  &python_module](auto package_data_tag_v) {
803  INFO("Testing boundary terms");
804  using boundary_correction_tag =
805  tmpl::type_from<decltype(package_data_tag_v)>;
806  using ResultType = typename boundary_correction_tag::type;
807  try {
808  CAPTURE(python_module);
809  const std::string& python_function =
810  gsl::at(python_dg_boundary_terms_functions, function_name_index);
811  CAPTURE(python_function);
812  // Make explicitly depend on tag type to avoid type deduction issues with
813  // GCC7
814  const typename boundary_correction_tag::type python_result =
815  pypp::call<ResultType>(
816  python_module, python_function,
817  get<DgPackageDataTags>(interior_package_data)...,
818  get<DgPackageDataTags>(exterior_package_data)...,
819  dg_formulation == ::dg::Formulation::StrongInertial);
821  get<boundary_correction_tag>(boundary_corrections), python_result,
822  Approx::custom().epsilon(epsilon).scale(1.0));
823  } catch (const std::exception& e) {
824  INFO("On line " << __LINE__ << " Python call to "
825  << gsl::at(python_dg_boundary_terms_functions,
826  function_name_index)
827  << " in module " << python_module
828  << " failed: " << e.what());
829  REQUIRE(false);
830  }
831  ++function_name_index;
832  });
833 }
834 } // namespace detail
835 
836 /*!
837  * \ingroup TestingFrameworkGroup
838  * \brief Tests that the `dg_package_data` and `dg_boundary_terms` functions
839  * agree with the python implementation.
840  *
841  * The variables are filled with random numbers between zero and one before
842  * being passed to the implementations. If in the future we need support for
843  * negative numbers we can add the ability to specify a single range for all
844  * random numbers or each individually.
845  *
846  * Please note the following:
847  * - The `pypp::SetupLocalPythonEnvironment` must be created before the
848  * `test_boundary_correction_with_python` can be called.
849  * - The `dg_formulation` is passed as a bool `use_strong_form` to the python
850  * functions since we don't want to rely on python bindings for the enum.
851  * - You can convert custom types using the `ConversionClassList` template
852  * parameter, which is then passed to `pypp::call()`. This allows you to,
853  * e.g., convert an equation of state into an array locally in a test file.
854  * - There must be one python function to compute the packaged data for each tag
855  * in `dg_package_field_tags`
856  * - There must be one python function to compute the boundary correction for
857  * each tag in `System::variables_tag`
858  * - The arguments to the python functions for computing the packaged data are
859  * the same as the arguments for the C++ `dg_package_data` function, excluding
860  * the `gsl::not_null` arguments.
861  * - The arguments to the python functions for computing the boundary
862  * corrections are the same as the arguments for the C++ `dg_boundary_terms`
863  * function, excluding the `gsl::not_null` arguments.
864  * - By default, each input tensor for `dg_package_data` is randomly generated
865  * from the interval `[-1,1)` (except the metric, which for systems with a
866  * metric is generated to be close to flat space). The argument `ranges` is a
867  * `TaggedTuple` of `TestHelpers::evolution::dg::Tags::Range<tag>` that
868  * enables the caller to pick a custom interval for generating the input
869  * `tag`. This is useful, for example, for tensors that need positive values.
870  * Each `tag` in `ranges` must be an argument of `dg_package_data`.
871  */
872 template <typename System, typename ConversionClassList = tmpl::list<>,
873  typename BoundaryCorrection, size_t FaceDim, typename... VolumeTags,
874  typename... RangeTags>
876  const gsl::not_null<std::mt19937*> generator,
877  const std::string& python_module,
878  const std::array<
879  std::string,
880  tmpl::size<typename BoundaryCorrection::dg_package_field_tags>::value>&
881  python_dg_package_data_functions,
882  const std::array<
883  std::string,
884  tmpl::size<typename System::variables_tag::tags_list>::value>&
885  python_dg_boundary_terms_functions,
886  const BoundaryCorrection& correction, const Mesh<FaceDim>& face_mesh,
887  const tuples::TaggedTuple<VolumeTags...>& volume_data,
889  const double epsilon = 1.0e-12) {
890  static_assert(std::is_final_v<std::decay_t<BoundaryCorrection>>,
891  "All boundary correction classes must be marked `final`.");
892  using package_temporary_tags =
893  typename BoundaryCorrection::dg_package_data_temporary_tags;
894  using package_primitive_tags = detail::get_correction_primitive_vars<
895  System::has_primitive_and_conservative_vars, BoundaryCorrection>;
896  using variables_tags = typename System::variables_tag::tags_list;
897  using flux_variables = typename System::flux_variables;
898  using flux_tags =
901 
902  for (const auto use_moving_mesh : {false, true}) {
903  for (const auto dg_formulation :
904  {::dg::Formulation::StrongInertial, ::dg::Formulation::WeakInertial}) {
905  detail::test_with_python<System, ConversionClassList, variables_tags>(
906  generator, python_module, python_dg_package_data_functions,
907  python_dg_boundary_terms_functions, correction, face_mesh,
908  volume_data, ranges, use_moving_mesh, dg_formulation, epsilon,
909  tmpl::append<variables_tags, flux_tags, package_temporary_tags,
910  package_primitive_tags>{},
911  typename BoundaryCorrection::dg_package_field_tags{});
912  }
913  }
914 }
915 } // namespace TestHelpers::evolution::dg
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
TestHelpers::evolution::dg::test_boundary_correction_conservation
void test_boundary_correction_conservation(const gsl::not_null< std::mt19937 * > generator, const BoundaryCorrection &correction, const Mesh< FaceDim > &face_mesh, const tuples::TaggedTuple< VolumeTags... > &volume_data, const tuples::TaggedTuple< Tags::Range< RangeTags >... > &ranges, const ZeroOnSmoothSolution zero_on_smooth_solution=ZeroOnSmoothSolution::Yes, const double eps=1.0e-12)
Checks that the boundary correction is conservative and that for smooth solutions the strong-form cor...
Definition: BoundaryCorrections.hpp:505
std::string
Frame::Inertial
Definition: IndexType.hpp:44
std::exception
gr::Tags::SpatialMetric
Definition: Tags.hpp:26
TestingFramework.hpp
CHECK_VARIABLES_CUSTOM_APPROX
#define CHECK_VARIABLES_CUSTOM_APPROX(a, b, appx)
Same as CHECK_VARIABLES_APPROX, but with a user-defined Approx. The third argument should be of type ...
Definition: TestHelpers.hpp:462
db::tag_name
std::string tag_name() noexcept
Get the name of a DataBox tag, including prefixes.
Definition: TagName.hpp:41
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
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:791
Tags::NormalDotFlux
Prefix indicating a boundary unit normal vector dotted into the flux.
Definition: Prefixes.hpp:96
CHECK_ITERABLE_CUSTOM_APPROX
#define CHECK_ITERABLE_CUSTOM_APPROX(a, b, appx)
Same as CHECK_ITERABLE_APPROX with user-defined Approx. The third argument should be of type Approx.
Definition: TestingFramework.hpp:134
TestHelpers::evolution::dg::Tags::Range
Tag for a TaggedTuple that holds the range of validity for the variable associated with Tag.
Definition: Range.hpp:13
std::uniform_real_distribution
TestHelpers.hpp
dot_product
void dot_product(const gsl::not_null< Scalar< DataType > * > dot_product, const Tensor< DataType, Symmetry< 1 >, index_list< Index >> &vector_a, const Tensor< DataType, Symmetry< 1 >, index_list< Index >> &vector_b) noexcept
Compute the Euclidean dot product of two vectors or one forms.
Definition: DotProduct.hpp:25
domain::push_back
CoordinateMap< SourceFrame, TargetFrame, Maps..., NewMap > push_back(CoordinateMap< SourceFrame, TargetFrame, Maps... > old_map, NewMap new_map) noexcept
Creates a CoordinateMap by appending the new map to the end of the old maps.
ERROR
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:37
Pypp.hpp
DotProduct.hpp
cstddef
std::array< double, 2 >
Mesh::number_of_grid_points
size_t number_of_grid_points() const noexcept
The total number of grid points in all dimensions.
Definition: Mesh.hpp:165
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:46
serialize_and_deserialize
T serialize_and_deserialize(const T &t)
Serializes and deserializes an object t of type T
Definition: TestHelpers.hpp:45
std::decay_t
Variables.hpp
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:49
Tags::dt
Prefix indicating a time derivative.
Definition: Prefixes.hpp:29
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
make_with_random_values
ReturnType make_with_random_values(const gsl::not_null< UniformRandomBitGenerator * > generator, const gsl::not_null< RandomNumberDistribution * > distribution, const T &used_for_size) noexcept
Make a data structure and fill it with random values.
Definition: MakeWithRandomValues.hpp:185
Gsl.hpp
Tensor.hpp
db::wrap_tags_in
tmpl::transform< TagList, tmpl::bind< Wrapper, tmpl::_1, tmpl::pin< Args >... > > wrap_tags_in
Create a new tmpl::list of tags by wrapping each tag in TagList in Wrapper<_, Args....
Definition: PrefixHelpers.hpp:30
magnitude
Scalar< DataType > magnitude(const Tensor< DataType, Symmetry< 1 >, index_list< Index >> &vector) noexcept
Compute the Euclidean magnitude of a rank-1 tensor.
Definition: Magnitude.hpp:27
TestHelpers::evolution::dg::test_boundary_correction_with_python
void test_boundary_correction_with_python(const gsl::not_null< std::mt19937 * > generator, const std::string &python_module, const std::array< std::string, tmpl::size< typename BoundaryCorrection::dg_package_field_tags >::value > &python_dg_package_data_functions, const std::array< std::string, tmpl::size< typename System::variables_tag::tags_list >::value > &python_dg_boundary_terms_functions, const BoundaryCorrection &correction, const Mesh< FaceDim > &face_mesh, const tuples::TaggedTuple< VolumeTags... > &volume_data, const tuples::TaggedTuple< Tags::Range< RangeTags >... > &ranges, const double epsilon=1.0e-12)
Tests that the dg_package_data and dg_boundary_terms functions agree with the python implementation.
Definition: BoundaryCorrections.hpp:875
gr::spatial_metric
tnsr::ii< DataType, SpatialDim, Frame > spatial_metric(const tnsr::aa< DataType, SpatialDim, Frame > &spacetime_metric) noexcept
Compute spatial metric from spacetime metric.
optional
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
dg::Formulation
Formulation
The DG formulation to use.
Definition: Formulation.hpp:34
Prefixes.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:625
type_traits
TMPL.hpp
Mesh.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13
string