Support weak form in elliptic DG

sxs-collaboration · Dec 17, 2023 · 2a1ea7b · 2a1ea7b
1 parent f100dd8
commit 2a1ea7b
Show file tree

Hide file tree

Showing 11 changed files with 189 additions and 36 deletions.
diff --git a/src/Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp b/src/Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp
@@ -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<
@@ -597,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),

diff --git a/src/Elliptic/DiscontinuousGalerkin/DgOperator.hpp b/src/Elliptic/DiscontinuousGalerkin/DgOperator.hpp
@@ -596,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,
@@ -662,28 +665,48 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
     if (local_data_is_zero) {
       operator_applied_to_vars->initialize(num_points, 0.);
     } else {
-      ASSERT(formulation == ::dg::Formulation::StrongInertial,
-             "Only the strong formulation is currently implemented.");
-      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
@@ -802,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.
@@ -913,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,
@@ -965,10 +993,10 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
     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, formulation,
-        temporal_id, sources_args);
+        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/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<

diff --git a/src/Elliptic/DiscontinuousGalerkin/Tags.hpp b/src/Elliptic/DiscontinuousGalerkin/Tags.hpp
@@ -108,12 +108,7 @@ 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) {
-    if (not(value == ::dg::Formulation::StrongInertial)) {
-      ERROR_NO_TRACE("Only the strong formulation is currently supported.");
-    }
-    return value;
-  }
+  static type create_from_options(const type value) { return value; }
 };
 
 }  // namespace Tags

diff --git a/support/Pipelines/Bbh/InitialData.yaml b/support/Pipelines/Bbh/InitialData.yaml
@@ -85,7 +85,7 @@ Discretization:
     PenaltyParameter: 1.
     Massive: True
     Quadrature: GaussLobatto
-    Formulation: StrongInertial
+    Formulation: WeakInertial
 
 Observers:
   VolumeFileName: "BbhVolume"

diff --git a/tests/InputFiles/Poisson/Lorentzian.yaml b/tests/InputFiles/Poisson/Lorentzian.yaml
@@ -0,0 +1,117 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+Executable: SolvePoisson3D
+Testing:
+  Check: parse;execute_check_output
+  Timeout: 40
+  Priority: High
+ExpectedOutput:
+  - LorentzianReductions.h5
+  - LorentzianVolume0.h5
+OutputFileChecks:
+  - Label: Discretization error
+    Subfile: /ErrorNorms.dat
+    FileGlob: LorentzianReductions.h5
+    SkipColumns: [0, 1, 2]
+    AbsoluteTolerance: 0.05
+
+---
+
+Background: &solution Lorentzian
+
+InitialGuess: Zero
+
+RandomizeInitialGuess: None
+
+DomainCreator:
+  Sphere:
+    InnerRadius: 10
+    OuterRadius: &outer_radius 1e9
+    Interior:
+      FillWithSphericity: 0.
+    InitialRefinement: 0
+    InitialGridPoints: 6
+    UseEquiangularMap: True
+    EquatorialCompression: None
+    RadialPartitioning: [&outer_shell_inner_radius 20.]
+    RadialDistribution: [Linear, &outer_shell_distribution Inverse]
+    WhichWedges: All
+    TimeDependentMaps: None
+    OuterBoundaryCondition:
+      AnalyticSolution:
+        Solution: *solution
+        Field: Dirichlet
+
+Discretization:
+  DiscontinuousGalerkin:
+    PenaltyParameter: 1.
+    Massive: True
+    Quadrature: GaussLobatto
+    Formulation: WeakInertial
+
+Observers:
+  VolumeFileName: "LorentzianVolume"
+  ReductionFileName: "LorentzianReductions"
+
+LinearSolver:
+  ConvergenceCriteria:
+    MaxIterations: 10
+    RelativeResidual: 0.
+    AbsoluteResidual: 1.e-3
+  Verbosity: Quiet
+
+  Multigrid:
+    Iterations: 1
+    MaxLevels: Auto
+    PreSmoothing: True
+    PostSmoothingAtBottom: False
+    Verbosity: Silent
+    OutputVolumeData: True
+
+  SchwarzSmoother:
+    Iterations: 3
+    MaxOverlap: 2
+    Verbosity: Silent
+    SubdomainSolver:
+      ExplicitInverse:
+        WriteMatrixToFile: None
+    ObservePerCoreReductions: False
+
+RadiallyCompressedCoordinates:
+  InnerRadius: *outer_shell_inner_radius
+  OuterRadius: *outer_radius
+  Compression: *outer_shell_distribution
+
+EventsAndTriggers:
+  - Trigger: HasConverged
+    Events:
+      - ObserveNorms:
+          SubfileName: ErrorNorms
+          TensorsToObserve:
+            - Name: Error(Field)
+              NormType: L2Norm
+              Components: Sum
+      - ObserveFields:
+          SubfileName: VolumeData
+          VariablesToObserve:
+            - Field
+            - Error(Field)
+            - RadiallyCompressedCoordinates
+            - FixedSource(Field)
+          InterpolateToMesh: None
+          CoordinatesFloatingPointType: Double
+          FloatingPointTypes: [Double]
+
+PhaseChangeAndTriggers:
+
+BuildMatrix:
+  MatrixSubfileName: Matrix
+  Verbosity: Verbose
+
+Parallelization:
+  ElementDistribution: NumGridPoints
+
+ResourceInfo:
+  AvoidGlobalProc0: false
+  Singletons: Auto
diff --git a/tests/InputFiles/Poisson/ProductOfSinusoids3D.yaml b/tests/InputFiles/Poisson/ProductOfSinusoids3D.yaml
@@ -59,7 +59,7 @@ Discretization:
     PenaltyParameter: 1.
     Massive: False
     Quadrature: GaussLobatto
-    Formulation: StrongInertial
+    Formulation: WeakInertial
 
 Observers:
   VolumeFileName: "PoissonProductOfSinusoids3DVolume"

diff --git a/tests/InputFiles/Xcts/BinaryBlackHole.yaml b/tests/InputFiles/Xcts/BinaryBlackHole.yaml
@@ -117,7 +117,7 @@ Discretization:
     PenaltyParameter: 1.
     Massive: True
     Quadrature: GaussLobatto
-    Formulation: StrongInertial
+    Formulation: WeakInertial
 
 Observers:
   VolumeFileName: "BbhVolume"

diff --git a/tests/InputFiles/Xcts/HeadOnBns.yaml b/tests/InputFiles/Xcts/HeadOnBns.yaml
@@ -108,7 +108,7 @@ Discretization:
     PenaltyParameter: 1.
     Massive: True
     Quadrature: GaussLobatto
-    Formulation: StrongInertial
+    Formulation: WeakInertial
 
 Observers:
   VolumeFileName: "BnsVolume"

diff --git a/tests/InputFiles/Xcts/KerrSchild.yaml b/tests/InputFiles/Xcts/KerrSchild.yaml
@@ -72,7 +72,7 @@ Discretization:
     PenaltyParameter: 1.
     Massive: True
     Quadrature: GaussLobatto
-    Formulation: StrongInertial
+    Formulation: WeakInertial
 
 Observers:
   VolumeFileName: "KerrSchildVolume"

diff --git a/tests/InputFiles/Xcts/TovStar.yaml b/tests/InputFiles/Xcts/TovStar.yaml
@@ -62,7 +62,7 @@ Discretization:
     PenaltyParameter: 1.
     Massive: True
     Quadrature: GaussLobatto
-    Formulation: StrongInertial
+    Formulation: WeakInertial
 
 Observers:
   VolumeFileName: "TovStarVolume"