14 #include "DataStructures/DataBox/PrefixHelpers.hpp"
16 #include "DataStructures/DataVector.hpp"
18 #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
24 #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
28 #include "Utilities/TaggedTuple.hpp"
30 namespace TestHelpers::evolution::dg {
33 enum class ZeroOnSmoothSolution { Yes, No };
36 template <
bool HasPrimitiveVars = false>
37 struct get_correction_primitive_vars_impl {
38 template <
typename BoundaryCorrection>
39 using f = tmpl::list<>;
43 struct get_correction_primitive_vars_impl<true> {
44 template <
typename BoundaryCorrection>
45 using f =
typename BoundaryCorrection::dg_package_data_primitive_tags;
48 template <
bool HasPrimitiveVars,
typename BoundaryCorrection>
49 using get_correction_primitive_vars =
50 typename get_correction_primitive_vars_impl<HasPrimitiveVars>::template f<
53 template <
bool HasPrimitiveVars = false>
54 struct get_system_primitive_vars_impl {
55 template <
typename System>
56 using f = tmpl::list<>;
60 struct get_system_primitive_vars_impl<true> {
61 template <
typename System>
62 using f =
typename System::primitive_variables_tag::tags_list;
65 template <
bool HasPrimitiveVars,
typename System>
66 using get_system_primitive_vars =
typename get_system_primitive_vars_impl<
67 HasPrimitiveVars>::template f<System>;
69 template <
typename BoundaryCorrection,
typename... PackageTags,
70 typename... FaceTags,
typename... VolumeTags,
size_t Dim>
71 void call_dg_package_data(
72 const gsl::not_null<Variables<tmpl::list<PackageTags...>>*> package_data,
73 const BoundaryCorrection& correction,
74 const Variables<tmpl::list<FaceTags...>>& face_variables,
76 const tnsr::i<DataVector, Dim, Frame::Inertial>& unit_normal_covector,
77 const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
80 if (mesh_velocity.has_value()) {
81 normal_dot_mesh_velocity =
84 correction.dg_package_data(
86 get<FaceTags>(face_variables)..., unit_normal_covector, mesh_velocity,
87 normal_dot_mesh_velocity, get<VolumeTags>(volume_data)...);
90 template <
typename BoundaryCorrection,
typename... BoundaryCorrectionTags,
91 typename... PackageTags>
92 void call_dg_boundary_terms(
93 const gsl::not_null<Variables<tmpl::list<BoundaryCorrectionTags...>>*>
95 const BoundaryCorrection& correction,
96 const Variables<tmpl::list<PackageTags...>>& interior_package_data,
97 const Variables<tmpl::list<PackageTags...>>& exterior_package_data,
99 correction.dg_boundary_terms(
100 make_not_null(&get<BoundaryCorrectionTags>(*boundary_corrections))...,
101 get<PackageTags>(interior_package_data)...,
102 get<PackageTags>(exterior_package_data)..., dg_formulation);
105 template <
typename System,
typename BoundaryCorrection,
size_t FaceDim,
106 typename... VolumeTags>
107 void test_boundary_correction_impl(
108 const BoundaryCorrection& correction,
const Mesh<FaceDim>& face_mesh,
111 const ZeroOnSmoothSolution zero_on_smooth_solution) {
112 CAPTURE(use_moving_mesh);
113 CAPTURE(dg_formulation);
115 using variables_tags =
typename System::variables_tag::tags_list;
116 using primitive_variables_tags = detail::get_system_primitive_vars<
117 System::has_primitive_and_conservative_vars, System>;
118 using flux_variables =
typename System::flux_variables;
122 using temporary_tags =
123 typename System::compute_volume_time_derivative_terms::temporary_tags;
126 using dg_package_field_tags =
127 typename BoundaryCorrection::dg_package_field_tags;
128 using package_temporary_tags =
129 typename BoundaryCorrection::dg_package_data_temporary_tags;
130 using package_primitive_tags = detail::get_correction_primitive_vars<
131 System::has_primitive_and_conservative_vars, BoundaryCorrection>;
138 tmpl::list_difference<package_temporary_tags, temporary_tags>,
140 "There are temporary tags needed by the boundary correction that are not "
141 "computed as temporary tags by the system");
147 std::is_same_v<tmpl::list_difference<package_primitive_tags,
148 primitive_variables_tags>,
150 "There are primitive tags needed by the boundary correction that are not "
151 "listed in the system as being primitive variables");
160 if (use_moving_mesh) {
162 tnsr::I<DataVector, FaceDim + 1, Frame::Inertial>>(
166 Variables<dg_package_field_tags> interior_package_data{used_for_size.size()};
168 Variables<tmpl::append<variables_tags, flux_tags, package_temporary_tags,
169 package_primitive_tags>>>(
172 Variables<dg_package_field_tags> exterior_package_data{used_for_size.size()};
174 Variables<tmpl::append<variables_tags, flux_tags, package_temporary_tags,
175 package_primitive_tags>>>(
181 tnsr::i<DataVector, FaceDim + 1, Frame::Inertial>>(
184 magnitude(interior_unit_normal_covector);
185 for (
auto& t : interior_unit_normal_covector) {
186 t /=
get(interior_normal_magnitude);
188 auto exterior_unit_normal_covector = interior_unit_normal_covector;
189 for (
auto& t : exterior_unit_normal_covector) {
193 call_dg_package_data(
make_not_null(&interior_package_data), correction,
194 interior_fields_on_face, volume_data,
195 interior_unit_normal_covector, mesh_velocity);
196 call_dg_package_data(
make_not_null(&exterior_package_data), correction,
197 exterior_fields_on_face, volume_data,
198 exterior_unit_normal_covector, mesh_velocity);
200 Variables<dt_variables_tags> boundary_corrections{
202 call_dg_boundary_terms(
make_not_null(&boundary_corrections), correction,
203 interior_package_data, exterior_package_data,
206 if (dg_formulation == ::dg::Formulation::StrongInertial) {
211 Variables<dt_variables_tags> expected_boundary_corrections{
213 call_dg_boundary_terms(
make_not_null(&expected_boundary_corrections),
214 correction, interior_package_data,
215 exterior_package_data,
216 ::dg::Formulation::WeakInertial);
218 tmpl::for_each<flux_variables>([&interior_package_data,
219 &expected_boundary_corrections](
220 auto flux_variable_tag_v) noexcept {
221 using flux_variable_tag = tmpl::type_from<decltype(flux_variable_tag_v)>;
224 const auto& normal_dot_flux =
225 get<normal_dot_flux_tag>(interior_package_data);
226 auto& expected_boundary_correction =
227 get<dt_tag>(expected_boundary_corrections);
228 for (
size_t tensor_index = 0;
229 tensor_index < expected_boundary_correction.size(); ++tensor_index) {
230 expected_boundary_correction[tensor_index] -=
231 normal_dot_flux[tensor_index];
235 INFO(
"Check weak and strong boundary terms match.");
236 CHECK(boundary_corrections == expected_boundary_corrections);
239 if (zero_on_smooth_solution == ZeroOnSmoothSolution::Yes) {
241 "Testing that if the solution is the same on both sides the "
242 "StrongInertial correction is identically zero.");
243 Variables<dg_package_field_tags> interior_package_data_opposite_signs{
244 used_for_size.size()};
245 call_dg_package_data(
make_not_null(&interior_package_data_opposite_signs),
246 correction, interior_fields_on_face, volume_data,
247 exterior_unit_normal_covector, mesh_velocity);
248 Variables<dt_variables_tags> zero_boundary_correction{
250 call_dg_boundary_terms(
make_not_null(&zero_boundary_correction),
251 correction, interior_package_data,
252 interior_package_data_opposite_signs,
253 ::dg::Formulation::StrongInertial);
254 Variables<dt_variables_tags> expected_zero_boundary_correction{
256 tmpl::for_each<dt_variables_tags>([&expected_zero_boundary_correction,
257 &zero_boundary_correction](
258 auto dt_variables_tag_v) noexcept {
259 using dt_variables_tag = tmpl::type_from<decltype(dt_variables_tag_v)>;
262 Approx custom_approx = Approx::custom().epsilon(1.e-12).scale(1.0);
264 get<dt_variables_tag>(zero_boundary_correction),
265 get<dt_variables_tag>(expected_zero_boundary_correction),
269 }
else if (dg_formulation == ::dg::Formulation::WeakInertial) {
271 "Checking that swapping the two sides results in an overall minus "
273 Variables<dt_variables_tags> reverse_side_boundary_corrections{
275 call_dg_boundary_terms(
make_not_null(&reverse_side_boundary_corrections),
276 correction, exterior_package_data,
277 interior_package_data, dg_formulation);
280 reverse_side_boundary_corrections *= -1.0;
281 tmpl::for_each<flux_variables>([&boundary_corrections,
282 &reverse_side_boundary_corrections](
283 auto flux_variable_tag_v) {
284 using flux_variable_tag = tmpl::type_from<decltype(flux_variable_tag_v)>;
292 ERROR(
"DG formulation must be StrongInertial or WeakInertial, not "
303 template <
typename System,
typename BoundaryCorrection,
size_t FaceDim,
304 typename... VolumeTags>
306 const BoundaryCorrection& correction,
const Mesh<FaceDim>& face_mesh,
308 const ZeroOnSmoothSolution zero_on_smooth_solution =
309 ZeroOnSmoothSolution::Yes) {
310 for (
const auto use_moving_mesh : {
true,
false}) {
311 for (
const auto& dg_formulation :
312 {::dg::Formulation::StrongInertial, ::dg::Formulation::WeakInertial}) {
313 detail::test_boundary_correction_impl<System>(
314 correction, face_mesh, volume_data, use_moving_mesh, dg_formulation,
315 zero_on_smooth_solution);
321 template <
typename ConversionClassList,
typename VariablesTags,
322 typename BoundaryCorrection,
size_t FaceDim,
typename... FaceTags,
323 typename... VolumeTags,
typename... DgPackageDataTags>
324 void test_with_python(
328 tmpl::size<typename BoundaryCorrection::dg_package_field_tags>::value>&
329 python_dg_package_data_functions,
331 python_dg_boundary_terms_functions,
332 const BoundaryCorrection& correction,
const Mesh<FaceDim>& face_mesh,
335 const double epsilon, tmpl::list<FaceTags...> ,
336 tmpl::list<DgPackageDataTags...> ) {
338 CAPTURE(dg_formulation);
339 CAPTURE(use_moving_mesh);
341 using dg_package_field_tags =
342 typename BoundaryCorrection::dg_package_field_tags;
348 Variables<dg_package_field_tags> package_data{used_for_size.size()};
349 const auto fields_on_face =
353 tnsr::i<DataVector, FaceDim + 1, Frame::Inertial>>(
355 const auto normal_magnitude =
magnitude(unit_normal_covector);
356 for (
auto& component : unit_normal_covector) {
357 component /=
get(normal_magnitude);
361 if (use_moving_mesh) {
363 tnsr::I<DataVector, FaceDim + 1, Frame::Inertial>>(
367 if (mesh_velocity.has_value()) {
368 normal_dot_mesh_velocity =
373 call_dg_package_data(
make_not_null(&package_data), correction, fields_on_face,
374 volume_data, unit_normal_covector, mesh_velocity);
377 size_t function_name_index = 0;
378 tmpl::for_each<dg_package_field_tags>(
379 [epsilon, &fields_on_face, &function_name_index, &mesh_velocity,
380 &normal_dot_mesh_velocity, &package_data,
381 &python_dg_package_data_functions, &python_module, &unit_normal_covector,
382 &volume_data](
auto package_data_tag_v) {
385 INFO(
"Testing package data");
386 using package_data_tag = tmpl::type_from<decltype(package_data_tag_v)>;
387 using ResultType =
typename package_data_tag::type;
389 CAPTURE(python_module);
391 gsl::at(python_dg_package_data_functions, function_name_index));
392 const auto python_result =
393 pypp::call<ResultType, ConversionClassList>(
395 gsl::at(python_dg_package_data_functions,
396 function_name_index),
397 get<FaceTags>(fields_on_face)..., unit_normal_covector,
398 mesh_velocity, normal_dot_mesh_velocity,
399 get<VolumeTags>(volume_data)...);
401 get<package_data_tag>(package_data), python_result,
402 Approx::custom().epsilon(epsilon).scale(1.0));
404 INFO(
"On line " << __LINE__ <<
" Python call to "
405 <<
gsl::at(python_dg_package_data_functions,
407 <<
" in module " << python_module
408 <<
" failed: " << e.what());
411 ++function_name_index;
415 const auto interior_package_data =
416 make_with_random_values<Variables<dg_package_field_tags>>(
418 const auto exterior_package_data =
419 make_with_random_values<Variables<dg_package_field_tags>>(
424 Variables<VariablesTags> boundary_corrections{
428 call_dg_boundary_terms(
make_not_null(&boundary_corrections), correction,
429 interior_package_data, exterior_package_data,
433 function_name_index = 0;
434 tmpl::for_each<VariablesTags>([&boundary_corrections, &dg_formulation,
435 epsilon, &exterior_package_data,
436 &function_name_index, &interior_package_data,
437 &python_dg_boundary_terms_functions,
438 &python_module](
auto package_data_tag_v) {
439 INFO(
"Testing boundary terms");
440 using boundary_correction_tag =
441 tmpl::type_from<decltype(package_data_tag_v)>;
442 using ResultType =
typename boundary_correction_tag::type;
444 CAPTURE(python_module);
446 gsl::at(python_dg_boundary_terms_functions, function_name_index);
447 CAPTURE(python_function);
450 const typename boundary_correction_tag::type python_result =
451 pypp::call<ResultType>(
452 python_module, python_function,
453 get<DgPackageDataTags>(interior_package_data)...,
454 get<DgPackageDataTags>(exterior_package_data)...,
455 dg_formulation == ::dg::Formulation::StrongInertial);
457 get<boundary_correction_tag>(boundary_corrections), python_result,
458 Approx::custom().epsilon(epsilon).scale(1.0));
460 INFO(
"On line " << __LINE__ <<
" Python call to "
461 <<
gsl::at(python_dg_boundary_terms_functions,
463 <<
" in module " << python_module
464 <<
" failed: " << e.what());
467 ++function_name_index;
501 template <
typename System,
typename ConversionClassList = tmpl::list<>,
502 typename BoundaryCorrection,
size_t FaceDim,
typename... VolumeTags>
507 tmpl::size<typename BoundaryCorrection::dg_package_field_tags>::value>&
508 python_dg_package_data_functions,
511 tmpl::size<typename System::variables_tag::tags_list>::value>&
512 python_dg_boundary_terms_functions,
513 const BoundaryCorrection& correction,
const Mesh<FaceDim>& face_mesh,
515 const double epsilon = 1.0e-12) {
517 "All boundary correction classes must be marked `final`.");
518 using package_temporary_tags =
519 typename BoundaryCorrection::dg_package_data_temporary_tags;
520 using package_primitive_tags = detail::get_correction_primitive_vars<
521 System::has_primitive_and_conservative_vars, BoundaryCorrection>;
522 using variables_tags =
typename System::variables_tag::tags_list;
523 using flux_variables =
typename System::flux_variables;
528 for (
const auto use_moving_mesh : {
false,
true}) {
529 for (
const auto dg_formulation :
530 {::dg::Formulation::StrongInertial, ::dg::Formulation::WeakInertial}) {
531 detail::test_with_python<ConversionClassList, variables_tags>(
532 python_module, python_dg_package_data_functions,
533 python_dg_boundary_terms_functions, correction, face_mesh,
534 volume_data, use_moving_mesh, dg_formulation, epsilon,
535 tmpl::append<variables_tags, flux_tags, package_temporary_tags,
536 package_primitive_tags>{},
537 typename BoundaryCorrection::dg_package_field_tags{});