diff --git a/tests/Unit/Evolution/DiscontinuousGalerkin/BoundaryCorrectionsHelper.py b/tests/Unit/Evolution/DiscontinuousGalerkin/BoundaryCorrectionsHelper.py
new file mode 100644
index 000000000000..9994fb1fafe6
--- /dev/null
+++ b/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))
diff --git a/tests/Unit/Evolution/DiscontinuousGalerkin/Test_BoundaryCorrectionsHelper.cpp b/tests/Unit/Evolution/DiscontinuousGalerkin/Test_BoundaryCorrectionsHelper.cpp
index 3ec5dc33eaa9..19280fcd0079 100644
--- a/tests/Unit/Evolution/DiscontinuousGalerkin/Test_BoundaryCorrectionsHelper.cpp
+++ b/tests/Unit/Evolution/DiscontinuousGalerkin/Test_BoundaryCorrectionsHelper.cpp
@@ -11,6 +11,7 @@
 #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"
@@ -18,6 +19,33 @@
 #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>();
 }
diff --git a/tests/Unit/Helpers/Evolution/DiscontinuousGalerkin/BoundaryCorrections.hpp b/tests/Unit/Helpers/Evolution/DiscontinuousGalerkin/BoundaryCorrections.hpp
index 4646ebec0a1b..2fc2667f028b 100644
--- a/tests/Unit/Helpers/Evolution/DiscontinuousGalerkin/BoundaryCorrections.hpp
+++ b/tests/Unit/Helpers/Evolution/DiscontinuousGalerkin/BoundaryCorrections.hpp
@@ -5,8 +5,8 @@
 
 #include "Framework/TestingFramework.hpp"
 
-#include <boost/optional.hpp>
 #include <cstddef>
+#include <optional>
 #include <random>
 #include <string>
 #include <type_traits>
@@ -18,6 +18,7 @@
 #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
 #include "DataStructures/Tensor/Tensor.hpp"
 #include "DataStructures/Variables.hpp"
+#include "Framework/Pypp.hpp"
 #include "Framework/TestHelpers.hpp"
 #include "Helpers/DataStructures/MakeWithRandomValues.hpp"
 #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
