Merge pull request #5679 from nilsvu/elliptic_weak_form

Support weak form in elliptic DG
sxs-collaboration · Dec 17, 2023 · e8c48be · e8c48be
2 parents 6585ad3 + 2a1ea7b
commit e8c48be
Show file tree

Hide file tree

Showing 22 changed files with 276 additions and 59 deletions.
diff --git a/src/Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp b/src/Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp
@@ -356,7 +356,7 @@ struct ReceiveMortarDataAndApplyOperator<
  public:
   using const_global_cache_tags =
       tmpl::list<elliptic::dg::Tags::PenaltyParameter,
-                 elliptic::dg::Tags::Massive>;
+                 elliptic::dg::Tags::Massive, elliptic::dg::Tags::Formulation>;
   using inbox_tags = tmpl::list<mortar_data_inbox_tag>;
 
   template <typename DbTags, typename... InboxTags, typename Metavariables,
@@ -403,6 +403,11 @@ struct ReceiveMortarDataAndApplyOperator<
                                               Frame::Inertial>>(box),
         db::get<domain::Tags::DetInvJacobian<Frame::ElementLogical,
                                              Frame::Inertial>>(box),
+        db::get<
+            domain::Tags::DetJacobian<Frame::ElementLogical, Frame::Inertial>>(
+            box),
+        db::get<domain::Tags::DetTimesInvJacobian<Dim, Frame::ElementLogical,
+                                                  Frame::Inertial>>(box),
         db::get<domain::Tags::Faces<
             Dim, domain::Tags::UnnormalizedFaceNormalMagnitude<Dim>>>(box),
         db::get<domain::Tags::Faces<
@@ -417,7 +422,8 @@ struct ReceiveMortarDataAndApplyOperator<
                                     Frame::ElementLogical, Frame::Inertial>,
                                 Dim>>(box),
         db::get<elliptic::dg::Tags::PenaltyParameter>(box),
-        db::get<elliptic::dg::Tags::Massive>(box), temporal_id,
+        db::get<elliptic::dg::Tags::Massive>(box),
+        db::get<elliptic::dg::Tags::Formulation>(box), temporal_id,
         std::forward_as_tuple(db::get<SourcesArgsTags>(box)...));
 
     return {Parallel::AlgorithmExecution::Continue, std::nullopt};
@@ -526,7 +532,7 @@ struct ImposeInhomogeneousBoundaryConditionsOnSource<
  public:
   using const_global_cache_tags =
       tmpl::list<elliptic::dg::Tags::PenaltyParameter,
-                 elliptic::dg::Tags::Massive,
+                 elliptic::dg::Tags::Massive, elliptic::dg::Tags::Formulation,
                  domain::Tags::ExternalBoundaryConditions<Dim>>;
 
   template <typename DbTags, typename... InboxTags, typename Metavariables,
@@ -596,6 +602,11 @@ struct ImposeInhomogeneousBoundaryConditionsOnSource<
                                               Frame::Inertial>>(box),
         db::get<domain::Tags::DetInvJacobian<Frame::ElementLogical,
                                              Frame::Inertial>>(box),
+        db::get<
+            domain::Tags::DetJacobian<Frame::ElementLogical, Frame::Inertial>>(
+            box),
+        db::get<domain::Tags::DetTimesInvJacobian<Dim, Frame::ElementLogical,
+                                                  Frame::Inertial>>(box),
         db::get<domain::Tags::Faces<Dim, domain::Tags::FaceNormal<Dim>>>(box),
         db::get<domain::Tags::Faces<Dim, domain::Tags::FaceNormalVector<Dim>>>(
             box),
@@ -610,7 +621,8 @@ struct ImposeInhomogeneousBoundaryConditionsOnSource<
         db::get<::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>>(box),
         db::get<::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>>(box),
         db::get<elliptic::dg::Tags::PenaltyParameter>(box),
