diff --git a/src/Elliptic/DiscontinuousGalerkin/DgOperator.hpp b/src/Elliptic/DiscontinuousGalerkin/DgOperator.hpp
index bfc426fd5fb2..7a60689a3e19 100644
--- a/src/Elliptic/DiscontinuousGalerkin/DgOperator.hpp
+++ b/src/Elliptic/DiscontinuousGalerkin/DgOperator.hpp
@@ -626,7 +626,18 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
         // This is the sign flip that makes the operator _minus_ the Laplacian
         // for a Poisson system
         *operator_applied_to_vars *= -1.;
-      } else {
+      } else if (formulation == ::dg::Formulation::StrongLogical) {
+        // Strong divergence but with the Jacobian moved into the divergence:
+        //   div(F) = 1/J_p \sum_q (D_\hat{i})_pq J_q (J^\hat{i}_i)_q (F^i)_q.
+        const auto logical_fluxes = transform::first_index_to_different_frame(
+            primal_fluxes, det_times_inv_jacobian);
+        logical_divergence(operator_applied_to_vars, logical_fluxes, mesh);
+        if (massive) {
+          *operator_applied_to_vars *= -1.;
+        } else {
+          *operator_applied_to_vars *= -1. * get(det_inv_jacobian);
+        }
+      } else if (formulation == ::dg::Formulation::WeakInertial) {
         // Compute weak divergence:
         //   F^i \partial_i \phi = 1/w_p \sum_q
         //     (D^T_\hat{i})_pq (w det(J) J^\hat{i}_i F^i)_q
@@ -635,6 +646,11 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
         if (not massive) {
           *operator_applied_to_vars *= get(det_inv_jacobian);
         }
+      } else {
+        ERROR("Unsupported DG formulation: "
+              << formulation
+              << "\nSupported formulations are: StrongInertial, WeakInertial, "
+                 "StrongLogical.");
       }
       if constexpr (not std::is_same_v<SourcesComputer, void>) {
         Variables<tmpl::list<OperatorTags...>> sources{num_points, 0.};
@@ -803,10 +819,13 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
           lhs[i] += 0.5 * rhs[i];
         }
       };
+      const bool is_strong_formulation =
+          formulation == ::dg::Formulation::StrongInertial or
+          formulation == ::dg::Formulation::StrongLogical;
       EXPAND_PACK_LEFT_TO_RIGHT(add_avg_contribution(
           get<::Tags::NormalDotFlux<PrimalMortarVars>>(local_data.field_data),
           get<::Tags::NormalDotFlux<PrimalMortarVars>>(remote_data.field_data),
-          formulation == ::dg::Formulation::StrongInertial ? 0.5 : -0.5));
+          is_strong_formulation ? 0.5 : -0.5));
 
       // Project from the mortar back down to the face if needed, lift and add
       // to operator. See auxiliary boundary corrections above for details.
