Merge pull request #2700 from nilsdeppe/boundary_correction_random

Add python test function for boundary corrections

tests/Unit/Evolution/DiscontinuousGalerkin/BoundaryCorrectionsHelper.py

@@ -0,0 +1,71 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+import numpy as np
+
+
+def dg_package_data_var1(var1, var2, flux_var1, flux_var2, normal_covector,
+                         mesh_velocity, normal_dot_mesh_velocity,
+                         volume_double):
+    return var1
+
+
+def dg_package_data_var1_normal_dot_flux(var1, var2, flux_var1, flux_var2,
+                                         normal_covector, mesh_velocity,
+                                         normal_dot_mesh_velocity,
+                                         volume_double):
+    return np.einsum("i,i->", flux_var1, normal_covector)
+
+
+def dg_package_data_var2(var1, var2, flux_var1, flux_var2, normal_covector,
+                         mesh_velocity, normal_dot_mesh_velocity,
+                         volume_double):
+    return var2
+
+
+def dg_package_data_var2_normal_dot_flux(var1, var2, flux_var1, flux_var2,
+                                         normal_covector, mesh_velocity,
+                                         normal_dot_mesh_velocity,
+                                         volume_double):
+    return np.einsum("ij,j->i", flux_var2, normal_covector)
+
+
+def dg_package_data_abs_char_speed(var1, var2, flux_var1, flux_var2,
+                                   normal_covector, mesh_velocity,
+                                   normal_dot_mesh_velocity, volume_double):
+    if not isinstance(volume_double, float):
+        volume_double = volume_double[0]
+    if normal_dot_mesh_velocity is None:
+        return np.abs(volume_double * var1)
+    else:
+        return np.abs(volume_double * var1 - normal_dot_mesh_velocity)
+
+
+def dg_boundary_terms_var1(var1_int, normal_dot_flux_var1_int, var2_int,
+                           normal_dot_flux_var2_int, abs_char_speed_int,
+                           var1_ext, normal_dot_flux_var1_ext, var2_ext,
+                           normal_dot_flux_var2_ext, abs_char_speed_ext,
+                           use_strong_form):
+    if use_strong_form:
+        return (-0.5 * (normal_dot_flux_var1_int + normal_dot_flux_var1_ext) -
+                0.5 * np.maximum(abs_char_speed_int, abs_char_speed_ext) *
+                (var1_ext - var1_int))
+    else:
+        return (0.5 * (normal_dot_flux_var1_int - normal_dot_flux_var1_ext) -
+                0.5 * np.maximum(abs_char_speed_int, abs_char_speed_ext) *
+                (var1_ext - var1_int))
+
+
+def dg_boundary_terms_var2(var1_int, normal_dot_flux_var1_int, var2_int,
+                           normal_dot_flux_var2_int, abs_char_speed_int,
+                           var1_ext, normal_dot_flux_var1_ext, var2_ext,
+                           normal_dot_flux_var2_ext, abs_char_speed_ext,
+                           use_strong_form):
+    if use_strong_form:
+        return (-0.5 * (normal_dot_flux_var2_int + normal_dot_flux_var2_ext) -
+                0.5 * np.maximum(abs_char_speed_int, abs_char_speed_ext) *
+                (var2_ext - var2_int))
+    else:
+        return (0.5 * (normal_dot_flux_var2_int - normal_dot_flux_var2_ext) -
+                0.5 * np.maximum(abs_char_speed_int, abs_char_speed_ext) *
+                (var2_ext - var2_int))

tests/Unit/Evolution/DiscontinuousGalerkin/Test_BoundaryCorrectionsHelper.cpp

@@ -11,13 +11,41 @@
 #include "DataStructures/Tensor/Tensor.hpp"
 #include "DataStructures/Variables.hpp"
 #include "DataStructures/VariablesTag.hpp"