-        db::get<elliptic::dg::Tags::Massive>(box), apply_boundary_condition,
+        db::get<elliptic::dg::Tags::Massive>(box),
+        db::get<elliptic::dg::Tags::Formulation>(box), apply_boundary_condition,
         std::forward_as_tuple(db::get<FluxesArgsTags>(box)...),
         std::forward_as_tuple(db::get<SourcesArgsTags>(box)...),
         fluxes_args_on_faces);

diff --git a/src/Elliptic/DiscontinuousGalerkin/DgOperator.hpp b/src/Elliptic/DiscontinuousGalerkin/DgOperator.hpp
@@ -24,6 +24,7 @@
 #include "Elliptic/Protocols/FirstOrderSystem.hpp"
 #include "Elliptic/Systems/GetSourcesComputer.hpp"
 #include "NumericalAlgorithms/DiscontinuousGalerkin/ApplyMassMatrix.hpp"
+#include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
 #include "NumericalAlgorithms/DiscontinuousGalerkin/LiftFlux.hpp"
 #include "NumericalAlgorithms/DiscontinuousGalerkin/LiftFromBoundary.hpp"
 #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
@@ -595,6 +596,9 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
       const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
                             Frame::Inertial>& inv_jacobian,
       const Scalar<DataVector>& det_inv_jacobian,
+      const Scalar<DataVector>& det_jacobian,
+      const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
+                            Frame::Inertial>& det_times_inv_jacobian,
       const DirectionMap<Dim, Scalar<DataVector>>& face_normal_magnitudes,
       const DirectionMap<Dim, Scalar<DataVector>>& face_jacobians,
       const DirectionMap<Dim,
@@ -605,7 +609,7 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
       const ::dg::MortarMap<Dim, ::dg::MortarSize<Dim - 1>>& all_mortar_sizes,
       const ::dg::MortarMap<Dim, Scalar<DataVector>>& mortar_jacobians,
       const double penalty_parameter, const bool massive,
-      const TemporalId& /*temporal_id*/,
+      const ::dg::Formulation formulation, const TemporalId& /*temporal_id*/,
       const std::tuple<SourcesArgs...>& sources_args,
       const DataIsZero& data_is_zero = NoDataIsZero{},
       const DirectionsPredicate& directions_predicate = AllDirections{}) {
@@ -661,26 +665,48 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
     if (local_data_is_zero) {
       operator_applied_to_vars->initialize(num_points, 0.);
     } else {
-      divergence(operator_applied_to_vars, primal_fluxes, mesh, inv_jacobian);
-      // This is the sign flip that makes the operator _minus_ the Laplacian for
-      // a Poisson system
-      *operator_applied_to_vars *= -1.;
+      // "Massive" operators retain the factors from the volume integral:
+      //   \int_volume div(F) \phi_p = w_p det(J)_p div(F)_p
+      // Here, `w` are the quadrature weights (the diagonal logical mass matrix
+      // with mass-lumping) and det(J) is the Jacobian determinant. The
+      // quantities are evaluated at the grid point `p`.
+      if (formulation == ::dg::Formulation::StrongInertial) {
+        // Compute strong divergence:
+        //   div(F) = (J^\hat{i}_i)_p \sum_q (D_\hat{i})_pq (F^i)_q.
+        divergence(operator_applied_to_vars, primal_fluxes, mesh,
+                   massive ? det_times_inv_jacobian : inv_jacobian);
+        // This is the sign flip that makes the operator _minus_ the Laplacian
+        // for a Poisson system
+        *operator_applied_to_vars *= -1.;
+      } else {
+        // 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
+        weak_divergence(operator_applied_to_vars, primal_fluxes, mesh,
+                        det_times_inv_jacobian);
+        if (not massive) {
+          *operator_applied_to_vars *= get(det_inv_jacobian);
+        }
+      }
       if constexpr (not std::is_same_v<SourcesComputer, void>) {
+        Variables<tmpl::list<OperatorTags...>> sources{num_points, 0.};
         std::apply(
-            [&operator_applied_to_vars, &primal_vars,
+            [&sources, &primal_vars,
              &primal_fluxes](const auto&... expanded_sources_args) {
               SourcesComputer::apply(
-                  make_not_null(
-                      &get<OperatorTags>(*operator_applied_to_vars))...,
+                  make_not_null(&get<OperatorTags>(sources))...,
                   expanded_sources_args..., get<PrimalVars>(primal_vars)...,
                   get<PrimalFluxesVars>(primal_fluxes)...);
             },
             sources_args);
+        if (massive) {
+          sources *= get(det_jacobian);
+        }
+        *operator_applied_to_vars += sources;
       }
-      if (massive) {
-        *operator_applied_to_vars /= get(det_inv_jacobian);
-        ::dg::apply_mass_matrix(operator_applied_to_vars, mesh);
-      }
+    }
+    if (massive) {
+      ::dg::apply_mass_matrix(operator_applied_to_vars, mesh);
     }
 
     // Add boundary corrections
@@ -799,18 +825,20 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
           get<Tags::NormalDotFluxForJump<PrimalMortarVars>>(
               remote_data.field_data)));
       primal_boundary_corrections_on_mortar *= penalty;
