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  std::uniform_real_distribution<> dist(-1.0, 1.0);
210  DataVector used_for_size{face_mesh.number_of_grid_points()};
211 
212  // Fill all fields with random values in [-1,1), then, for each tag with a
213  // specified range, overwrite with new random values in [min,max)
214  Variables<dg_package_field_tags> interior_package_data{used_for_size.size()};
215  auto interior_fields_on_face =
216  make_with_random_values<Variables<face_tags_with_curved_background>>(
217  generator, make_not_null(&dist), used_for_size);
218  tmpl::for_each<tmpl::list<RangeTags...>>([&generator,
219  &interior_fields_on_face,
220  &ranges](auto tag_v) {
221  using tag = tmpl::type_from<decltype(tag_v)>;
222  const std::array<double, 2>& range = tuples::get<Tags::Range<tag>>(ranges);
223  std::uniform_real_distribution<> local_dist(range[0], range[1]);
224  fill_with_random_values(make_not_null(&get<tag>(interior_fields_on_face)),
225  generator, make_not_null(&local_dist));
226  });
227 
228  // Same as above but now for external data
229  Variables<dg_package_field_tags> exterior_package_data{used_for_size.size()};
230  auto exterior_fields_on_face =
231  make_with_random_values<Variables<face_tags_with_curved_background>>(
232  generator, make_not_null(&dist), used_for_size);
233  tmpl::for_each<tmpl::list<RangeTags...>>([&generator,
234  &exterior_fields_on_face,
235  &ranges](auto tag_v) {
236  using tag = tmpl::type_from<decltype(tag_v)>;
237  const std::array<double, 2>& range = tuples::get<Tags::Range<tag>>(ranges);
238  std::uniform_real_distribution<> local_dist(range[0], range[1]);
239  fill_with_random_values(make_not_null(&get<tag>(exterior_fields_on_face)),
240  generator, make_not_null(&local_dist));
241  });
242 
243  // Compute the interior and exterior normal vectors so they are pointing in
244  // opposite directions.
245  auto interior_unit_normal_covector = make_with_random_values<
246  tnsr::i<DataVector, FaceDim + 1, Frame::Inertial>>(
247  generator, make_not_null(&dist), used_for_size);
248  tnsr::I<DataVector, FaceDim + 1, Frame::Inertial>
249  interior_unit_normal_vector{};
250 
251  tnsr::i<DataVector, FaceDim + 1, Frame::Inertial>
252  exterior_unit_normal_covector;
253  tnsr::I<DataVector, FaceDim + 1, Frame::Inertial>
254  exterior_unit_normal_vector{};
255  if constexpr (not curved_background) {
256  const Scalar<DataVector> interior_normal_magnitude =
257  magnitude(interior_unit_normal_covector);
258  for (auto& t : interior_unit_normal_covector) {
259  t /= get(interior_normal_magnitude);
260  }
261  exterior_unit_normal_covector = interior_unit_normal_covector;
262  for (auto& t : exterior_unit_normal_covector) {
263  t *= -1.0;
264  }
265  } else {
266  using inv_spatial_metric = typename detail::inverse_spatial_metric_tag<
267  curved_background>::template f<System>;
268  exterior_unit_normal_covector = interior_unit_normal_covector;
269  for (auto& t : exterior_unit_normal_covector) {
270  t *= -1.0;
271  }
272  detail::adjust_spatial_metric_or_inverse(
273  make_not_null(&get<inv_spatial_metric>(interior_fields_on_face)));
274 
275  // make the exterior inverse spatial metric be close to the interior one
276  // since the solution should be smooth. We can't change too much because we
277  // want the spatial metric to have diagonal terms much larger than the
278  // off-diagonal terms.
279  std::uniform_real_distribution<> inv_metric_change_dist(0.999, 1.0);
280  for (size_t i = 0; i < FaceDim + 1; ++i) {
281  for (size_t j = i; j < FaceDim + 1; ++j) {
282  get<inv_spatial_metric>(exterior_fields_on_face).get(i, j) =
283  inv_metric_change_dist(*generator) *
284  get<inv_spatial_metric>(interior_fields_on_face).get(i, j);
285  }
286  }
287  detail::normalize_vector_and_covector(
288  make_not_null(&interior_unit_normal_covector),
289  make_not_null(&interior_unit_normal_vector),
290  get<inv_spatial_metric>(interior_fields_on_face));
291  detail::normalize_vector_and_covector(
292  make_not_null(&exterior_unit_normal_covector),
293  make_not_null(&exterior_unit_normal_vector),
294  get<inv_spatial_metric>(exterior_fields_on_face));
295 
296  // if dg_package_data depends on the (not inverted) spatial metric, we also
297  // adjust the spatial metric. Note that for testing purposes, the metric and
298  // inverse metric don't need to be inverses of each other. However, we do
299  // want both to be physically reasonable, i.e., close to flat space.
300  using spatial_metric =
302  if constexpr (tmpl::list_contains_v<face_tags_with_curved_background,
303  spatial_metric>) {
304  detail::adjust_spatial_metric_or_inverse(
305  make_not_null(&get<spatial_metric>(interior_fields_on_face)));
306  // as for external inverse metric: we get the external metric by slightly
307  // tweaking the internal metric
308  for (size_t i = 0; i < FaceDim + 1; ++i) {
309  for (size_t j = i; j < FaceDim + 1; ++j) {
310  get<spatial_metric>(exterior_fields_on_face).get(i, j) =
311  inv_metric_change_dist(*generator) *
312  get<spatial_metric>(interior_fields_on_face).get(i, j);
313  }
314  }
315  }
316  }
317 
319  mesh_velocity{};
320  if (use_moving_mesh) {
321  mesh_velocity = make_with_random_values<
322  tnsr::I<DataVector, FaceDim + 1, Frame::Inertial>>(
323  generator, make_not_null(&dist), used_for_size);
324  }
325 
326  if constexpr (curved_background) {
327  call_dg_package_data(
328  make_not_null(&interior_package_data), correction,
329  interior_fields_on_face, volume_data, interior_unit_normal_covector,
330  interior_unit_normal_vector, mesh_velocity, face_tags{});
331  call_dg_package_data(
332  make_not_null(&exterior_package_data), correction,
333  exterior_fields_on_face, volume_data, exterior_unit_normal_covector,
334  exterior_unit_normal_vector, mesh_velocity, face_tags{});
335  } else {
336  call_dg_package_data(make_not_null(&interior_package_data), correction,
337  interior_fields_on_face, volume_data,
338  interior_unit_normal_covector, mesh_velocity,
339  face_tags{});
340  call_dg_package_data(make_not_null(&exterior_package_data), correction,
341  exterior_fields_on_face, volume_data,
342  exterior_unit_normal_covector, mesh_velocity,
343  face_tags{});
344  }
345 
346  Variables<dt_variables_tags> boundary_corrections{
347  face_mesh.number_of_grid_points()};
348  call_dg_boundary_terms(make_not_null(&boundary_corrections), correction,
349  interior_package_data, exterior_package_data,
350  dg_formulation);
351 
352  if (dg_formulation == ::dg::Formulation::StrongInertial) {
353  // The strong form should be (WeakForm - (n_i F^i)_{interior}).
354  // Since we also test conservation for the weak form we just need to test
355  // that the strong form satisfies the above definition.
356 
357  if constexpr (curved_background) {
358  call_dg_package_data(
359  make_not_null(&interior_package_data), correction,
360  interior_fields_on_face, volume_data, interior_unit_normal_covector,
361  interior_unit_normal_vector, mesh_velocity, face_tags{});
362  call_dg_package_data(
363  make_not_null(&exterior_package_data), correction,
364  exterior_fields_on_face, volume_data, exterior_unit_normal_covector,
365  exterior_unit_normal_vector, mesh_velocity, face_tags{});
366  } else {
367  call_dg_package_data(make_not_null(&interior_package_data), correction,
368  interior_fields_on_face, volume_data,
369  interior_unit_normal_covector, mesh_velocity,
370  face_tags{});
371  call_dg_package_data(make_not_null(&exterior_package_data), correction,
372  exterior_fields_on_face, volume_data,
373  exterior_unit_normal_covector, mesh_velocity,
374  face_tags{});
375  }
376 
377  Variables<dt_variables_tags> expected_boundary_corrections{
378  face_mesh.number_of_grid_points()};
379  call_dg_boundary_terms(make_not_null(&expected_boundary_corrections),
380  correction, interior_package_data,
381  exterior_package_data,
382  ::dg::Formulation::WeakInertial);
383 
384  tmpl::for_each<flux_variables>([&interior_package_data,
385  &expected_boundary_corrections](
386  auto flux_variable_tag_v) noexcept {
387  using flux_variable_tag = tmpl::type_from<decltype(flux_variable_tag_v)>;
388  using normal_dot_flux_tag = ::Tags::NormalDotFlux<flux_variable_tag>;
389  using dt_tag = ::Tags::dt<flux_variable_tag>;
390  const auto& normal_dot_flux =
391  get<normal_dot_flux_tag>(interior_package_data);
392  auto& expected_boundary_correction =
393  get<dt_tag>(expected_boundary_corrections);
394  for (size_t tensor_index = 0;
395  tensor_index < expected_boundary_correction.size(); ++tensor_index) {
396  expected_boundary_correction[tensor_index] -=
397  normal_dot_flux[tensor_index];
398  }
399  });
400  {
401  INFO("Check weak and strong boundary terms match.");
403  boundary_corrections, expected_boundary_corrections, custom_approx);
404  }
405 
406  if (zero_on_smooth_solution == ZeroOnSmoothSolution::Yes) {
407  INFO(
408  "Testing that if the solution is the same on both sides the "
409  "StrongInertial correction is identically zero.");
410  Variables<dg_package_field_tags> interior_package_data_opposite_signs{
411  used_for_size.size()};
412  if constexpr (curved_background) {
413  // If the solution is the same on both sides then the interior and
414  // exterior normal vectors only differ by a sign.
415  for (size_t i = 0; i < FaceDim + 1; ++i) {
416  exterior_unit_normal_vector.get(i) =
417  -interior_unit_normal_vector.get(i);
418  exterior_unit_normal_covector.get(i) =
419  -interior_unit_normal_covector.get(i);
420  }
421  call_dg_package_data(
422  make_not_null(&interior_package_data_opposite_signs), correction,
423  interior_fields_on_face, volume_data, exterior_unit_normal_covector,
424  exterior_unit_normal_vector, mesh_velocity, face_tags{});
425  } else {
426  call_dg_package_data(
427  make_not_null(&interior_package_data_opposite_signs), correction,
428  interior_fields_on_face, volume_data, exterior_unit_normal_covector,
429  mesh_velocity, face_tags{});
430  }
431  Variables<dt_variables_tags> zero_boundary_correction{
432  face_mesh.number_of_grid_points()};
433  call_dg_boundary_terms(make_not_null(&zero_boundary_correction),
434  correction, interior_package_data,
435  interior_package_data_opposite_signs,
436  ::dg::Formulation::StrongInertial);
437  Variables<dt_variables_tags> expected_zero_boundary_correction{
438  face_mesh.number_of_grid_points(), 0.0};
439  tmpl::for_each<dt_variables_tags>([&custom_approx,
440  &expected_zero_boundary_correction,
441  &zero_boundary_correction](
442  auto dt_variables_tag_v) noexcept {
443  using dt_variables_tag = tmpl::type_from<decltype(dt_variables_tag_v)>;
444  const std::string tag_name = db::tag_name<dt_variables_tag>();
445  CAPTURE(tag_name);
447  get<dt_variables_tag>(zero_boundary_correction),
448  get<dt_variables_tag>(expected_zero_boundary_correction),
449  custom_approx);
450  });
451  }
452  } else if (dg_formulation == ::dg::Formulation::WeakInertial) {
453  INFO(
454  "Checking that swapping the two sides results in an overall minus "
455  "sign.");
456  Variables<dt_variables_tags> reverse_side_boundary_corrections{
457  face_mesh.number_of_grid_points()};
458  call_dg_boundary_terms(make_not_null(&reverse_side_boundary_corrections),
459  correction, exterior_package_data,
460  interior_package_data, dg_formulation);
461  // Check that the flux leaving one element equals the flux entering its
462  // neighbor, i.e., F*(interior, exterior) == -F*(exterior, interior)
463  reverse_side_boundary_corrections *= -1.0;
464  tmpl::for_each<flux_variables>([&boundary_corrections, &custom_approx,
465  &reverse_side_boundary_corrections](
466  auto flux_variable_tag_v) {
467  using flux_variable_tag = tmpl::type_from<decltype(flux_variable_tag_v)>;
468  const std::string tag_name = db::tag_name<flux_variable_tag>();
469  CAPTURE(tag_name);
471  get<::Tags::dt<flux_variable_tag>>(reverse_side_boundary_corrections),
472  get<::Tags::dt<flux_variable_tag>>(boundary_corrections),
473  custom_approx);
474  });
475  } else {
476  ERROR("DG formulation must be StrongInertial or WeakInertial, not "
477  << dg_formulation);
478  }
479 }
480 } // namespace detail
481 
482 /*!
483  * \ingroup TestingFrameworkGroup
484  * \brief Checks that the boundary correction is conservative and that for
485  * smooth solutions the strong-form correction is zero.
486  */
487 template <typename System, typename BoundaryCorrection, size_t FaceDim,
488  typename... VolumeTags, typename... RangeTags>
490  const gsl::not_null<std::mt19937*> generator,
491  const BoundaryCorrection& correction, const Mesh<FaceDim>& face_mesh,
492  const tuples::TaggedTuple<VolumeTags...>& volume_data,
494  const ZeroOnSmoothSolution zero_on_smooth_solution =
495  ZeroOnSmoothSolution::Yes,
496  const double eps = 1.0e-12) {
497  for (const auto use_moving_mesh : {true, false}) {
498  for (const auto& dg_formulation :
499  {::dg::Formulation::StrongInertial, ::dg::Formulation::WeakInertial}) {
500  detail::test_boundary_correction_conservation_impl<System>(
501  generator, correction, face_mesh, volume_data, ranges,
502  use_moving_mesh, dg_formulation, zero_on_smooth_solution, eps);
503  }
504  }
505 }
506 
507 namespace detail {
508 template <typename System, typename ConversionClassList, typename VariablesTags,
509  typename BoundaryCorrection, size_t FaceDim, typename... FaceTags,
510  typename... VolumeTags, typename... RangeTags,
511  typename... DgPackageDataTags>
512 void test_with_python(
513  const gsl::not_null<std::mt19937*> generator,
514  const std::string& python_module,
515  const std::array<
516  std::string,
517  tmpl::size<typename BoundaryCorrection::dg_package_field_tags>::value>&
518  python_dg_package_data_functions,
519  const std::array<std::string, tmpl::size<VariablesTags>::value>&
520  python_dg_boundary_terms_functions,
521  const BoundaryCorrection& correction, const Mesh<FaceDim>& face_mesh,
522  const tuples::TaggedTuple<VolumeTags...>& volume_data,
523  const tuples::TaggedTuple<Tags::Range<RangeTags>...>& ranges,
524  const bool use_moving_mesh, const ::dg::Formulation dg_formulation,
525  const double epsilon, tmpl::list<FaceTags...> /*meta*/,
526  tmpl::list<DgPackageDataTags...> /*meta*/) {
527  CAPTURE(face_mesh);
528  CAPTURE(dg_formulation);
529  CAPTURE(use_moving_mesh);
530  REQUIRE(face_mesh.number_of_grid_points() >= 1);
531  constexpr bool curved_background =
532  detail::has_inverse_spatial_metric_tag_v<System>;
533  using dg_package_field_tags =
534  typename BoundaryCorrection::dg_package_field_tags;
535 
536  using face_tags = tmpl::list<FaceTags...>;
537  using face_tags_with_curved_background = tmpl::conditional_t<
538  curved_background,
539  tmpl::remove_duplicates<tmpl::push_back<
540  face_tags, typename TestHelpers::evolution::dg::detail::
541  inverse_spatial_metric_tag<
542  curved_background>::template f<System>>>,
543  face_tags>;
544 
545  // Sanity check: we apply the same ranges to the random inputs of the
546  // `dg_package_data` function and the `dg_boundary_terms` function. So
547  // first we check that each range tag appears in at least one of these
548  // function's arguments.
549  using range_tags_list = tmpl::list<RangeTags...>;
550  using ranges_for_dg_boundary_terms_only =
552  using ranges_unused = tmpl::list_difference<ranges_for_dg_boundary_terms_only,
553  dg_package_field_tags>;
554  static_assert(std::is_same_v<tmpl::list<>, ranges_unused>,
555  "Received Tags::Range for Tags that are neither arguments to "
556  "dg_package_data nor dg_boundary_terms");
557 
558  std::uniform_real_distribution<> dist(-1.0, 1.0);
559  DataVector used_for_size{face_mesh.number_of_grid_points()};
560 
561  // Fill all fields with random values in [-1,1), then, for each tag with a
562  // specified range, overwrite with new random values in [min,max)
563  auto fields_on_face =
564  make_with_random_values<Variables<face_tags_with_curved_background>>(
565  generator, make_not_null(&dist), used_for_size);
566  tmpl::for_each<tmpl::list<RangeTags...>>([&generator, &fields_on_face,
567  &ranges](auto tag_v) {
568  using tag = tmpl::type_from<decltype(tag_v)>;
569  const std::array<double, 2>& range = tuples::get<Tags::Range<tag>>(ranges);
570  std::uniform_real_distribution<> local_dist(range[0], range[1]);
571  fill_with_random_values(make_not_null(&get<tag>(fields_on_face)), generator,
572  make_not_null(&local_dist));
573  });
574 
575  auto unit_normal_covector = make_with_random_values<
576  tnsr::i<DataVector, FaceDim + 1, Frame::Inertial>>(
577  generator, make_not_null(&dist), used_for_size);
578  tnsr::I<DataVector, FaceDim + 1, Frame::Inertial> unit_normal_vector{};
579 
580  if constexpr (not curved_background) {
581  const Scalar<DataVector> normal_magnitude = magnitude(unit_normal_covector);
582  for (auto& t : unit_normal_covector) {
583  t /= get(normal_magnitude);
584  }
585  } else {
586  using inv_spatial_metric = typename detail::inverse_spatial_metric_tag<
587  curved_background>::template f<System>;
588  detail::adjust_spatial_metric_or_inverse(
589  make_not_null(&get<inv_spatial_metric>(fields_on_face)));
590  detail::normalize_vector_and_covector(
591  make_not_null(&unit_normal_covector),
592  make_not_null(&unit_normal_vector),
593  get<inv_spatial_metric>(fields_on_face));
594 
595  using spatial_metric =
597  if constexpr (tmpl::list_contains_v<face_tags_with_curved_background,
598  spatial_metric>) {
599  detail::adjust_spatial_metric_or_inverse(
600  make_not_null(&get<spatial_metric>(fields_on_face)));
601  }
602  }
603 
605  mesh_velocity{};
606  if (use_moving_mesh) {
607  mesh_velocity = make_with_random_values<
608  tnsr::I<DataVector, FaceDim + 1, Frame::Inertial>>(
609  generator, make_not_null(&dist), used_for_size);
610  }
611  std::optional<Scalar<DataVector>> normal_dot_mesh_velocity{};
612  if (mesh_velocity.has_value()) {
613  normal_dot_mesh_velocity =
614  dot_product(*mesh_velocity, unit_normal_covector);
615  }
616 
617  // Call C++ implementation of dg_package_data
618  Variables<dg_package_field_tags> package_data{used_for_size.size()};
619  if constexpr (curved_background) {
620  call_dg_package_data(make_not_null(&package_data), correction,
621  fields_on_face, volume_data, unit_normal_covector,
622  unit_normal_vector, mesh_velocity,
623  tmpl::list<FaceTags...>{});
624  } else {
625  call_dg_package_data(make_not_null(&package_data), correction,
626  fields_on_face, volume_data, unit_normal_covector,
627  mesh_velocity, tmpl::list<FaceTags...>{});
628  }
629 
630  // Call python implementation of dg_package_data
631  size_t function_name_index = 0;
632  tmpl::for_each<dg_package_field_tags>(
633  [epsilon, &fields_on_face, &function_name_index, &mesh_velocity,
634  &normal_dot_mesh_velocity, &package_data,
635  &python_dg_package_data_functions, &python_module, &unit_normal_covector,
636  &unit_normal_vector, &volume_data](auto package_data_tag_v) {
637  // avoid compiler warnings if there isn't any volume data
638  (void)volume_data;
639  INFO("Testing package data");
640  using package_data_tag = tmpl::type_from<decltype(package_data_tag_v)>;
641  using ResultType = typename package_data_tag::type;
642  const std::string tag_name = db::tag_name<package_data_tag>();
643  CAPTURE(tag_name);
644  try {
645  CAPTURE(python_module);
646  CAPTURE(
647  gsl::at(python_dg_package_data_functions, function_name_index));
648  CAPTURE(pretty_type::short_name<ResultType>());
649  if constexpr (curved_background) {
650  const auto python_result =
651  pypp::call<ResultType, ConversionClassList>(
652  python_module,
653  gsl::at(python_dg_package_data_functions,
654  function_name_index),
655  get<FaceTags>(fields_on_face)..., unit_normal_covector,
656  unit_normal_vector, mesh_velocity, normal_dot_mesh_velocity,
657  get<VolumeTags>(volume_data)...);
659  get<package_data_tag>(package_data), python_result,
660  Approx::custom().epsilon(epsilon).scale(1.0));
661  } else {
662  (void)unit_normal_vector;
663  const auto python_result =
664  pypp::call<ResultType, ConversionClassList>(
665  python_module,
666  gsl::at(python_dg_package_data_functions,
667  function_name_index),
668  get<FaceTags>(fields_on_face)..., unit_normal_covector,
669  mesh_velocity, normal_dot_mesh_velocity,
670  get<VolumeTags>(volume_data)...);
672  get<package_data_tag>(package_data), python_result,
673  Approx::custom().epsilon(epsilon).scale(1.0));
674  }
675  } catch (const std::exception& e) {
676  INFO("On line " << __LINE__ << " Python call to "
677  << gsl::at(python_dg_package_data_functions,
678  function_name_index)
679  << " in module " << python_module
680  << " failed: " << e.what());
681  REQUIRE(false);
682  }
683  ++function_name_index;
684  });
685 
686  // Now we need to check the dg_boundary_terms function.
687  //
688  // In order for the package data to have a self-consistent state (not all
689  // things being packaged are independent), we compute the packaged data.
690  // We reuse the previously computed packaged data as the interior packaged
691  // data.
692  const auto& interior_package_data = package_data;
693 
694  auto exterior_fields_on_face =
695  make_with_random_values<Variables<face_tags_with_curved_background>>(
696  generator, make_not_null(&dist), used_for_size);
697  tmpl::for_each<tmpl::list<RangeTags...>>([&generator,
698  &exterior_fields_on_face,
699  &ranges](auto tag_v) {
700  using tag = tmpl::type_from<decltype(tag_v)>;
701  const std::array<double, 2>& range = tuples::get<Tags::Range<tag>>(ranges);
702  std::uniform_real_distribution<> local_dist(range[0], range[1]);
703  fill_with_random_values(make_not_null(&get<tag>(exterior_fields_on_face)),
704  generator, make_not_null(&local_dist));
705  });
706 
707  tnsr::i<DataVector, FaceDim + 1, Frame::Inertial>
708  exterior_unit_normal_covector;
709  tnsr::I<DataVector, FaceDim + 1, Frame::Inertial>
710  exterior_unit_normal_vector{};
711  if constexpr (not curved_background) {
712  exterior_unit_normal_covector = unit_normal_covector;
713  for (auto& t : exterior_unit_normal_covector) {
714  t *= -1.0;
715  }
716  } else {
717  using inv_spatial_metric = typename detail::inverse_spatial_metric_tag<
718  curved_background>::template f<System>;
719  exterior_unit_normal_covector = unit_normal_covector;
720  for (auto& t : exterior_unit_normal_covector) {
721  t *= -1.0;
722  }
723 
724  // make the exterior inverse spatial metric be close to the interior one
725  // since the solution should be smooth. We can't change too much because we
726  // want the spatial metric to have diagonal terms much larger than the
727  // off-diagonal terms.
728  std::uniform_real_distribution<> inv_metric_change_dist(0.999, 1.0);
729  for (size_t i = 0; i < FaceDim + 1; ++i) {
730  for (size_t j = i; j < FaceDim + 1; ++j) {
731  get<inv_spatial_metric>(exterior_fields_on_face).get(i, j) =
732  inv_metric_change_dist(*generator) *
733  get<inv_spatial_metric>(fields_on_face).get(i, j);
734  }
735  }
736  detail::normalize_vector_and_covector(
737  make_not_null(&exterior_unit_normal_covector),
738  make_not_null(&exterior_unit_normal_vector),
739  get<inv_spatial_metric>(exterior_fields_on_face));
740 
741  // if dg_package_data depends on the (not inverted) spatial metric, we also
742  // adjust the spatial metric. Note that for testing purposes, the metric and
743  // inverse metric don't need to be inverses of each other. However, we do
744  // want both to be physically reasonable, i.e., close to flat space.
745  using spatial_metric =
747  if constexpr (tmpl::list_contains_v<face_tags_with_curved_background,
748  spatial_metric>) {
749  // as for external inverse metric: we get the external metric by slightly
750  // tweaking the internal metric
751  for (size_t i = 0; i < FaceDim + 1; ++i) {
752  for (size_t j = i; j < FaceDim + 1; ++j) {
753  get<spatial_metric>(exterior_fields_on_face).get(i, j) =
754  inv_metric_change_dist(*generator) *
755  get<spatial_metric>(fields_on_face).get(i, j);
756  }
757  }
758  }
759  }
760 
761  Variables<dg_package_field_tags> exterior_package_data{used_for_size.size()};
762  if constexpr (curved_background) {
763  call_dg_package_data(
764  make_not_null(&exterior_package_data), correction,
765  exterior_fields_on_face, volume_data, exterior_unit_normal_covector,
766  exterior_unit_normal_vector, mesh_velocity, face_tags{});
767  } else {
768  call_dg_package_data(make_not_null(&exterior_package_data), correction,
769  exterior_fields_on_face, volume_data,
770  exterior_unit_normal_covector, mesh_velocity,
771  face_tags{});
772  }
773 
774  // We don't need to prefix the VariablesTags with anything because we are not
775  // interacting with any code that cares about what the tags are, just that the
776  // types matched the evolved variables.
777  Variables<VariablesTags> boundary_corrections{
778  face_mesh.number_of_grid_points()};
779 
780  // Call C++ implementation of dg_boundary_terms
781  call_dg_boundary_terms(make_not_null(&boundary_corrections), correction,
782  interior_package_data, exterior_package_data,
783  dg_formulation);
784 
785  // Call python implementation of dg_boundary_terms
786  function_name_index = 0;
787  tmpl::for_each<VariablesTags>([&boundary_corrections, &dg_formulation,
788  epsilon, &exterior_package_data,
789  &function_name_index, &interior_package_data,
790  &python_dg_boundary_terms_functions,
791  &python_module](auto package_data_tag_v) {
792  INFO("Testing boundary terms");
793  using boundary_correction_tag =
794  tmpl::type_from<decltype(package_data_tag_v)>;
795  using ResultType = typename boundary_correction_tag::type;
796  try {
797  CAPTURE(python_module);
798  const std::string& python_function =
799  gsl::at(python_dg_boundary_terms_functions, function_name_index);
800  CAPTURE(python_function);
801  // Make explicitly depend on tag type to avoid type deduction issues with
802  // GCC7
803  const typename boundary_correction_tag::type python_result =
804  pypp::call<ResultType>(
805  python_module, python_function,
806  get<DgPackageDataTags>(interior_package_data)...,
807  get<DgPackageDataTags>(exterior_package_data)...,
808  dg_formulation == ::dg::Formulation::StrongInertial);
810  get<boundary_correction_tag>(boundary_corrections), python_result,
811  Approx::custom().epsilon(epsilon).scale(1.0));
812  } catch (const std::exception& e) {
813  INFO("On line " << __LINE__ << " Python call to "
814  << gsl::at(python_dg_boundary_terms_functions,
815  function_name_index)
816  << " in module " << python_module
817  << " failed: " << e.what());
818  REQUIRE(false);
819  }
820  ++function_name_index;
821  });
822 }
823 } // namespace detail
824 
825 /*!
826  * \ingroup TestingFrameworkGroup
827  * \brief Tests that the `dg_package_data` and `dg_boundary_terms` functions
828  * agree with the python implementation.
829  *
830  * The variables are filled with random numbers between zero and one before
831  * being passed to the implementations. If in the future we need support for
832  * negative numbers we can add the ability to specify a single range for all
833  * random numbers or each individually.
834  *
835  * Please note the following:
836  * - The `pypp::SetupLocalPythonEnvironment` must be created before the
837  * `test_boundary_correction_with_python` can be called.
838  * - The `dg_formulation` is passed as a bool `use_strong_form` to the python
839  * functions since we don't want to rely on python bindings for the enum.
840  * - You can convert custom types using the `ConversionClassList` template
841  * parameter, which is then passed to `pypp::call()`. This allows you to,
842  * e.g., convert an equation of state into an array locally in a test file.
843  * - There must be one python function to compute the packaged data for each tag
844  * in `dg_package_field_tags`
845  * - There must be one python function to compute the boundary correction for
846  * each tag in `System::variables_tag`
847  * - The arguments to the python functions for computing the packaged data are
848  * the same as the arguments for the C++ `dg_package_data` function, excluding
849  * the `gsl::not_null` arguments.
850  * - The arguments to the python functions for computing the boundary
851  * corrections are the same as the arguments for the C++ `dg_boundary_terms`
852  * function, excluding the `gsl::not_null` arguments.
853  * - `ranges` is a `TaggedTuple` of
854  * `TestHelpers::evolution::dg::Tags::Range<tag>` specifying a custom range in
855  * which to generate the random values. This can be used for ensuring that
856  * positive quantities are randomly generated on the interval
857  * `[lower_bound,upper_bound)`, choosing `lower_bound` to be `0` or some small
858  * number. The default interval if a tag is not listed is `[-1,1)`. The range
859  * is used for setting the random inputs to `dg_package_data`.
860  */
861 template <typename System, typename ConversionClassList = tmpl::list<>,
862  typename BoundaryCorrection, size_t FaceDim, typename... VolumeTags,
863  typename... RangeTags>
865  const gsl::not_null<std::mt19937*> generator,
866  const std::string& python_module,
867  const std::array<
868  std::string,
869  tmpl::size<typename BoundaryCorrection::dg_package_field_tags>::value>&
870  python_dg_package_data_functions,
871  const std::array<
872  std::string,
873  tmpl::size<typename System::variables_tag::tags_list>::value>&
874  python_dg_boundary_terms_functions,
875  const BoundaryCorrection& correction, const Mesh<FaceDim>& face_mesh,
876  const tuples::TaggedTuple<VolumeTags...>& volume_data,
878  const double epsilon = 1.0e-12) {
879  static_assert(std::is_final_v<std::decay_t<BoundaryCorrection>>,
880  "All boundary correction classes must be marked `final`.");
881  using package_temporary_tags =
882  typename BoundaryCorrection::dg_package_data_temporary_tags;
883  using package_primitive_tags = detail::get_correction_primitive_vars<
884  System::has_primitive_and_conservative_vars, BoundaryCorrection>;
885  using variables_tags = typename System::variables_tag::tags_list;
886  using flux_variables = typename System::flux_variables;
887  using flux_tags =
890 
891  for (const auto use_moving_mesh : {false, true}) {
892  for (const auto dg_formulation :
893  {::dg::Formulation::StrongInertial, ::dg::Formulation::WeakInertial}) {
894  detail::test_with_python<System, ConversionClassList, variables_tags>(
895  generator, python_module, python_dg_package_data_functions,
896  python_dg_boundary_terms_functions, correction, face_mesh,
897  volume_data, ranges, use_moving_mesh, dg_formulation, epsilon,
898  tmpl::append<variables_tags, flux_tags, package_temporary_tags,
899  package_primitive_tags>{},
900  typename BoundaryCorrection::dg_package_field_tags{});
901  }
902  }
903 }
904 } // 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:489
std::string
Frame::Inertial
Definition: IndexType.hpp:44
std::exception
gr::Tags::SpatialMetric
Definition: Tags.hpp:26
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:660
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
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:156
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:48
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:864
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:589
type_traits
TMPL.hpp
Mesh.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecThirdOrderPiecewisePolynomial.hpp:13
string