Elliptic DG: support modifying boundary data

This is used to solve for different variables in
some regions of the domain to absorb singularities.

src/Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp

@@ -28,6 +28,7 @@
 #include "Elliptic/DiscontinuousGalerkin/Initialization.hpp"
 #include "Elliptic/DiscontinuousGalerkin/Tags.hpp"
 #include "Elliptic/Systems/GetFluxesComputer.hpp"
+#include "Elliptic/Systems/GetModifyBoundaryData.hpp"
 #include "Elliptic/Systems/GetSourcesComputer.hpp"
 #include "Elliptic/Utilities/ApplyAt.hpp"
 #include "NumericalAlgorithms/DiscontinuousGalerkin/HasReceivedFromAllMortars.hpp"
@@ -498,6 +499,9 @@ template <typename System, bool Linearized, typename TemporalIdTag,
           typename PrimalMortarFieldsTag = PrimalFieldsTag,
           typename PrimalMortarFluxesTag = PrimalFluxesTag>
 struct DgOperator {
+  using system = System;
+  using temporal_id_tag = TemporalIdTag;
+
  private:
   static constexpr size_t Dim = System::volume_dim;
 
@@ -530,15 +534,17 @@ template <typename System, typename FixedSourcesTag,
           typename FluxesArgsTags =
               elliptic::get_fluxes_argument_tags<System, false>,
           typename SourcesArgsTags =
-              elliptic::get_sources_argument_tags<System, false>>
+              elliptic::get_sources_argument_tags<System, false>,
+          typename ModifyBoundaryDataArgsTags =
+              elliptic::get_modify_boundary_data_args_tags<System>>
 struct ImposeInhomogeneousBoundaryConditionsOnSource;
 
 /// \cond
 template <typename System, typename FixedSourcesTag, typename... FluxesArgsTags,
-          typename... SourcesArgsTags>
+          typename... SourcesArgsTags, typename... ModifyBoundaryDataArgsTags>
 struct ImposeInhomogeneousBoundaryConditionsOnSource<
     System, FixedSourcesTag, tmpl::list<FluxesArgsTags...>,
-    tmpl::list<SourcesArgsTags...>> {
+    tmpl::list<SourcesArgsTags...>, tmpl::list<ModifyBoundaryDataArgsTags...>> {
  private:
   static constexpr size_t Dim = System::volume_dim;
   using BoundaryConditionsBase = typename System::boundary_conditions_base;
@@ -636,12 +642,16 @@ struct ImposeInhomogeneousBoundaryConditionsOnSource<
                                                    Frame::Inertial>>>(box),
         db::get<::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>>(box),
         db::get<::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>>(box),
+        db::get<::Tags::Mortars<domain::Tags::DetSurfaceJacobian<
+                                    Frame::ElementLogical, Frame::Inertial>,
+                                Dim>>(box),
         db::get<::Tags::Mortars<elliptic::dg::Tags::PenaltyFactor, Dim>>(box),
         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);
+        fluxes_args_on_faces,
+        std::forward_as_tuple(db::get<ModifyBoundaryDataArgsTags>(box)...));
 
     // Move the mutated data back into the DataBox
     db::mutate<FixedSourcesTag>(

src/Elliptic/DiscontinuousGalerkin/DgOperator.hpp

@@ -517,11 +517,11 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
 
   // --- This is essentially a break to communicate the mortar data ---
 
-  template <bool AllDataIsZero, typename... OperatorTags,
-            typename... PrimalVars, typename... PrimalFluxesVars,
-            typename... PrimalMortarVars, typename... PrimalMortarFluxes,
-            typename TemporalId, typename... FluxesArgs,
-            typename... SourcesArgs, typename DataIsZero = NoDataIsZero,
+  template <typename... OperatorTags, typename... PrimalVars,
+            typename... PrimalFluxesVars, typename... PrimalMortarVars,
+            typename... PrimalMortarFluxes, typename TemporalId,
+            typename... FluxesArgs, typename... SourcesArgs,
+            typename DataIsZero = NoDataIsZero,
             typename DirectionsPredicate = AllDirections>
   static void apply_operator(
       const gsl::not_null<Variables<tmpl::list<OperatorTags...>>*>
@@ -686,11 +686,6 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
       const auto& neighbor_id = mortar_id.id();
       const bool is_internal =
           (neighbor_id != ElementId<Dim>::external_boundary_id());
-      if constexpr (AllDataIsZero) {
-        if (is_internal) {
-          continue;
-        }
-      }
       if (not directions_predicate(direction)) {
         continue;
       }
@@ -992,6 +987,7 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
 
   template <typename... FixedSourcesTags, typename ApplyBoundaryCondition,
             typename... FluxesArgs, typename... SourcesArgs,
+            typename... ModifyBoundaryDataArgs,
             bool LocalLinearized = Linearized,
             // This function adds nothing to the fixed sources if the operator
             // is linearized, so it shouldn't be used in that case
@@ -1016,23 +1012,23 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
           face_jacobian_times_inv_jacobians,
       const ::dg::MortarMap<Dim, Mesh<Dim - 1>>& all_mortar_meshes,
       const ::dg::MortarMap<Dim, ::dg::MortarSize<Dim - 1>>& all_mortar_sizes,
+      const ::dg::MortarMap<Dim, Scalar<DataVector>>& mortar_jacobians,
       const ::dg::MortarMap<Dim, Scalar<DataVector>>& penalty_factors,
       const bool massive, const ::dg::Formulation formulation,
       const ApplyBoundaryCondition& apply_boundary_condition,
       const std::tuple<FluxesArgs...>& fluxes_args,
       const std::tuple<SourcesArgs...>& sources_args,
-      const DirectionMap<Dim, std::tuple<FluxesArgs...>>&
-          fluxes_args_on_faces) {
+      const DirectionMap<Dim, std::tuple<FluxesArgs...>>& fluxes_args_on_faces,
+      const std::tuple<ModifyBoundaryDataArgs...>& modify_boundary_data_args) {
     // We just feed zero variables through the nonlinear operator to extract the
     // constant contribution at external boundaries. Since the variables are
     // zero the operator simplifies quite a lot. The simplification is probably
     // not very important for performance because this function will only be
     // called when solving a linear elliptic system and only once during
     // initialization, but we specialize the operator for zero data nonetheless
-    // just so we can ignore internal boundaries. For internal boundaries we
-    // would unnecessarily have to copy mortar data around to emulate the
-    // communication step, so by just skipping internal boundaries we avoid
-    // that.
+    // just so we can ignore internal boundaries. Only when the system modifies
+    // boundary data do we need to handle internal boundaries (see below),
+    // otherwise we can skip them.
     const size_t num_points = mesh.number_of_grid_points();
     const Variables<tmpl::list<PrimalFields...>> zero_primal_vars{num_points,
                                                                   0.};
@@ -1042,7 +1038,7 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
         unused_deriv_vars_buffer{};
     Variables<tmpl::list<FixedSourcesTags...>> operator_applied_to_zero_vars{
         num_points};
-    // Set up data on mortars. We only need them at external boundaries.
+    // Set up data on mortars
     ::dg::MortarMap<Dim, MortarData<size_t, tmpl::list<PrimalFields...>,
                                     tmpl::list<PrimalFluxes...>>>
         all_mortar_data{};
@@ -1054,15 +1050,53 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
                               element, mesh, inv_jacobian, face_normals,
                               all_mortar_meshes, all_mortar_sizes, temporal_id,
                               apply_boundary_condition, fluxes_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_normals, face_normal_vectors,
-        face_normal_magnitudes, face_jacobians,
-        face_jacobian_times_inv_jacobians, all_mortar_meshes, all_mortar_sizes,
-        {}, penalty_factors, massive, formulation, temporal_id,
-        fluxes_args_on_faces, sources_args);
+    // Modify internal boundary data if needed, e.g. to transform from one
+    // variable to another when crossing element boundaries. This is a nonlinear
+    // operation in the sense that feeding zero through the operator is nonzero,
+    // so we evaluate it here and add the contribution to the fixed sources,
+    // just like inhomogeneous external boundary conditions.
+    if constexpr (not std::is_same_v<typename System::modify_boundary_data,
+                                     void>) {
+      for (const auto& [direction, neighbors] : element.neighbors()) {
+        for (const auto& neighbor_id : neighbors) {
+          const ::dg::MortarId<Dim> mortar_id{direction, neighbor_id};
+          ASSERT(not all_mortar_data.contains(mortar_id),
+                 "Mortar data with ID " << mortar_id << " already exists.");
+          auto& mortar_data = all_mortar_data[mortar_id];
+          // Manufacture zero mortar data, store as local data, then treat as
+          // received from neighbor and apply modifications
+          BoundaryData<tmpl::list<PrimalFields...>, tmpl::list<PrimalFluxes...>>
+              remote_boundary_data_on_mortar{};
+          remote_boundary_data_on_mortar.field_data.initialize(
+              all_mortar_meshes.at(mortar_id).number_of_grid_points(), 0.);
+          mortar_data.local_insert(temporal_id, remote_boundary_data_on_mortar);
+          // Modify "received" mortar data
+          std::apply(
+              [&remote_boundary_data_on_mortar,
+               &mortar_id](const auto&... args) {
+                System::modify_boundary_data::apply(
+                    make_not_null(&get<PrimalFields>(
+                        remote_boundary_data_on_mortar.field_data))...,
+                    make_not_null(&get<::Tags::NormalDotFlux<PrimalFields>>(
+                        remote_boundary_data_on_mortar.field_data))...,
+                    mortar_id, args...);
+              },
+              modify_boundary_data_args);
+          // Insert as remote mortar data
+          mortar_data.remote_insert(temporal_id,
+                                    std::move(remote_boundary_data_on_mortar));
+        }
+      }
+    }
+    apply_operator(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_normals, face_normal_vectors, face_normal_magnitudes,
+                   face_jacobians, face_jacobian_times_inv_jacobians,
+                   all_mortar_meshes, all_mortar_sizes, mortar_jacobians,
+                   penalty_factors, massive, formulation, temporal_id,
+                   fluxes_args_on_faces, sources_args);
     // Impose the nonlinear (constant) boundary contribution as fixed sources on
     // the RHS of the equations
     *fixed_sources -= operator_applied_to_zero_vars;
@@ -1094,7 +1128,7 @@ void prepare_mortar_data(Args&&... args) {
  */
 template <typename System, bool Linearized, typename... Args>
 void apply_operator(Args&&... args) {
-  detail::DgOperatorImpl<System, Linearized>::template apply_operator<false>(
+  detail::DgOperatorImpl<System, Linearized>::apply_operator(
       std::forward<Args>(args)...);
 }
 

src/Elliptic/DiscontinuousGalerkin/Initialization.hpp

@@ -274,12 +274,14 @@ struct InitializeFacesAndMortars : tt::ConformsTo<::amr::protocols::Projector> {
                      // external boundaries.
                      ::Tags::deriv<domain::Tags::UnnormalizedFaceNormal<Dim>,
                                    tmpl::size_t<Dim>, Frame::Inertial>>>,
-      tmpl::list<::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>,
-                 ::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>,
-                 ::Tags::Mortars<domain::Tags::DetSurfaceJacobian<
-                                     Frame::ElementLogical, Frame::Inertial>,
-                                 Dim>,
-                 ::Tags::Mortars<elliptic::dg::Tags::PenaltyFactor, Dim>>>;
+      tmpl::list<
+          ::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>,
+          ::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>,
+          ::Tags::Mortars<domain::Tags::Coordinates<Dim, Frame::Inertial>, Dim>,
+          ::Tags::Mortars<domain::Tags::DetSurfaceJacobian<
+                              Frame::ElementLogical, Frame::Inertial>,
+                          Dim>,
+          ::Tags::Mortars<elliptic::dg::Tags::PenaltyFactor, Dim>>>;
   using argument_tags = tmpl::append<
       tmpl::list<domain::Tags::Mesh<Dim>, domain::Tags::Element<Dim>,
                  domain::Tags::NeighborMesh<Dim>, domain::Tags::ElementMap<Dim>,
@@ -318,6 +320,8 @@ struct InitializeFacesAndMortars : tt::ConformsTo<::amr::protocols::Projector> {
       const gsl::not_null<::dg::MortarMap<Dim, Mesh<Dim - 1>>*> mortar_meshes,
       const gsl::not_null<::dg::MortarMap<Dim, ::dg::MortarSize<Dim - 1>>*>
           mortar_sizes,
+      const gsl::not_null<::dg::MortarMap<Dim, tnsr::I<DataVector, Dim>>*>
+          all_mortar_inertial_coords,
       const gsl::not_null<::dg::MortarMap<Dim, Scalar<DataVector>>*>
           mortar_jacobians,
       const gsl::not_null<::dg::MortarMap<Dim, Scalar<DataVector>>*>
@@ -333,9 +337,10 @@ struct InitializeFacesAndMortars : tt::ConformsTo<::amr::protocols::Projector> {
     apply(face_directions, faces_inertial_coords, face_normals,
           face_normal_vectors, face_normal_magnitudes, face_jacobians,
           face_jacobian_times_inv_jacobian, deriv_unnormalized_face_normals,
-          mortar_meshes, mortar_sizes, mortar_jacobians, penalty_factors, mesh,
-          element, neighbor_meshes, element_map, inv_jacobian, domain,
-          functions_of_time, penalty_parameter, nullptr, nullptr, amr_data...);
+          mortar_meshes, mortar_sizes, all_mortar_inertial_coords,
+          mortar_jacobians, penalty_factors, mesh, element, neighbor_meshes,
+          element_map, inv_jacobian, domain, functions_of_time,
+          penalty_parameter, nullptr, nullptr, amr_data...);
   }
   template <typename Background, typename Metavariables, typename... AmrData>
   static void apply(
@@ -359,6 +364,8 @@ struct InitializeFacesAndMortars : tt::ConformsTo<::amr::protocols::Projector> {
       const gsl::not_null<::dg::MortarMap<Dim, Mesh<Dim - 1>>*> mortar_meshes,
       const gsl::not_null<::dg::MortarMap<Dim, ::dg::MortarSize<Dim - 1>>*>
           mortar_sizes,
+      const gsl::not_null<::dg::MortarMap<Dim, tnsr::I<DataVector, Dim>>*>
+          all_mortar_inertial_coords,
       const gsl::not_null<::dg::MortarMap<Dim, Scalar<DataVector>>*>
           mortar_jacobians,
       const gsl::not_null<::dg::MortarMap<Dim, Scalar<DataVector>>*>
@@ -450,6 +457,7 @@ struct InitializeFacesAndMortars : tt::ConformsTo<::amr::protocols::Projector> {
     // Mortars (internal directions)
     mortar_meshes->clear();
     mortar_sizes->clear();
+    all_mortar_inertial_coords->clear();
     mortar_jacobians->clear();
     penalty_factors->clear();
     const auto& element_id = element.id();
@@ -471,8 +479,11 @@ struct InitializeFacesAndMortars : tt::ConformsTo<::amr::protocols::Projector> {
         const auto& mortar_size = mortar_sizes->at(mortar_id);
         const auto mortar_logical_coords = detail::mortar_logical_coordinates(
             mortar_mesh, mortar_size, direction);
-        const auto mortar_inertial_coords =
-            element_map(mortar_logical_coords, 0., functions_of_time);
+        const auto& mortar_inertial_coords =
+            all_mortar_inertial_coords
+                ->emplace(mortar_id, element_map(mortar_logical_coords, 0.,
+                                                 functions_of_time))
+                .first->second;
         Scalar<DataVector> perpendicular_element_size{};
         if (Spectral::needs_projection(face_mesh, mortar_mesh, mortar_size)) {
           auto& mortar_jacobian = (*mortar_jacobians)[mortar_id];

src/Elliptic/Protocols/FirstOrderSystem.hpp

@@ -154,7 +154,7 @@ struct test_fields_and_fluxes<Dim, tmpl::list<PrimalFields...>,
  *   The function is expected to _add_ the sources \f$S_\alpha\f$ to the
  *   output buffers. It must also have the following alias:
  *
- * - `const_global_cache_tags`: the subset of `argument_tags` that can be
+ *   - `const_global_cache_tags`: the subset of `argument_tags` that can be
  *     retrieved from _any_ element's DataBox, because they are stored in the
  *     global cache.
  *
@@ -164,6 +164,35 @@ struct test_fields_and_fluxes<Dim, tmpl::list<PrimalFields...>,
  *   boundary conditions. Boundary conditions can be factory-created from this
  *   base class. Currently this should be a specialization of
  *   `elliptic::BoundaryConditions::BoundaryCondition`.
+ *
+ * - `modify_boundary_data` (optional, can be `void`): A class that can modify
+ *   boundary data received from a neighboring element, e.g. to transform from
+ *   one variable to another across element boundaries. This can be used to
+ *   solve for different variables in some regions of the domain to absorb
+ *   singularities. For example, when solving $-\Delta u = f$ we could define
+ *   $u=u_R + u_P$ in some region of the domain with a known (often singular)
+ *   field $u_P$ and solve only for the regular field $u_R$ in this region. So,
+ *   when receiving boundary data from outside this region, we subtract $u_P$,
+ *   and when receiving boundary data from inside this region, we add $u_P$. We
+ *   also add $\Delta u_P$ to the volume fixed sources $f$ inside the
+ *   regularized region.
+ *
+ *   The `modify_boundary_data` must have an `argument_tags` type alias and an
+ *   `apply` function that takes these arguments in this order:
+ *
+ *   1. The primal fields and the normal-dot-fluxes on the mortar as not-null
+ *      pointer. These hold the received data from the neighbor and can be
+ *      modified. Note: the normal-dot-fluxes have been computed with the
+ *      neighbor's normal, so from the perspective of the receiving element they
+ *      are $-n_i F^i$ where $n_i$ is the face normal of the receiving element.
+ *   2. The `const DirectionalId<Dim>& mortar_id` identifying the mortar.
+ *   3. The `argument_tags`.
+ *
+ *   Currently, modification made by this function must not depend on the
+ *   variables, meaning that the modification can only be adding or subtracting
+ *   a precomputed field. This is a simplification so the linearized operator is
+ *   not modified at all and can be relaxed if needed (then we also need
+ *   `modify_boundary_data_linearized`).
  */
 struct FirstOrderSystem {
   template <typename ConformingType>
@@ -191,6 +220,7 @@ struct FirstOrderSystem {
     static_assert(
         std::is_base_of_v<domain::BoundaryConditions::BoundaryCondition,
                           boundary_conditions_base>);
+    using modify_boundary_data = typename ConformingType::modify_boundary_data;
   };
 };
 

src/Elliptic/Systems/BnsInitialData/FirstOrderSystem.hpp

@@ -94,5 +94,6 @@ struct FirstOrderSystem
 
   using boundary_conditions_base =
       elliptic::BoundaryConditions::BoundaryCondition<3>;
+  using modify_boundary_data = void;
 };
 }  // namespace BnsInitialData

src/Elliptic/Systems/Elasticity/FirstOrderSystem.hpp

@@ -54,5 +54,6 @@ struct FirstOrderSystem
 
   using boundary_conditions_base =
       elliptic::BoundaryConditions::BoundaryCondition<Dim>;
+  using modify_boundary_data = void;
 };
 }  // namespace Elasticity

src/Elliptic/Systems/GetModifyBoundaryData.hpp

@@ -0,0 +1,24 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <type_traits>
+
+#include "Utilities/TMPL.hpp"
+
+namespace elliptic {
+namespace detail {
+struct NoModifyBoundaryData {
+  using argument_tags = tmpl::list<>;
+};
+}  // namespace detail
+
+/// The `argument_tags` of the `System::modify_boundary_data`, or an empty list
+/// if `System::modify_boundary_data` is `void`.
+template <typename System>
+using get_modify_boundary_data_args_tags = typename tmpl::conditional_t<
+    std::is_same_v<typename System::modify_boundary_data, void>,
+    detail::NoModifyBoundaryData,
+    typename System::modify_boundary_data>::argument_tags;
+}  // namespace elliptic