Merge pull request #6245 from nilsvu/dg_paper_tests

Add logical_divergence function, elliptic StrongLogical DG formulation

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.

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) {

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;
 

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;
 

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>;
 

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

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(

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;
 

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.");
 }

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

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

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

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));
+}

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}")