+#include "Framework/SetupLocalPythonEnvironment.hpp"
 #include "Helpers/DataStructures/MakeWithRandomValues.hpp"
 #include "Helpers/Evolution/DiscontinuousGalerkin/BoundaryCorrections.hpp"
 #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
 #include "Utilities/Gsl.hpp"
 #include "Utilities/TMPL.hpp"
 
 namespace {
+struct VolumeDouble {
+  double value;
+};
+
+struct VolumeDoubleConversion {
+  // convert to a std::array to test non-trivial type conversion
+  using unpacked_container = std::array<double, 1>;
+  using packed_container = VolumeDouble;
+  using packed_type = double;
+
+  static inline unpacked_container unpack(
+      const packed_container packed,
+      const size_t /*grid_point_index*/) noexcept {
+    return {{packed.value}};
+  }
+
+  static inline void pack(const gsl::not_null<packed_container*> packed,
+                          const unpacked_container unpacked,
+                          const size_t /*grid_point_index*/) {
+    packed->value = unpacked[0];
+  }
+
+  static inline size_t get_size(const packed_container& /*packed*/) noexcept {
+    return 1;
+  }
+};
+
 namespace Tags {
 struct Var1 : db::SimpleTag {
   using type = Scalar<DataVector>;
@@ -28,8 +56,9 @@ struct Var2 : db::SimpleTag {
   using type = tnsr::i<DataVector, Dim, Frame::Inertial>;
 };
 
+template <typename Type>
 struct VolumeDouble : db::SimpleTag {
-  using type = double;
+  using type = Type;
 };
 }  // namespace Tags
 
@@ -52,8 +81,8 @@ struct System {
   using compute_volume_time_derivative_terms = TimeDerivativeTerms;
 };
 
-template <size_t Dim>
-struct Correction {
+template <size_t Dim, typename VolumeDoubleType>
+struct Correction final {
  private:
   struct AbsCharSpeed : db::SimpleTag {
     using type = Scalar<DataVector>;
@@ -64,7 +93,8 @@ struct Correction {
       tmpl::list<Tags::Var1, ::Tags::NormalDotFlux<Tags::Var1>, Tags::Var2<Dim>,
                  ::Tags::NormalDotFlux<Tags::Var2<Dim>>, AbsCharSpeed>;
   using dg_package_data_temporary_tags = tmpl::list<>;
-  using dg_package_data_volume_tags = tmpl::list<Tags::VolumeDouble>;
+  using dg_package_data_volume_tags =
+      tmpl::list<Tags::VolumeDouble<VolumeDoubleType>>;
 
   double dg_package_data(
       const gsl::not_null<Scalar<DataVector>*> packaged_var1,
@@ -79,10 +109,16 @@ struct Correction {
       const tnsr::I<DataVector, Dim, Frame::Inertial>& flux_var1,
       const tnsr::Ij<DataVector, Dim, Frame::Inertial>& flux_var2,
       const tnsr::i<DataVector, Dim, Frame::Inertial>& normal_covector,
-      const boost::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
+      const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
       /*mesh_velocity*/,
-      const boost::optional<Scalar<DataVector>>& normal_dot_mesh_velocity,
-      const double volume_double) const noexcept {
+      const std::optional<Scalar<DataVector>>& normal_dot_mesh_velocity,
+      const VolumeDoubleType volume_double_in) const noexcept {
+    double volume_double = 0.0;
+    if constexpr (std::is_same_v<double, VolumeDoubleType>) {
+      volume_double = volume_double_in;
+    } else {
+      volume_double = volume_double_in.value;
+    }
     *packaged_var1 = var1;
     *packaged_var2 = var2;
     dot_product(packaged_normal_dot_flux_var1, flux_var1, normal_covector);
@@ -150,20 +186,37 @@ struct Correction {
   }
 };
 
-template <size_t Dim>
+template <size_t Dim, typename VolumeDoubleType>
 void test() {
-  const Correction<Dim> correction{};
+  const Correction<Dim, VolumeDoubleType> correction{};
   const Mesh<Dim - 1> face_mesh{Dim * Dim, Spectral::Basis::Legendre,
                                 Spectral::Quadrature::Gauss};
   TestHelpers::evolution::dg::test_boundary_correction_conservation<
       System<Dim>>(correction, face_mesh,
-                   tuples::TaggedTuple<Tags::VolumeDouble>{2.3});
+                   tuples::TaggedTuple<Tags::VolumeDouble<VolumeDoubleType>>{
+                       VolumeDoubleType{2.3}});
+  TestHelpers::evolution::dg::test_boundary_correction_with_python<
+      System<Dim>, tmpl::list<VolumeDoubleConversion>>(
+      "BoundaryCorrectionsHelper",
+      {{"dg_package_data_var1", "dg_package_data_var1_normal_dot_flux",
+        "dg_package_data_var2", "dg_package_data_var2_normal_dot_flux",
+        "dg_package_data_abs_char_speed"}},
+      {{"dg_boundary_terms_var1", "dg_boundary_terms_var2"}}, correction,
+      face_mesh,
+      tuples::TaggedTuple<Tags::VolumeDouble<VolumeDoubleType>>{
+          VolumeDoubleType{2.3}});
 }
 }  // namespace
 
 SPECTRE_TEST_CASE("Unit.Evolution.DG.BoundaryCorrectionsHelper",
                   "[Unit][Evolution]") {
-  test<1>();
-  test<2>();
-  test<3>();
+  pypp::SetupLocalPythonEnvironment local_python_env{
+      "Evolution/DiscontinuousGalerkin/"};
+  test<1, double>();
+  test<2, double>();
+  test<3, double>();
+
+  test<1, VolumeDouble>();
+  test<2, VolumeDouble>();
+  test<3, VolumeDouble>();
 }