Add support for StrongLogical elliptic DG formulation

sxs-collaboration · Aug 24, 2024 · 4ccdca1 · 4ccdca1
1 parent 99e3b5e
commit 4ccdca1
Show file tree

Hide file tree

Showing 4 changed files with 325 additions and 61 deletions.
diff --git a/src/Elliptic/DiscontinuousGalerkin/DgOperator.hpp b/src/Elliptic/DiscontinuousGalerkin/DgOperator.hpp
@@ -626,6 +626,16 @@ 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 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);
+        *operator_applied_to_vars *= -1.;
+        if (not massive) {
+          *operator_applied_to_vars *= get(det_inv_jacobian);
+        }
       } else if (formulation == ::dg::Formulation::WeakInertial) {
         // Compute weak divergence:
         //   F^i \partial_i \phi = 1/w_p \sum_q
@@ -638,7 +648,8 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
       } else {
         ERROR("Unsupported DG formulation: "
               << formulation
-              << "\nSupported formulations are: StrongInertial, WeakInertial.");
+              << "\nSupported formulations are: StrongInertial, WeakInertial, "
+                 "StrongLogical.");
       }
       if constexpr (not std::is_same_v<SourcesComputer, void>) {
         Variables<tmpl::list<OperatorTags...>> sources{num_points, 0.};
@@ -807,10 +818,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/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}")