@@ -73,10 +74,10 @@ void call_dg_package_data(
     const Variables<tmpl::list<FaceTags...>>& face_variables,
     const tuples::TaggedTuple<VolumeTags...>& volume_data,
     const tnsr::i<DataVector, Dim, Frame::Inertial>& unit_normal_covector,
-    const boost::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
+    const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
         mesh_velocity) {
-  boost::optional<Scalar<DataVector>> normal_dot_mesh_velocity{};
-  if (static_cast<bool>(mesh_velocity)) {
+  std::optional<Scalar<DataVector>> normal_dot_mesh_velocity{};
+  if (mesh_velocity.has_value()) {
     normal_dot_mesh_velocity =
         dot_product(*mesh_velocity, unit_normal_covector);
   }
@@ -154,7 +155,7 @@ void test_boundary_correction_impl(
   std::uniform_real_distribution<> pos_neg_dist(-1.0, 1.0);
   DataVector used_for_size{face_mesh.number_of_grid_points()};
 
-  boost::optional<tnsr::I<DataVector, FaceDim + 1, Frame::Inertial>>
+  std::optional<tnsr::I<DataVector, FaceDim + 1, Frame::Inertial>>
       mesh_velocity{};
   if (use_moving_mesh) {
     mesh_velocity = make_with_random_values<
@@ -315,4 +316,226 @@ void test_boundary_correction_conservation(
     }
   }
 }
+
+namespace detail {
+template <typename ConversionClassList, typename VariablesTags,
+          typename BoundaryCorrection, size_t FaceDim, typename... FaceTags,
+          typename... VolumeTags, typename... DgPackageDataTags>
+void test_with_python(
+    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<VariablesTags>::value>&
+        python_dg_boundary_terms_functions,
+    const BoundaryCorrection& correction, const Mesh<FaceDim>& face_mesh,
+    const tuples::TaggedTuple<VolumeTags...>& volume_data,
+    const bool use_moving_mesh, const ::dg::Formulation dg_formulation,
+    const double epsilon, tmpl::list<FaceTags...> /*meta*/,
+    tmpl::list<DgPackageDataTags...> /*meta*/) {
+  CAPTURE(face_mesh);
+  CAPTURE(dg_formulation);
+  CAPTURE(use_moving_mesh);
+  REQUIRE(face_mesh.number_of_grid_points() >= 1);
+  using dg_package_field_tags =
+      typename BoundaryCorrection::dg_package_field_tags;
+
+  MAKE_GENERATOR(gen);
+  std::uniform_real_distribution<> dist(0.0, 1.0);
+  DataVector used_for_size{face_mesh.number_of_grid_points()};
+
+  Variables<dg_package_field_tags> package_data{used_for_size.size()};
+  const auto fields_on_face =
+      make_with_random_values<Variables<tmpl::list<FaceTags...>>>(
+          make_not_null(&gen), make_not_null(&dist), used_for_size);
+  auto unit_normal_covector = make_with_random_values<
+      tnsr::i<DataVector, FaceDim + 1, Frame::Inertial>>(
+      make_not_null(&gen), make_not_null(&dist), used_for_size);
+  const auto normal_magnitude = magnitude(unit_normal_covector);
+  for (auto& component : unit_normal_covector) {
+    component /= get(normal_magnitude);
+  }
+  std::optional<tnsr::I<DataVector, FaceDim + 1, Frame::Inertial>>
+      mesh_velocity{};
+  if (use_moving_mesh) {
+    mesh_velocity = make_with_random_values<
+        tnsr::I<DataVector, FaceDim + 1, Frame::Inertial>>(
+        make_not_null(&gen), make_not_null(&dist), used_for_size);
+  }
+  std::optional<Scalar<DataVector>> normal_dot_mesh_velocity{};
+  if (mesh_velocity.has_value()) {
+    normal_dot_mesh_velocity =
+        dot_product(*mesh_velocity, unit_normal_covector);
+  }
+
+  // Call C++ implementation of dg_package_data
+  call_dg_package_data(make_not_null(&package_data), correction, fields_on_face,
+                       volume_data, unit_normal_covector, mesh_velocity);
+
+  // Call python implementation of dg_package_data
+  size_t function_name_index = 0;
+  tmpl::for_each<dg_package_field_tags>(
+      [epsilon, &fields_on_face, &function_name_index, &mesh_velocity,
+       &normal_dot_mesh_velocity, &package_data,
+       &python_dg_package_data_functions, &python_module, &unit_normal_covector,
+       &volume_data](auto package_data_tag_v) {
+        // avoid compiler warnings if there isn't any volume data
+        (void)volume_data;
+        INFO("Testing package data");
+        using package_data_tag = tmpl::type_from<decltype(package_data_tag_v)>;
+        using ResultType = typename package_data_tag::type;
+        try {
+          CAPTURE(python_module);
+          CAPTURE(
+              gsl::at(python_dg_package_data_functions, function_name_index));
+          const auto python_result =
+              pypp::call<ResultType, ConversionClassList>(
+                  python_module,
+                  gsl::at(python_dg_package_data_functions,
+                          function_name_index),
+                  get<FaceTags>(fields_on_face)..., unit_normal_covector,
+                  mesh_velocity, normal_dot_mesh_velocity,
+                  get<VolumeTags>(volume_data)...);
+          CHECK_ITERABLE_CUSTOM_APPROX(
+              get<package_data_tag>(package_data), python_result,
+              Approx::custom().epsilon(epsilon).scale(1.0));
+        } catch (const std::exception& e) {
+          INFO("On line " << __LINE__ << " Python call to "
+                          << gsl::at(python_dg_package_data_functions,
+                                     function_name_index)
+                          << " in module " << python_module
+                          << " failed: " << e.what());
+          REQUIRE(false);
+        }
+        ++function_name_index;
+      });
+
+  // Now we need to check the dg_boundary_terms function.
+  const auto interior_package_data =
+      make_with_random_values<Variables<dg_package_field_tags>>(
+          make_not_null(&gen), make_not_null(&dist), used_for_size);
+  const auto exterior_package_data =
+      make_with_random_values<Variables<dg_package_field_tags>>(
+          make_not_null(&gen), make_not_null(&dist), used_for_size);
+  // We don't need to prefix the VariablesTags with anything because we are not
+  // interacting with any code that cares about what the tags are, just that the
+  // types matched the evolved variables.
+  Variables<VariablesTags> boundary_corrections{
+      face_mesh.number_of_grid_points()};
+
+  // Call C++ implementation of dg_boundary_terms
+  call_dg_boundary_terms(make_not_null(&boundary_corrections), correction,
+                         interior_package_data, exterior_package_data,
+                         dg_formulation);
+
+  // Call python implementation of dg_boundary_terms
+  function_name_index = 0;
+  tmpl::for_each<VariablesTags>([&boundary_corrections, &dg_formulation,
+                                 epsilon, &exterior_package_data,
+                                 &function_name_index, &interior_package_data,
+                                 &python_dg_boundary_terms_functions,
+                                 &python_module](auto package_data_tag_v) {
+    INFO("Testing boundary terms");
+    using boundary_correction_tag =
+        tmpl::type_from<decltype(package_data_tag_v)>;
+    using ResultType = typename boundary_correction_tag::type;
+    try {
+      CAPTURE(python_module);
+      const std::string& python_function =
+          gsl::at(python_dg_boundary_terms_functions, function_name_index);
+      CAPTURE(python_function);
+      // Make explicitly depend on tag type to avoid type deduction issues with
+      // GCC7
+      const typename boundary_correction_tag::type python_result =
+          pypp::call<ResultType>(
+              python_module, python_function,
+              get<DgPackageDataTags>(interior_package_data)...,
+              get<DgPackageDataTags>(exterior_package_data)...,
+              dg_formulation == ::dg::Formulation::StrongInertial);
+      CHECK_ITERABLE_CUSTOM_APPROX(
+          get<boundary_correction_tag>(boundary_corrections), python_result,
+          Approx::custom().epsilon(epsilon).scale(1.0));
+    } catch (const std::exception& e) {
+      INFO("On line " << __LINE__ << " Python call to "
+                      << gsl::at(python_dg_boundary_terms_functions,
+                                 function_name_index)
+                      << " in module " << python_module
+                      << " failed: " << e.what());
+      REQUIRE(false);
+    }
+    ++function_name_index;
+  });
+}
+}  // namespace detail
+
+/*!
+ * \ingroup TestingFrameworkGroup
+ * \brief Tests that the `dg_package_data` and `dg_boundary_terms` functions
+ * agree with the python implementation.
+ *
+ * The variables are filled with random numbers between zero and one before
+ * being passed to the implementations. If in the future we need support for
+ * negative numbers we can add the ability to specify a single range for all
+ * random numbers or each individually.
+ *
+ * Please note the following:
+ * - The `pypp::SetupLocalPythonEnvironment` must be created before the
+ *   `test_boundary_correction_with_python` can be called.
+ * - The `dg_formulation` is passed as a bool `use_strong_form` to the python
+ *   functions since we don't want to rely on python bindings for the enum.
+ * - You can convert custom types using the `ConversionClassList` template
+ *   parameter, which is then passed to `pypp::call()`. This allows you to,
+ *   e.g., convert an equation of state into an array locally in a test file.
+ * - There must be one python function to compute the packaged data for each tag
+ *   in `dg_package_field_tags`
+ * - There must be one python function to compute the boundary correction for
+ *   each tag in `System::variables_tag`
+ * - The arguments to the python functions for computing the packaged data are
+ *   the same as the arguments for the C++ `dg_package_data` function, excluding
+ *   the `gsl::not_null` arguments.
+ * - The arguments to the python functions for computing the boundary
+ *   corrections are the same as the arguments for the C++ `dg_boundary_terms`
+ *   function, excluding the `gsl::not_null` arguments.
+ */
+template <typename System, typename ConversionClassList = tmpl::list<>,
+          typename BoundaryCorrection, size_t FaceDim, typename... VolumeTags>
+void test_boundary_correction_with_python(
+    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 double epsilon = 1.0e-12) {
+  static_assert(std::is_final_v<std::decay_t<BoundaryCorrection>>,
+                "All boundary correction classes must be marked `final`.");
+  using package_temporary_tags =
+      typename BoundaryCorrection::dg_package_data_temporary_tags;
+  using package_primitive_tags = detail::get_correction_primitive_vars<
+      System::has_primitive_and_conservative_vars, BoundaryCorrection>;
+  using variables_tags = typename System::variables_tag::tags_list;
+  using flux_variables = typename System::flux_variables;
+  using flux_tags =
+      db::wrap_tags_in<::Tags::Flux, flux_variables, tmpl::size_t<FaceDim + 1>,
+                       Frame::Inertial>;
+
+  for (const auto use_moving_mesh : {false, true}) {
+    for (const auto dg_formulation :
+         {::dg::Formulation::StrongInertial, ::dg::Formulation::WeakInertial}) {
+      detail::test_with_python<ConversionClassList, variables_tags>(
+          python_module, python_dg_package_data_functions,
+          python_dg_boundary_terms_functions, correction, face_mesh,
+          volume_data, use_moving_mesh, dg_formulation, epsilon,
+          tmpl::append<variables_tags, flux_tags, package_temporary_tags,
+                       package_primitive_tags>{},
+          typename BoundaryCorrection::dg_package_field_tags{});
+    }
+  }
+}
 }  // namespace TestHelpers::evolution::dg