diff --git a/src/Evolution/DiscontinuousGalerkin/Actions/VolumeTermsImpl.tpp b/src/Evolution/DiscontinuousGalerkin/Actions/VolumeTermsImpl.tpp
index 8ecea3628859..d7a7b55c2549 100644
--- a/src/Evolution/DiscontinuousGalerkin/Actions/VolumeTermsImpl.tpp
+++ b/src/Evolution/DiscontinuousGalerkin/Actions/VolumeTermsImpl.tpp
@@ -263,7 +263,12 @@ void volume_terms(
              "when using the weak form.");
       (*div_fluxes) *= get(*det_inverse_jacobian);
     } else {
-      ERROR("Unsupported DG formulation: " << dg_formulation);
+      // Note: when adding support for other DG formulations here, also check
+      // that the implementations of `BoundaryCorrection` classes support the
+      // added formulation.
+      ERROR("Unsupported DG formulation: "
+            << dg_formulation
+            << "\nSupported formulations are: StrongInertial, WeakInertial.");
     }
     tmpl::for_each<flux_variables>(
         [&dg_formulation, &dt_vars_ptr, &div_fluxes](auto var_tag_v) {
diff --git a/src/Evolution/Executables/GeneralizedHarmonic/EvolveGhBinaryBlackHole.hpp b/src/Evolution/Executables/GeneralizedHarmonic/EvolveGhBinaryBlackHole.hpp
index c5044cbb89b0..b098a2548761 100644
--- a/src/Evolution/Executables/GeneralizedHarmonic/EvolveGhBinaryBlackHole.hpp
+++ b/src/Evolution/Executables/GeneralizedHarmonic/EvolveGhBinaryBlackHole.hpp
@@ -239,8 +239,6 @@ struct EvolutionMetavars {
   static constexpr size_t volume_dim = 3;
   static constexpr bool use_damped_harmonic_rollon = false;
   using system = gh::System<volume_dim>;
-  static constexpr dg::Formulation dg_formulation =
-      dg::Formulation::StrongInertial;
   using temporal_id = Tags::TimeStepId;
   using TimeStepperBase = LtsTimeStepper;
 
diff --git a/src/Evolution/Executables/GrMhd/ValenciaDivClean/EvolveValenciaDivClean.hpp b/src/Evolution/Executables/GrMhd/ValenciaDivClean/EvolveValenciaDivClean.hpp
index a43c5e7959e9..2dd1d2acbb0e 100644
--- a/src/Evolution/Executables/GrMhd/ValenciaDivClean/EvolveValenciaDivClean.hpp
+++ b/src/Evolution/Executables/GrMhd/ValenciaDivClean/EvolveValenciaDivClean.hpp
@@ -211,8 +211,6 @@ struct EvolutionMetavars<tmpl::list<InterpolationTargetTags...>,
       UseParametrizedDeleptonization;
   static constexpr bool use_dg_subcell = true;
   static constexpr size_t volume_dim = 3;
-  static constexpr dg::Formulation dg_formulation =
-      dg::Formulation::StrongInertial;
   using initial_data_list =
       grmhd::ValenciaDivClean::InitialData::initial_data_list;
 
diff --git a/src/Evolution/Executables/NewtonianEuler/EvolveNewtonianEuler.hpp b/src/Evolution/Executables/NewtonianEuler/EvolveNewtonianEuler.hpp
index 26ae4d44826a..1e108d5b3a4b 100644
--- a/src/Evolution/Executables/NewtonianEuler/EvolveNewtonianEuler.hpp
+++ b/src/Evolution/Executables/NewtonianEuler/EvolveNewtonianEuler.hpp
@@ -133,8 +133,6 @@ struct EvolutionMetavars {
   // The use_dg_subcell flag controls whether to use "standard" limiting (false)
   // or a DG-FD hybrid scheme (true).
   static constexpr bool use_dg_subcell = true;
-  static constexpr dg::Formulation dg_formulation =
-      dg::Formulation::StrongInertial;
 
   using system = NewtonianEuler::System<Dim>;
 
diff --git a/src/Evolution/Executables/RadiationTransport/M1Grey/EvolveM1Grey.hpp b/src/Evolution/Executables/RadiationTransport/M1Grey/EvolveM1Grey.hpp
index ff36ebd37512..6943603c49c7 100644
--- a/src/Evolution/Executables/RadiationTransport/M1Grey/EvolveM1Grey.hpp
+++ b/src/Evolution/Executables/RadiationTransport/M1Grey/EvolveM1Grey.hpp
@@ -102,8 +102,6 @@ class CProxy_GlobalCache;
 
 struct EvolutionMetavars {
   static constexpr size_t volume_dim = 3;
-  static constexpr dg::Formulation dg_formulation =
-      dg::Formulation::StrongInertial;
 
   // To switch which initial data is evolved you only need to change the
   // line `using initial_data = ...;` and include the header file for the
diff --git a/src/Evolution/Executables/RelativisticEuler/Valencia/EvolveValencia.hpp b/src/Evolution/Executables/RelativisticEuler/Valencia/EvolveValencia.hpp
index 872dd22b7b65..f28a4efa7d0c 100644
--- a/src/Evolution/Executables/RelativisticEuler/Valencia/EvolveValencia.hpp
+++ b/src/Evolution/Executables/RelativisticEuler/Valencia/EvolveValencia.hpp
@@ -108,8 +108,6 @@ class CProxy_GlobalCache;
 template <size_t Dim, typename InitialData>
 struct EvolutionMetavars {
   static constexpr size_t volume_dim = Dim;
-  static constexpr dg::Formulation dg_formulation =
-      dg::Formulation::StrongInertial;
 
   using initial_data = InitialData;
   static_assert(
diff --git a/src/Evolution/Executables/ScalarWave/EvolveScalarWave.hpp b/src/Evolution/Executables/ScalarWave/EvolveScalarWave.hpp
index 5f8e9eb03597..762273d809b1 100644
--- a/src/Evolution/Executables/ScalarWave/EvolveScalarWave.hpp
+++ b/src/Evolution/Executables/ScalarWave/EvolveScalarWave.hpp
@@ -131,8 +131,6 @@ struct EvolutionMetavars {
   using initial_data_list = ScalarWave::Solutions::all_solutions<Dim>;
 
   using system = ScalarWave::System<Dim>;
-  static constexpr dg::Formulation dg_formulation =
-      dg::Formulation::StrongInertial;
   using temporal_id = Tags::TimeStepId;
   using TimeStepperBase = TimeStepper;
 
diff --git a/src/NumericalAlgorithms/DiscontinuousGalerkin/Formulation.cpp b/src/NumericalAlgorithms/DiscontinuousGalerkin/Formulation.cpp
index 17f699e7bdfb..111787ae3d27 100644
--- a/src/NumericalAlgorithms/DiscontinuousGalerkin/Formulation.cpp
+++ b/src/NumericalAlgorithms/DiscontinuousGalerkin/Formulation.cpp
@@ -16,6 +16,8 @@ std::ostream& operator<<(std::ostream& os, const Formulation t) {
       return os << "StrongInertial";
     case Formulation::WeakInertial:
       return os << "WeakInertial";
+    case Formulation::StrongLogical:
+      return os << "StrongLogical";
     default:
       ERROR("Unknown DG formulation.");
   }
@@ -30,9 +32,12 @@ dg::Formulation Options::create_from_yaml<dg::Formulation>::create<void>(
     return dg::Formulation::StrongInertial;
   } else if ("WeakInertial" == type_read) {
     return dg::Formulation::WeakInertial;
+  } else if ("StrongLogical" == type_read) {
+    return dg::Formulation::StrongLogical;
   }
-  PARSE_ERROR(options.context(), "Failed to convert \""
-                                     << type_read
-                                     << "\" to dg::Formulation. Must be one "
-                                        "of StrongInertial or WeakInertial.");
+  PARSE_ERROR(options.context(),
+              "Failed to convert \""
+                  << type_read
+                  << "\" to dg::Formulation. Must be one of "
+                     "StrongInertial, WeakInertial, or StrongLogical.");
 }
diff --git a/src/NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp b/src/NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp
index ffef79b1b9c7..adb63d3a4688 100644
--- a/src/NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp
+++ b/src/NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp
@@ -30,8 +30,28 @@ namespace dg {
  *   the "weak" part refers to the fact that the boundary correction terms are
  *   non-zero even if the solution is continuous at the interfaces.
  *   See \cite Teukolsky2015ega for an overview.
+ * - The `StrongLogical` formulation is also known as the transform then
+ *   integrate formulation. The "logical" part of the name refers to the fact
+ *   that the integration is done over the logical coordinates, while the
+ *   "strong" part refers to the fact that the boundary correction terms are
+ *   zero if the solution is continuous at the interfaces. This formulation
+ *   arises from the `StrongInertial` formulation by moving the Jacobians that
+ *   appear when computing the divergence of fluxes into the divergence so the
+ *   divergence is computed in logical coordinates:
+ *   \begin{equation}
+ *     \partial_i F^i = \frac{1}{J} \partial_\hat{\imath} J F^\hat{\imath}
+ *   \end{equation}
+ *   where $J$ is the Jacobian determinant and $\hat{\imath}$ are the logical
+ *   coordinates. This is possible because of the "metric identities"
+ *   \begin{equation}
+ *    \partial_\hat{\imath} \left(J \frac{\partial \xi^\hat{\imath}}
+ *    {\partial x^i}\right) = 0.
+ *   \end{equation}
+ *   See also `dg::metric_identity_det_jac_times_inv_jac` for details and for
+ *   functions that compute the Jacobians so they satisfy the metric identities
+ *   numerically (which may or may not be useful or necessary).
  */
-enum class Formulation { StrongInertial, WeakInertial };
+enum class Formulation { StrongInertial, WeakInertial, StrongLogical };
 
 std::ostream& operator<<(std::ostream& os, Formulation t);
 }  // namespace dg
diff --git a/src/NumericalAlgorithms/LinearOperators/Divergence.cpp b/src/NumericalAlgorithms/LinearOperators/Divergence.cpp
index b5cba1a54c47..09185380fcfa 100644
--- a/src/NumericalAlgorithms/LinearOperators/Divergence.cpp
+++ b/src/NumericalAlgorithms/LinearOperators/Divergence.cpp
@@ -6,6 +6,7 @@
 #include <array>
 #include <cstddef>
 
+#include "DataStructures/ApplyMatrices.hpp"
 #include "DataStructures/ComplexDataVector.hpp"
 #include "DataStructures/DataBox/Tag.hpp"
 #include "DataStructures/DataVector.hpp"
@@ -61,6 +62,26 @@ void divergence(const gsl::not_null<Scalar<DataType>*> div_input,
   }
 }
 
+template <typename ResultTensor, typename FluxTensor, size_t Dim>
+void logical_divergence(const gsl::not_null<ResultTensor*> div_flux,
+                        const FluxTensor& flux, const Mesh<Dim>& mesh) {
+  // Note: This function hasn't been optimized much at all. Feel free to
+  // optimize if needed!
+  static const Matrix identity_matrix{};
+  for (size_t d = 0; d < Dim; ++d) {
+    auto matrices = make_array<Dim>(std::cref(identity_matrix));
+    gsl::at(matrices, d) =
+        Spectral::differentiation_matrix(mesh.slice_through(d));
+    for (size_t storage_index = 0; storage_index < div_flux->size();
+         ++storage_index) {
+      const auto div_flux_index = div_flux->get_tensor_index(storage_index);
+      const auto flux_index = prepend(div_flux_index, d);
+      div_flux->get(div_flux_index) +=
+          apply_matrices(matrices, flux.get(flux_index), mesh.extents());
+    }
+  }
+}
+
 #define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data)
 #define DIM(data) BOOST_PP_TUPLE_ELEM(1, data)
 #define FRAME(data) BOOST_PP_TUPLE_ELEM(2, data)
@@ -85,3 +106,28 @@ GENERATE_INSTANTIATIONS(INSTANTIATE, (DataVector, ComplexDataVector), (1, 2, 3),
 #undef DIM
 #undef FRAME
 #undef INSTANTIATE
+
+#define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data)
+#define DIM(data) BOOST_PP_TUPLE_ELEM(1, data)
+#define TENSOR(data) BOOST_PP_TUPLE_ELEM(2, data)
+
+#define INSTANTIATION_SCALAR(r, data)                                     \
+  template void logical_divergence(                                       \
+      const gsl::not_null<Scalar<DTYPE(data)>*> div_flux,                 \
+      const tnsr::I<DTYPE(data), DIM(data), Frame::ElementLogical>& flux, \
+      const Mesh<DIM(data)>& mesh);
+#define INSTANTIATION_TENSOR(r, data)                                    \
+  template void logical_divergence(                                      \
+      const gsl::not_null<tnsr::TENSOR(data) < DTYPE(data), DIM(data),   \
+                          Frame::Inertial>* > div_flux,                  \
+      const TensorMetafunctions::prepend_spatial_index<                  \
+          tnsr::TENSOR(data) < DTYPE(data), DIM(data), Frame::Inertial>, \
+      DIM(data), UpLo::Up, Frame::ElementLogical > &flux,                \
+      const Mesh<DIM(data)>& mesh);
+
+GENERATE_INSTANTIATIONS(INSTANTIATION_SCALAR, (DataVector, ComplexDataVector),
+                        (1, 2, 3))
+GENERATE_INSTANTIATIONS(INSTANTIATION_TENSOR, (DataVector, ComplexDataVector),
+                        (1, 2, 3), (i, I))
+
+#undef INSTANTIATION
diff --git a/src/NumericalAlgorithms/LinearOperators/Divergence.hpp b/src/NumericalAlgorithms/LinearOperators/Divergence.hpp
index 88a6f759cd1e..ef10009e2764 100644
--- a/src/NumericalAlgorithms/LinearOperators/Divergence.hpp
+++ b/src/NumericalAlgorithms/LinearOperators/Divergence.hpp
@@ -89,6 +89,29 @@ void divergence(gsl::not_null<Scalar<DataType>*> div_input,
                                       DerivativeFrame>& inverse_jacobian);
 /// @}
 
+/// @{
+/*!
+ * \brief Compute the divergence of fluxes in logical coordinates
+ *
+ * Applies the logical differentiation matrix to the fluxes in each dimension
+ * and sums over dimensions.
+ *
+ * \see divergence
+ */
+template <typename ResultTensor, typename FluxTensor, size_t Dim>
+void logical_divergence(gsl::not_null<ResultTensor*> div_flux,
+                        const FluxTensor& flux, const Mesh<Dim>& mesh);
+
+template <typename FluxTags, size_t Dim>
+auto logical_divergence(const Variables<FluxTags>& flux, const Mesh<Dim>& mesh)
+    -> Variables<db::wrap_tags_in<Tags::div, FluxTags>>;
+
+template <typename... DivTags, typename... FluxTags, size_t Dim>
+void logical_divergence(
+    gsl::not_null<Variables<tmpl::list<DivTags...>>*> div_flux,
+    const Variables<tmpl::list<FluxTags...>>& flux, const Mesh<Dim>& mesh);
+/// @}
+
 namespace Tags {
 /*!
  * \ingroup DataBoxTagsGroup
diff --git a/src/NumericalAlgorithms/LinearOperators/Divergence.tpp b/src/NumericalAlgorithms/LinearOperators/Divergence.tpp
index 295a861996a7..14f94ee87f42 100644
--- a/src/NumericalAlgorithms/LinearOperators/Divergence.tpp
+++ b/src/NumericalAlgorithms/LinearOperators/Divergence.tpp
@@ -96,3 +96,36 @@ void divergence(
   };
   EXPAND_PACK_LEFT_TO_RIGHT(apply_div(FluxTags{}, DivTags{}));
 }
+
+template <typename FluxTags, size_t Dim>
+Variables<db::wrap_tags_in<Tags::div, FluxTags>> logical_divergence(
+    const Variables<FluxTags>& flux, const Mesh<Dim>& mesh) {
+  Variables<db::wrap_tags_in<Tags::div, FluxTags>> div_flux(
+      flux.number_of_grid_points());
+  logical_divergence(make_not_null(&div_flux), flux, mesh);
+  return div_flux;
+}
+
+template <typename... ResultTags, typename... FluxTags, size_t Dim>
+void logical_divergence(
+    const gsl::not_null<Variables<tmpl::list<ResultTags...>>*> div_flux,
+    const Variables<tmpl::list<FluxTags...>>& flux, const Mesh<Dim>& mesh) {
+  static_assert(
+      (... and
+       std::is_same_v<tmpl::front<typename FluxTags::type::index_list>,
+                      SpatialIndex<Dim, UpLo::Up, Frame::ElementLogical>>),
+      "The first index of each flux must be an upper spatial index in "
+      "element-logical coordinates.");
+  static_assert(
+      (... and
+       std::is_same_v<
+           typename ResultTags::type,
+           TensorMetafunctions::remove_first_index<typename FluxTags::type>>),
+      "The result tensors must have the same type as the flux tensors with "
+      "their first index removed.");
+  div_flux->initialize(mesh.number_of_grid_points(), 0.);
+  // Note: This function hasn't been optimized much at all. Feel free to
+  // optimize if needed!
+  EXPAND_PACK_LEFT_TO_RIGHT(logical_divergence(
+      make_not_null(&get<ResultTags>(*div_flux)), get<FluxTags>(flux), mesh));
+}
diff --git a/tests/Unit/Elliptic/DiscontinuousGalerkin/CMakeLists.txt b/tests/Unit/Elliptic/DiscontinuousGalerkin/CMakeLists.txt
index 3efe6fde9355..75fc3f772917 100644
--- a/tests/Unit/Elliptic/DiscontinuousGalerkin/CMakeLists.txt
+++ b/tests/Unit/Elliptic/DiscontinuousGalerkin/CMakeLists.txt
@@ -7,6 +7,7 @@ set(LIBRARY_SOURCES
   Test_DgOperator.cpp
   Test_Penalty.cpp
   Test_Tags.cpp
+  Test_LargeOuterRadius.cpp
   )
 
 add_test_library(${LIBRARY} "${LIBRARY_SOURCES}")
diff --git a/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_DgOperator.cpp b/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_DgOperator.cpp
index 4f0340f479b6..181812ec7f31 100644
--- a/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_DgOperator.cpp
+++ b/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_DgOperator.cpp
@@ -266,6 +266,7 @@ template <
 void test_dg_operator(
     const DomainCreator<Dim>& domain_creator, const double penalty_parameter,
     const bool use_massive_dg_operator, const Spectral::Quadrature quadrature,
+    const ::dg::Formulation dg_formulation,
     const AnalyticSolution& analytic_solution,
     // NOLINTNEXTLINE(google-runtime-references)
     Approx& analytic_solution_aux_approx,
@@ -283,6 +284,8 @@ void test_dg_operator(
     const bool test_amr = false) {
   CAPTURE(penalty_parameter);
   CAPTURE(use_massive_dg_operator);
+  CAPTURE(quadrature);
+  CAPTURE(dg_formulation);
 
   using element_array = ElementArray;
   using vars_tag = typename element_array::vars_tag;
@@ -329,8 +332,8 @@ void test_dg_operator(
       logging::Tags::Verbosity<::amr::OptionTags::AmrGroup>>{
       std::move(domain), domain_creator.functions_of_time(),
       std::move(boundary_conditions), penalty_parameter,
-      use_massive_dg_operator, quadrature, ::dg::Formulation::StrongInertial,
-      analytic_solution, std::move(amr_criteria),
+      use_massive_dg_operator, quadrature, dg_formulation, analytic_solution,
+      std::move(amr_criteria),
       ::amr::Policies{::amr::Isotropy::Anisotropic, ::amr::Limits{}, true},
       ::Verbosity::Debug}};
 
@@ -579,7 +582,7 @@ void test_dg_operator(
 
 }  // namespace
 
-// [[TimeOut, 10]]
+// [[TimeOut, 30]]
 SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") {
   domain::creators::register_derived_with_charm();
   // This is what the tests below check:
@@ -673,6 +676,19 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") {
           expected_operator_vars_rnd_right)) =
           DataVector{-215.34463496424186, -37.924582769790135,
                      2.1993037260176997, 51.85418137198471};
+      const std::vector<
+          std::tuple<std::unordered_map<ElementId<1>, Vars>,
+                     std::unordered_map<ElementId<1>, PrimalFluxes>,
+                     std::unordered_map<ElementId<1>, OperatorVars>>>
+          regression_test_data{
+              {{{left_id, std::move(vars_rnd_left)},
+                {midleft_id, std::move(vars_rnd_midleft)},
+                {midright_id, std::move(vars_rnd_midright)},
+                {right_id, std::move(vars_rnd_right)}},
+               {{left_id, std::move(expected_primal_fluxes_rnd_left)},
+                {right_id, std::move(expected_primal_fluxes_rnd_right)}},
+               {{left_id, std::move(expected_operator_vars_rnd_left)},
+                {right_id, std::move(expected_operator_vars_rnd_right)}}}};
       // Large tolerances for the comparison to the analytic solution because
       // this regression test runs at very low resolution. Below is another
       // analytic-solution test at higher resolution. The hard-coded numbers
@@ -682,18 +698,18 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") {
       Approx analytic_solution_operator_approx =
           Approx::custom().epsilon(3.e-2).scale(M_PI * penalty_parameter *
                                                 square(4) / 0.5);
-      test_dg_operator<system, true>(
-          domain_creator, penalty_parameter, false,
-          Spectral::Quadrature::GaussLobatto, analytic_solution,
-          analytic_solution_aux_approx, analytic_solution_operator_approx,
-          {{{{left_id, std::move(vars_rnd_left)},
-             {midleft_id, std::move(vars_rnd_midleft)},
-             {midright_id, std::move(vars_rnd_midright)},
-             {right_id, std::move(vars_rnd_right)}},
-            {{left_id, std::move(expected_primal_fluxes_rnd_left)},
-             {right_id, std::move(expected_primal_fluxes_rnd_right)}},
-            {{left_id, std::move(expected_operator_vars_rnd_left)},
-             {right_id, std::move(expected_operator_vars_rnd_right)}}}});
+      // Numbers should be identical for the different formulations on affine
+      // coordinate maps
+      for (const auto dg_formulation :
+           make_array(::dg::Formulation::StrongInertial,
+                      ::dg::Formulation::StrongLogical,
+                      ::dg::Formulation::WeakInertial)) {
+        test_dg_operator<system, true>(
+            domain_creator, penalty_parameter, false,
+            Spectral::Quadrature::GaussLobatto, dg_formulation,
+            analytic_solution, analytic_solution_aux_approx,
+            analytic_solution_operator_approx, regression_test_data);
+      }
     }
     {
       INFO("Higher-resolution analytic-solution tests");
@@ -721,8 +737,9 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") {
                                         Spectral::Quadrature::Gauss))) {
         test_dg_operator<system, true>(
             domain_creator, penalty_parameter, use_massive_dg_operator,
-            quadrature, analytic_solution, analytic_solution_aux_approx,
-            analytic_solution_operator_approx, {}, true);
+            quadrature, ::dg::Formulation::StrongInertial, analytic_solution,
+            analytic_solution_aux_approx, analytic_solution_operator_approx, {},
+            true);
       }
     }
   }
@@ -787,6 +804,16 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") {
           expected_operator_vars_rnd)) =
           DataVector{164.82044058319110, 9.68580366789113,  -2.370110768729976,
                      91.310963425225239, 30.31342257245155, 56.03237239864831};
+      const std::vector<
+          std::tuple<std::unordered_map<ElementId<2>, Vars>,
+                     std::unordered_map<ElementId<2>, PrimalFluxes>,
+                     std::unordered_map<ElementId<2>, OperatorVars>>>
+          regression_test_data{
+              {{{northwest_id, std::move(vars_rnd_northwest)},
+                {southwest_id, std::move(vars_rnd_southwest)},
+                {northeast_id, std::move(vars_rnd_northeast)}},
+               {{northwest_id, std::move(expected_primal_fluxes_rnd)}},
+               {{northwest_id, std::move(expected_operator_vars_rnd)}}}};
       // Large tolerances for the comparison to the analytic solution because
       // this regression test runs at very low resolution. Below is another
       // analytic-solution test at higher resolution. The hard-coded numbers
@@ -796,15 +823,18 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") {
       Approx analytic_solution_operator_approx =
           Approx::custom().epsilon(7.e-1).scale(M_PI * penalty_parameter *
                                                 square(3));
-      test_dg_operator<system, true>(
-          domain_creator, penalty_parameter, false,
-          Spectral::Quadrature::GaussLobatto, analytic_solution,
-          analytic_solution_aux_approx, analytic_solution_operator_approx,
-          {{{{northwest_id, std::move(vars_rnd_northwest)},
-             {southwest_id, std::move(vars_rnd_southwest)},
-             {northeast_id, std::move(vars_rnd_northeast)}},
-            {{northwest_id, std::move(expected_primal_fluxes_rnd)}},
-            {{northwest_id, std::move(expected_operator_vars_rnd)}}}});
+      // Numbers should be identical for the different formulations on affine
+      // coordinate maps
+      for (const auto dg_formulation :
+           make_array(::dg::Formulation::StrongInertial,
+                      ::dg::Formulation::StrongLogical,
+                      ::dg::Formulation::WeakInertial)) {
+        test_dg_operator<system, true>(
+            domain_creator, penalty_parameter, false,
+            Spectral::Quadrature::GaussLobatto, dg_formulation,
+            analytic_solution, analytic_solution_aux_approx,
+            analytic_solution_operator_approx, regression_test_data);
+      }
     }
     {
       INFO("Higher-resolution analytic-solution tests");
@@ -830,8 +860,9 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") {
                                         Spectral::Quadrature::Gauss))) {
         test_dg_operator<system, true>(
             domain_creator, penalty_parameter, use_massive_dg_operator,
-            quadrature, analytic_solution, analytic_solution_aux_approx,
-            analytic_solution_operator_approx, {}, true);
+            quadrature, ::dg::Formulation::StrongInertial, analytic_solution,
+            analytic_solution_aux_approx, analytic_solution_operator_approx, {},
+            true);
       }
     }
   }
@@ -860,6 +891,8 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") {
       const ElementId<3> neighbor_id_eta{0, {{{1, 0}, {1, 1}, {1, 0}}}};
       const ElementId<3> neighbor_id_zeta{0, {{{1, 0}, {1, 0}, {1, 1}}}};
       using Vars = Variables<tmpl::list<Var<Poisson::Tags::Field>>>;
+      using PrimalFluxes = Variables<tmpl::list<Var<::Tags::Flux<
+          Poisson::Tags::Field, tmpl::size_t<3>, Frame::Inertial>>>>;
       using OperatorVars =
           Variables<tmpl::list<DgOperatorAppliedTo<Var<Poisson::Tags::Field>>>>;
       Vars vars_rnd_self{24};
@@ -913,6 +946,17 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") {
           23.036282532293448, -149.02371403160321, -142.12061361346408,
           314.95846885876898, 329.52673640830761,  46.245499085300622,
           70.51403818198222,  45.264231133848220,  90.061981328724073};
+      const std::vector<
+          std::tuple<std::unordered_map<ElementId<3>, Vars>,
+                     std::unordered_map<ElementId<3>, PrimalFluxes>,
+                     std::unordered_map<ElementId<3>, OperatorVars>>>
+          regression_test_data{
+              {{{self_id, std::move(vars_rnd_self)},
+                {neighbor_id_xi, std::move(vars_rnd_neighbor_xi)},
+                {neighbor_id_eta, std::move(vars_rnd_neighbor_eta)},
+                {neighbor_id_zeta, std::move(vars_rnd_neighbor_zeta)}},
+               {},
+               {{self_id, std::move(expected_operator_vars_rnd)}}}};
       // Large tolerances for the comparison to the analytic solution because
       // this regression test runs at very low resolution. Below is another
       // analytic-solution test at higher resolution. The hard-coded numbers
@@ -922,16 +966,18 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") {
       Approx analytic_solution_operator_approx =
           Approx::custom().epsilon(8.e-1).scale(M_PI * penalty_parameter *
                                                 square(4) / 2.);
-      test_dg_operator<system, true>(
-          domain_creator, penalty_parameter, false,
-          Spectral::Quadrature::GaussLobatto, analytic_solution,
-          analytic_solution_aux_approx, analytic_solution_operator_approx,
-          {{{{self_id, std::move(vars_rnd_self)},
-             {neighbor_id_xi, std::move(vars_rnd_neighbor_xi)},
-             {neighbor_id_eta, std::move(vars_rnd_neighbor_eta)},
-             {neighbor_id_zeta, std::move(vars_rnd_neighbor_zeta)}},
-            {},
-            {{self_id, std::move(expected_operator_vars_rnd)}}}});
+      // Numbers should be identical for the different formulations on affine
+      // coordinate maps
+      for (const auto dg_formulation :
+           make_array(::dg::Formulation::StrongInertial,
+                      ::dg::Formulation::StrongLogical,
+                      ::dg::Formulation::WeakInertial)) {
+        test_dg_operator<system, true>(
+            domain_creator, penalty_parameter, false,
+            Spectral::Quadrature::GaussLobatto, dg_formulation,
+            analytic_solution, analytic_solution_aux_approx,
+            analytic_solution_operator_approx, regression_test_data);
+      }
     }
     {
       INFO("Higher-resolution analytic-solution tests");
@@ -958,8 +1004,9 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") {
                                         Spectral::Quadrature::Gauss))) {
         test_dg_operator<system, true>(
             domain_creator, penalty_parameter, use_massive_dg_operator,
-            quadrature, analytic_solution, analytic_solution_aux_approx,
-            analytic_solution_operator_approx, {}, true);
+            quadrature, ::dg::Formulation::StrongInertial, analytic_solution,
+            analytic_solution_aux_approx, analytic_solution_operator_approx, {},
+            true);
       }
     }
   }
@@ -1044,6 +1091,18 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") {
           expected_operator_vars_rnd_right)) = DataVector{
           1.60618450837350557, 4.71949148249203443, 15.72408761869980509,
           4.06376669069456398, 5.88578668691303086, 15.11012286654655945};
+      const std::vector<
+          std::tuple<std::unordered_map<ElementId<2>, Vars>,
+                     std::unordered_map<ElementId<2>, PrimalFluxes>,
+                     std::unordered_map<ElementId<2>, OperatorVars>>>
+          regression_test_data{
+              {{{lowerleft_id, std::move(vars_rnd_lowerleft)},
+                {upperleft_id, std::move(vars_rnd_upperleft)},
+                {right_id, std::move(vars_rnd_right)}},
+               {{lowerleft_id, std::move(expected_primal_fluxes_rnd_lowerleft)},
+                {right_id, std::move(expected_primal_fluxes_rnd_right)}},
+               {{lowerleft_id, std::move(expected_operator_vars_rnd_lowerleft)},
+                {right_id, std::move(expected_operator_vars_rnd_right)}}}};
       // Large tolerances for the comparison to the analytic solution because
       // this regression test runs at very low resolution. Below is another
       // analytic-solution test at higher resolution. The hard-coded numbers
@@ -1053,17 +1112,16 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") {
       Approx analytic_solution_operator_approx =
           Approx::custom().epsilon(7.e-1).scale(M_PI * penalty_parameter *
                                                 square(3));
-      test_dg_operator<system, true>(
-          domain_creator, penalty_parameter, true,
-          Spectral::Quadrature::GaussLobatto, analytic_solution,
-          analytic_solution_aux_approx, analytic_solution_operator_approx,
-          {{{{lowerleft_id, std::move(vars_rnd_lowerleft)},
-             {upperleft_id, std::move(vars_rnd_upperleft)},
-             {right_id, std::move(vars_rnd_right)}},
-            {{lowerleft_id, std::move(expected_primal_fluxes_rnd_lowerleft)},
-             {right_id, std::move(expected_primal_fluxes_rnd_right)}},
-            {{lowerleft_id, std::move(expected_operator_vars_rnd_lowerleft)},
-             {right_id, std::move(expected_operator_vars_rnd_right)}}}});
+      // Numbers should be identical for the strong formulations
+      for (const auto dg_formulation :
+           make_array(::dg::Formulation::StrongInertial,
+                      ::dg::Formulation::StrongLogical)) {
+        test_dg_operator<system, true>(
+            domain_creator, penalty_parameter, true,
+            Spectral::Quadrature::GaussLobatto, dg_formulation,
+            analytic_solution, analytic_solution_aux_approx,
+            analytic_solution_operator_approx, regression_test_data);
+      }
     }
     {
       INFO("Higher-resolution analytic-solution tests");
@@ -1085,13 +1143,15 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") {
       Approx analytic_solution_operator_approx =
           Approx::custom().epsilon(1.e-11).scale(M_PI * penalty_parameter *
                                                  square(12));
-      for (const auto& [massive, quadrature] :
+      for (const auto& [massive, quadrature, dg_formulation] :
            cartesian_product(make_array(true, false),
                              make_array(Spectral::Quadrature::GaussLobatto,
-                                        Spectral::Quadrature::Gauss))) {
+                                        Spectral::Quadrature::Gauss),
+                             make_array(::dg::Formulation::StrongInertial,
+                                        ::dg::Formulation::WeakInertial))) {
         test_dg_operator<system, true>(
             domain_creator, penalty_parameter, massive, quadrature,
-            analytic_solution, analytic_solution_aux_approx,
+            dg_formulation, analytic_solution, analytic_solution_aux_approx,
             analytic_solution_operator_approx, {}, true);
       }
     }
@@ -1190,8 +1250,9 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") {
           Approx::custom().epsilon(0.3).scale(1.);
       test_dg_operator<system, true>(
           domain_creator, penalty_parameter, true,
-          Spectral::Quadrature::GaussLobatto, analytic_solution,
-          analytic_solution_aux_approx, analytic_solution_operator_approx,
+          Spectral::Quadrature::GaussLobatto, ::dg::Formulation::StrongInertial,
+          analytic_solution, analytic_solution_aux_approx,
+          analytic_solution_operator_approx,
           {{{{center_id, std::move(vars_rnd_center)},
              {wedge_id, std::move(vars_rnd_wedge)}},
             {},
@@ -1220,10 +1281,14 @@ SPECTRE_TEST_CASE("Unit.Elliptic.DG.Operator", "[Unit][Elliptic]") {
           Approx::custom().epsilon(1.e-4).scale(1.);
       Approx analytic_solution_operator_approx =
           Approx::custom().epsilon(1.e-4).scale(1.);
-      for (const auto quadrature :
-           {Spectral::Quadrature::GaussLobatto, Spectral::Quadrature::Gauss}) {
+      for (const auto& [quadrature, dg_formulation] :
+           cartesian_product(make_array(Spectral::Quadrature::GaussLobatto,
+                                        Spectral::Quadrature::Gauss),
+                             make_array(::dg::Formulation::StrongInertial,
+                                        ::dg::Formulation::StrongLogical,
+                                        ::dg::Formulation::WeakInertial))) {
         test_dg_operator<system, true>(
-            domain_creator, penalty_parameter, true, quadrature,
+            domain_creator, penalty_parameter, true, quadrature, dg_formulation,
             analytic_solution, analytic_solution_aux_approx,
             analytic_solution_operator_approx, {}, true);
       }
diff --git a/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_LargeOuterRadius.cpp b/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_LargeOuterRadius.cpp
new file mode 100644
index 000000000000..f8054d21a586
--- /dev/null
+++ b/tests/Unit/Elliptic/DiscontinuousGalerkin/Test_LargeOuterRadius.cpp
@@ -0,0 +1,184 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Framework/TestingFramework.hpp"
+
+#include <cstddef>
+
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Tensor/EagerMath/Determinant.hpp"
+#include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
+#include "Domain/CoordinateMaps/CoordinateMap.tpp"
+#include "Domain/CoordinateMaps/Wedge.hpp"
+#include "Domain/ElementMap.hpp"
+#include "Domain/FaceNormal.hpp"
+#include "Domain/InterfaceLogicalCoordinates.hpp"
+#include "Domain/Structure/ElementId.hpp"
+#include "Elliptic/DiscontinuousGalerkin/DgOperator.hpp"
+#include "Elliptic/Systems/Poisson/FirstOrderSystem.hpp"
+#include "Elliptic/Systems/Poisson/Tags.hpp"
+#include "NumericalAlgorithms/DiscontinuousGalerkin/ApplyMassMatrix.hpp"
+#include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
+#include "NumericalAlgorithms/DiscontinuousGalerkin/MetricIdentityJacobian.hpp"
+#include "NumericalAlgorithms/DiscontinuousGalerkin/NormalDotFlux.hpp"
+#include "NumericalAlgorithms/LinearOperators/Divergence.hpp"
+#include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
+#include "NumericalAlgorithms/LinearOperators/WeakDivergence.hpp"
+#include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp"
+#include "NumericalAlgorithms/Spectral/Mesh.hpp"
+#include "PointwiseFunctions/AnalyticSolutions/Poisson/Lorentzian.hpp"
+#include "Utilities/TMPL.hpp"
+
+SPECTRE_TEST_CASE("Unit.Elliptic.DG.LargeOuterRadius", "[Unit][Elliptic]") {
+  // This test works for the StrongLogical and WeakInertial formulations, but
+  // breaks for the StrongInertial formulation.
+  const ::dg::Formulation formulation = ::dg::Formulation::StrongLogical;
+  // The test also breaks if the constant is set to 1. _and_ the numeric first
+  // derivative is used. It works fine if the analytic first derivative is used
+  // or if the constant is set to 0.
+  const Poisson::Solutions::Lorentzian<3> solution{/* plus_constant */ 1.};
+  const bool use_numeric_first_deriv = false;
+  // Set up the grid
+  using Wedge3D = domain::CoordinateMaps::Wedge<3>;
+  const double inner_radius = 1e2;
+  CAPTURE(inner_radius);
+  const double outer_radius = 1e9;
+  CAPTURE(outer_radius);
+  auto block_map =
+      domain::make_coordinate_map_base<Frame::BlockLogical, Frame::Inertial>(
+          Wedge3D{inner_radius, outer_radius, 1., 1.,
+                  OrientationMap<3>::create_aligned(), true,
+                  Wedge3D::WedgeHalves::Both,
+                  domain::CoordinateMaps::Distribution::Inverse});
+  const ElementId<3> element_id{0, {{{2, 0}, {2, 0}, {0, 0}}}};
+  CAPTURE(element_id);
+  const Mesh<3> mesh{12, Spectral::Basis::Legendre,
+                     Spectral::Quadrature::GaussLobatto};
+  CAPTURE(mesh);
+
+  // Quantify grid stretching by looking at the midpoint of the element
+  const ElementMap<3, Frame::Inertial> element_map{element_id,
+                                                   std::move(block_map)};
+  const auto midpoint_radius = magnitude(
+      element_map(tnsr::I<double, 3, Frame::ElementLogical>{{{0., 0., 0.}}}));
+  CAPTURE(midpoint_radius);
+  Approx custom_approx = Approx::custom().epsilon(1.0e-6).scale(1.0);
+  const size_t num_points = mesh.number_of_grid_points();
+  const auto logical_coords = logical_coordinates(mesh);
+  const auto inertial_coords = element_map(logical_coords);
+  const auto inv_jacobian = element_map.inv_jacobian(logical_coords);
+  const auto det_jacobian = determinant(element_map.jacobian(logical_coords));
+
+  // Take first derivative
+  using Field = Poisson::Tags::Field;
+  using FieldDeriv =
+      ::Tags::deriv<Poisson::Tags::Field, tmpl::size_t<3>, Frame::Inertial>;
+  using Flux =
+      ::Tags::Flux<Poisson::Tags::Field, tmpl::size_t<3>, Frame::Inertial>;
+  using FixedSource = ::Tags::FixedSource<Field>;
+  const auto vars = solution.variables(
+      inertial_coords, tmpl::list<Field, FieldDeriv, FixedSource>{});
+  const auto& u = get<Field>(vars);
+  const auto& du_analytic = get<FieldDeriv>(vars);
+  const auto& f = get<FixedSource>(vars);
+  const auto du_numeric = partial_derivative(u, mesh, inv_jacobian);
+  CHECK_ITERABLE_CUSTOM_APPROX(du_numeric, du_analytic, custom_approx);
+
+  // Take second derivative
+  CAPTURE(formulation);
+  const bool massive = true;
+  CAPTURE(massive);
+  tnsr::I<DataVector, 3> flux{num_points};
+  Poisson::flat_cartesian_fluxes(
+      make_not_null(&flux), use_numeric_first_deriv ? du_numeric : du_analytic);
+  if (formulation == ::dg::Formulation::StrongInertial) {
+    Scalar<DataVector> ddu{num_points, 0.};
+    divergence(make_not_null(&ddu), flux, mesh, inv_jacobian);
+    DataVector lhs = -get(ddu);
+    DataVector rhs = get(f);
+    // Massive scheme multiplies by Jacobian determinant
+    if (massive) {
+      lhs *= get(det_jacobian);
+      rhs *= get(det_jacobian);
+    }
+    CHECK_ITERABLE_CUSTOM_APPROX(lhs, rhs, custom_approx);
+  } else if (formulation == ::dg::Formulation::StrongLogical) {
+    Scalar<DataVector> ddu{num_points, 0.};
+    tnsr::I<DataVector, 3, Frame::ElementLogical> logical_flux =
+        transform::first_index_to_different_frame(flux, inv_jacobian);
+    for (size_t d = 0; d < 3; ++d) {
+      logical_flux.get(d) *= get(det_jacobian);
+    }
+    logical_divergence(make_not_null(&ddu), logical_flux, mesh);
+    DataVector lhs = -get(ddu);
+    DataVector rhs = get(f);
+    // Massive scheme multiplies by Jacobian determinant
+    if (massive) {
+      ::dg::apply_mass_matrix(make_not_null(&lhs), mesh);
+      ::dg::apply_mass_matrix(make_not_null(&rhs), mesh);
+      rhs *= get(det_jacobian);
+    } else {
+      lhs /= get(det_jacobian);
+    }
+    CHECK_ITERABLE_CUSTOM_APPROX(lhs, rhs, custom_approx);
+  } else {
+    auto det_times_inv_jacobian = inv_jacobian;
+    for (auto& component : det_times_inv_jacobian) {
+      component *= get(det_jacobian);
+    }
+    Variables<tmpl::list<Flux>> fluxes{mesh.number_of_grid_points()};
+    get<Flux>(fluxes) = flux;
+    Variables<tmpl::list<Field>> lhs{mesh.number_of_grid_points()};
+    weak_divergence(make_not_null(&lhs), fluxes, mesh, det_times_inv_jacobian);
+    DataVector rhs = get(f);
+    if (massive) {
+      ::dg::apply_mass_matrix(make_not_null(&lhs), mesh);
+      ::dg::apply_mass_matrix(make_not_null(&rhs), mesh);
+      rhs *= get(det_jacobian);
+    } else {
+      lhs /= get(det_jacobian);
+    }
+    // Add boundary corrections
+    for (const auto direction : Direction<3>::all_directions()) {
+      CAPTURE(direction);
+      // Compute face geometry
+      const auto face_mesh = mesh.slice_away(direction.dimension());
+      const size_t face_num_points = face_mesh.number_of_grid_points();
+      const auto face_logical_coords =
+          interface_logical_coordinates(face_mesh, direction);
+      const auto face_inertial_coords = element_map(face_logical_coords);
+      auto face_normal =
+          unnormalized_face_normal(face_mesh, element_map, direction);
+      const auto face_normal_magnitude = magnitude(face_normal);
+      for (size_t d = 0; d < 3; ++d) {
+        face_normal.get(d) /= get(face_normal_magnitude);
+      }
+      auto face_jacobian = face_normal_magnitude;
+      const auto det_jacobian_face =
+          determinant(element_map.jacobian(face_logical_coords));
+      get(face_jacobian) *= get(det_jacobian_face);
+      // Compute boundary correction
+      const auto vars_on_face =
+          solution.variables(face_inertial_coords, tmpl::list<FieldDeriv>{});
+      const auto& du_on_face = get<FieldDeriv>(vars_on_face);
+      Variables<tmpl::list<Flux>> fluxes_on_face{face_num_points};
+      Poisson::flat_cartesian_fluxes(make_not_null(&get<Flux>(fluxes_on_face)),
+                                     du_on_face);
+      auto boundary_correction =
+          normal_dot_flux<tmpl::list<Field>>(face_normal, fluxes_on_face);
+      ASSERT(massive, "Only massive scheme supported here");
+      ::dg::apply_mass_matrix(make_not_null(&boundary_correction), face_mesh);
+      boundary_correction *= get(face_jacobian);
+      boundary_correction *= -1.;
+      if (mesh.quadrature(0) == Spectral::Quadrature::GaussLobatto) {
+        add_slice_to_data(make_not_null(&lhs), boundary_correction,
+                          mesh.extents(), direction.dimension(),
+                          index_to_slice_at(mesh.extents(), direction));
+      } else {
+        ::dg::lift_boundary_terms_gauss_points(
+            make_not_null(&lhs), boundary_correction, mesh, direction);
+      }
+    }
+    CHECK_ITERABLE_CUSTOM_APPROX(get(get<Field>(lhs)), rhs, custom_approx);
+  }
+}
diff --git a/tests/Unit/NumericalAlgorithms/DiscontinuousGalerkin/Tags/Test_Formulation.cpp b/tests/Unit/NumericalAlgorithms/DiscontinuousGalerkin/Tags/Test_Formulation.cpp
index 341eea796a44..257ca9bd554d 100644
--- a/tests/Unit/NumericalAlgorithms/DiscontinuousGalerkin/Tags/Test_Formulation.cpp
+++ b/tests/Unit/NumericalAlgorithms/DiscontinuousGalerkin/Tags/Test_Formulation.cpp
@@ -17,4 +17,6 @@ SPECTRE_TEST_CASE("Unit.DiscontinuousGalerkin.Tags.Formulation",
             "StrongInertial") == dg::Formulation::StrongInertial);
   CHECK(TestHelpers::test_option_tag<dg::OptionTags::Formulation>(
             "WeakInertial") == dg::Formulation::WeakInertial);
+  CHECK(TestHelpers::test_option_tag<dg::OptionTags::Formulation>(
+            "StrongLogical") == dg::Formulation::StrongLogical);
 }
diff --git a/tests/Unit/NumericalAlgorithms/DiscontinuousGalerkin/Test_Formulation.cpp b/tests/Unit/NumericalAlgorithms/DiscontinuousGalerkin/Test_Formulation.cpp
index 3a442caa4761..13a1940d7ac5 100644
--- a/tests/Unit/NumericalAlgorithms/DiscontinuousGalerkin/Test_Formulation.cpp
+++ b/tests/Unit/NumericalAlgorithms/DiscontinuousGalerkin/Test_Formulation.cpp
@@ -12,4 +12,5 @@ SPECTRE_TEST_CASE("Unit.DiscontinuousGalerkin.Formulation",
                   "[Unit][NumericalAlgorithms]") {
   CHECK(get_output(dg::Formulation::StrongInertial) == "StrongInertial");
   CHECK(get_output(dg::Formulation::WeakInertial) == "WeakInertial");
+  CHECK(get_output(dg::Formulation::StrongInertial) == "StrongInertial");
 }
diff --git a/tests/Unit/NumericalAlgorithms/LinearOperators/Test_Divergence.cpp b/tests/Unit/NumericalAlgorithms/LinearOperators/Test_Divergence.cpp
index f0bb4165a833..3b571fb3402b 100644
--- a/tests/Unit/NumericalAlgorithms/LinearOperators/Test_Divergence.cpp
+++ b/tests/Unit/NumericalAlgorithms/LinearOperators/Test_Divergence.cpp
@@ -17,8 +17,10 @@
 #include "DataStructures/DataBox/Prefixes.hpp"
 #include "DataStructures/DataBox/Tag.hpp"
 #include "DataStructures/DataVector.hpp"
+#include "DataStructures/Tensor/EagerMath/Determinant.hpp"
 #include "DataStructures/Tensor/Tensor.hpp"
 #include "DataStructures/Variables.hpp"
+#include "DataStructures/Variables/FrameTransform.hpp"
 #include "DataStructures/VariablesTag.hpp"
 #include "Domain/CoordinateMaps/Affine.hpp"
 #include "Domain/CoordinateMaps/CoordinateMap.hpp"
@@ -143,6 +145,7 @@ void test_divergence_impl(
   const auto xi = logical_coordinates(mesh);
   const auto x = coordinate_map(xi);
   const auto inv_jacobian = coordinate_map.inv_jacobian(xi);
+  const auto det_jacobian = determinant(coordinate_map.jacobian(xi));
   MathFunctions::TensorProduct<Dim> f(1.0, std::move(functions));
   using flux_tags = two_fluxes<DataType, Dim, Frame>;
   Variables<flux_tags> fluxes(num_grid_points);
@@ -166,6 +169,15 @@ void test_divergence_impl(
   const auto& expected =
       get<Tags::div<Flux1<DataType, Dim, Frame>>>(div_fluxes);
   CHECK_ITERABLE_CUSTOM_APPROX(expected, div_vector, local_approx);
+
+  // Test logical divergence
+  auto logical_fluxes =
+      transform::first_index_to_different_frame(fluxes, inv_jacobian);
+  Variables<db::wrap_tags_in<Tags::div, flux_tags>> logical_div_fluxes(
+      num_grid_points);
+  logical_divergence(make_not_null(&logical_div_fluxes), logical_fluxes, mesh);
+  CHECK_VARIABLES_CUSTOM_APPROX(logical_div_fluxes, expected_div_fluxes,
+                                local_approx);
 }
 
 template <typename DataType>
@@ -314,7 +326,7 @@ void test_divergence_compute() {
 }
 }  // namespace
 
-// [[Timeout, 10]]
+// [[Timeout, 20]]
 SPECTRE_TEST_CASE("Unit.Numerical.LinearOperators.Divergence",
                   "[NumericalAlgorithms][LinearOperators][Unit]") {
   test_divergence<DataVector>();