-      const auto add_remote_avg_contribution = [](auto& lhs, const auto& rhs) {
+      const auto add_avg_contribution = [](auto& lhs, const auto& rhs,
+                                           const double prefactor) {
         for (size_t i = 0; i < lhs.size(); ++i) {
-          lhs[i] += 0.5 * rhs[i];
+          lhs[i] += prefactor * rhs[i];
         }
       };
-      EXPAND_PACK_LEFT_TO_RIGHT(add_remote_avg_contribution(
+      EXPAND_PACK_LEFT_TO_RIGHT(add_avg_contribution(
           get<PrimalMortarVars>(primal_boundary_corrections_on_mortar),
-          get<::Tags::NormalDotFlux<PrimalMortarVars>>(local_data.field_data)));
-      EXPAND_PACK_LEFT_TO_RIGHT(add_remote_avg_contribution(
+          get<::Tags::NormalDotFlux<PrimalMortarVars>>(local_data.field_data),
+          formulation == ::dg::Formulation::StrongInertial ? 0.5 : -0.5));
+      EXPAND_PACK_LEFT_TO_RIGHT(add_avg_contribution(
           get<PrimalMortarVars>(primal_boundary_corrections_on_mortar),
-          get<::Tags::NormalDotFlux<PrimalMortarVars>>(
-              remote_data.field_data)));
+          get<::Tags::NormalDotFlux<PrimalMortarVars>>(remote_data.field_data),
+          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.
@@ -910,6 +938,9 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
       const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
                             Frame::Inertial>& inv_jacobian,
       const Scalar<DataVector>& det_inv_jacobian,
+      const Scalar<DataVector>& det_jacobian,
+      const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
+                            Frame::Inertial>& det_times_inv_jacobian,
       const DirectionMap<Dim, tnsr::i<DataVector, Dim>>& face_normals,
       const DirectionMap<Dim, tnsr::I<DataVector, Dim>>& face_normal_vectors,
       const DirectionMap<Dim, Scalar<DataVector>>& face_normal_magnitudes,
@@ -921,6 +952,7 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
       const ::dg::MortarMap<Dim, Mesh<Dim - 1>>& all_mortar_meshes,
       const ::dg::MortarMap<Dim, ::dg::MortarSize<Dim - 1>>& all_mortar_sizes,
       const double penalty_parameter, const bool massive,
+      const ::dg::Formulation formulation,
       const ApplyBoundaryCondition& apply_boundary_condition,
       const std::tuple<FluxesArgs...>& fluxes_args,
       const std::tuple<SourcesArgs...>& sources_args,
@@ -958,13 +990,13 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
         face_normal_vectors, face_normal_magnitudes, all_mortar_meshes,
         all_mortar_sizes, temporal_id, apply_boundary_condition, fluxes_args,
         fluxes_args_on_faces);
-    apply_operator<true>(make_not_null(&operator_applied_to_zero_vars),
-                         make_not_null(&all_mortar_data), zero_primal_vars,
-                         primal_fluxes_buffer, element, mesh, inv_jacobian,
-                         det_inv_jacobian, face_normal_magnitudes,
-                         face_jacobians, face_jacobian_times_inv_jacobians,
-                         all_mortar_meshes, all_mortar_sizes, {},
-                         penalty_parameter, massive, temporal_id, sources_args);
+    apply_operator<true>(
+        make_not_null(&operator_applied_to_zero_vars),
+        make_not_null(&all_mortar_data), zero_primal_vars, primal_fluxes_buffer,
+        element, mesh, inv_jacobian, det_inv_jacobian, det_jacobian,
+        det_times_inv_jacobian, face_normal_magnitudes, face_jacobians,
+        face_jacobian_times_inv_jacobians, all_mortar_meshes, all_mortar_sizes,
+        {}, penalty_parameter, massive, formulation, temporal_id, sources_args);
     // Impose the nonlinear (constant) boundary contribution as fixed sources on
     // the RHS of the equations
     *fixed_sources -= operator_applied_to_zero_vars;

diff --git a/src/Elliptic/DiscontinuousGalerkin/Initialization.cpp b/src/Elliptic/DiscontinuousGalerkin/Initialization.cpp
@@ -49,6 +49,10 @@ void InitializeGeometry<Dim>::operator()(
                                         Frame::Inertial>*>
         inv_jacobian,
     const gsl::not_null<Scalar<DataVector>*> det_inv_jacobian,
+    const gsl::not_null<Scalar<DataVector>*> det_jacobian,
+    const gsl::not_null<InverseJacobian<DataVector, Dim, Frame::ElementLogical,
+                                        Frame::Inertial>*>
+        det_times_inv_jacobian,
     const std::vector<std::array<size_t, Dim>>& initial_extents,
     const std::vector<std::array<size_t, Dim>>& initial_refinement,
     const Domain<Dim>& domain,
@@ -74,9 +78,19 @@ void InitializeGeometry<Dim>::operator()(
   *inertial_coords =
       element_map->operator()(*logical_coords, 0., functions_of_time);
   // Jacobian
+  // Note: we can try to use `::dg::metric_identity_jacobian_quantities` here.
+  // When I tried (NV, Dec 2023) the DG residuals diverged on a sphere domain
+  // with a large outer boundary (1e9). This was possibly because no
+  // metric-identity Jacobians were used on faces, though I also tried slicing
+  // the metric-identity Jacobian to the faces and that didn't help.
   *inv_jacobian =
       element_map->inv_jacobian(*logical_coords, 0., functions_of_time);
   *det_inv_jacobian = determinant(*inv_jacobian);
+  get(*det_jacobian) = 1. / get(*det_inv_jacobian);
+  *det_times_inv_jacobian = *inv_jacobian;
+  for (auto& component : *det_times_inv_jacobian) {
+    component *= get(*det_jacobian);
+  }
 }
 
 namespace detail {

diff --git a/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp b/src/Elliptic/DiscontinuousGalerkin/Initialization.hpp
@@ -72,7 +72,10 @@ struct InitializeGeometry {
       domain::Tags::Coordinates<Dim, Frame::Inertial>,
       domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
                                     Frame::Inertial>,
-      domain::Tags::DetInvJacobian<Frame::ElementLogical, Frame::Inertial>>;
+      domain::Tags::DetInvJacobian<Frame::ElementLogical, Frame::Inertial>,
+      domain::Tags::DetJacobian<Frame::ElementLogical, Frame::Inertial>,
+      domain::Tags::DetTimesInvJacobian<Dim, Frame::ElementLogical,
+                                        Frame::Inertial>>;
   using argument_tags =
       tmpl::list<domain::Tags::InitialExtents<Dim>,
                  domain::Tags::InitialRefinementLevels<Dim>,
@@ -88,6 +91,10 @@ struct InitializeGeometry {
                                     Frame::Inertial>*>
           inv_jacobian,
       gsl::not_null<Scalar<DataVector>*> det_inv_jacobian,
+      gsl::not_null<Scalar<DataVector>*> det_jacobian,
+      gsl::not_null<InverseJacobian<DataVector, Dim, Frame::ElementLogical,
+                                    Frame::Inertial>*>
+          det_times_inv_jacobian,
       const std::vector<std::array<size_t, Dim>>& initial_extents,
       const std::vector<std::array<size_t, Dim>>& initial_refinement,
       const Domain<Dim>& domain,

diff --git a/src/Elliptic/DiscontinuousGalerkin/SubdomainOperator/SubdomainOperator.hpp b/src/Elliptic/DiscontinuousGalerkin/SubdomainOperator/SubdomainOperator.hpp
@@ -158,6 +158,9 @@ struct SubdomainOperator
       domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
                                     Frame::Inertial>,
       domain::Tags::DetInvJacobian<Frame::ElementLogical, Frame::Inertial>,
+      domain::Tags::DetJacobian<Frame::ElementLogical, Frame::Inertial>,
+      domain::Tags::DetTimesInvJacobian<Dim, Frame::ElementLogical,
+                                        Frame::Inertial>,
       domain::Tags::Faces<Dim,
                           domain::Tags::UnnormalizedFaceNormalMagnitude<Dim>>,
       domain::Tags::Faces<Dim, domain::Tags::DetSurfaceJacobian<
@@ -170,7 +173,8 @@ struct SubdomainOperator
       ::Tags::Mortars<domain::Tags::DetSurfaceJacobian<Frame::ElementLogical,
                                                        Frame::Inertial>,
                       Dim>,
-      elliptic::dg::Tags::PenaltyParameter, elliptic::dg::Tags::Massive>;
+      elliptic::dg::Tags::PenaltyParameter, elliptic::dg::Tags::Massive,
+      elliptic::dg::Tags::Formulation>;
   using fluxes_args_tags = typename System::fluxes_computer::argument_tags;
   using sources_args_tags =
       elliptic::get_sources_argument_tags<System, linearized>;
@@ -181,9 +185,9 @@ struct SubdomainOperator
 
   // These tags can be taken directly from the central element's DataBox, even
   // when evaluating neighbors
-  using args_tags_from_center = tmpl::remove_duplicates<
-      tmpl::push_back<ArgsTagsFromCenter, elliptic::dg::Tags::PenaltyParameter,
-                      elliptic::dg::Tags::Massive>>;
+  using args_tags_from_center = tmpl::remove_duplicates<tmpl::push_back<
+      ArgsTagsFromCenter, elliptic::dg::Tags::PenaltyParameter,
+      elliptic::dg::Tags::Massive, elliptic::dg::Tags::Formulation>>;
 
   // Data on neighbors is stored in the central element's DataBox in
   // `LinearSolver::Schwarz::Tags::Overlaps` maps, so we wrap the argument tags

diff --git a/src/Elliptic/DiscontinuousGalerkin/Tags.hpp b/src/Elliptic/DiscontinuousGalerkin/Tags.hpp
@@ -4,6 +4,7 @@
 #pragma once
 
 #include "DataStructures/DataBox/Tag.hpp"
+#include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
 #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
 #include "Options/String.hpp"
 #include "Utilities/ErrorHandling/Error.hpp"
@@ -52,6 +53,13 @@ struct Quadrature {
   using group = DiscontinuousGalerkin;
 };
 
+struct Formulation {
+  using type = ::dg::Formulation;
+  static constexpr Options::String help =
+      "The DG formulation to use (strong or weak).";
+  using group = DiscontinuousGalerkin;
+};
+
 }  // namespace OptionTags
 
 /// DataBox tags related to elliptic discontinuous Galerkin schemes
@@ -95,5 +103,13 @@ struct Quadrature : db::SimpleTag {
   }
 };
 
+/// The DG formulation to use (strong or weak)
+struct Formulation : db::SimpleTag {
+  using type = ::dg::Formulation;
+  static constexpr bool pass_metavariables = false;
+  using option_tags = tmpl::list<OptionTags::Formulation>;
+  static type create_from_options(const type value) { return value; }
+};
+
 }  // namespace Tags
 }  // namespace elliptic::dg