Add vtu slicing to self-force executables

src/Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp

@@ -298,18 +298,21 @@ template <typename System, bool Linearized, typename TemporalIdTag,
           typename FluxesArgsTags =
               elliptic::get_fluxes_argument_tags<System, Linearized>,
           typename SourcesArgsTags =
-              elliptic::get_sources_argument_tags<System, Linearized>>
+              elliptic::get_sources_argument_tags<System, Linearized>,
+          typename ModifyBoundaryDataArgsTags =
+              elliptic::get_modify_boundary_data_args_tags<System, Linearized>>
 struct ReceiveMortarDataAndApplyOperator;
 
 template <typename System, bool Linearized, typename TemporalIdTag,
           typename PrimalFieldsTag, typename PrimalFluxesTag,
           typename OperatorAppliedToFieldsTag, typename PrimalMortarFieldsTag,
           typename PrimalMortarFluxesTag, typename... FluxesArgsTags,
-          typename... SourcesArgsTags>
+          typename... SourcesArgsTags, typename... ModifyBoundaryDataArgsTags>
 struct ReceiveMortarDataAndApplyOperator<
     System, Linearized, TemporalIdTag, PrimalFieldsTag, PrimalFluxesTag,
     OperatorAppliedToFieldsTag, PrimalMortarFieldsTag, PrimalMortarFluxesTag,
-    tmpl::list<FluxesArgsTags...>, tmpl::list<SourcesArgsTags...>> {
+    tmpl::list<FluxesArgsTags...>, tmpl::list<SourcesArgsTags...>,
+    tmpl::list<ModifyBoundaryDataArgsTags...>> {
  private:
   static constexpr size_t Dim = System::volume_dim;
   using all_mortar_data_tag = ::Tags::Mortars<
@@ -416,7 +419,8 @@ struct ReceiveMortarDataAndApplyOperator<
         db::get<elliptic::dg::Tags::Massive>(box),
         db::get<elliptic::dg::Tags::Formulation>(box), temporal_id,
         fluxes_args_on_faces,
-        std::forward_as_tuple(db::get<SourcesArgsTags>(box)...));
+        std::forward_as_tuple(db::get<SourcesArgsTags>(box)...),
+        std::forward_as_tuple(db::get<ModifyBoundaryDataArgsTags>(box)...));
 
     // Increment temporal ID
     db::mutate<TemporalIdTag>(
@@ -552,7 +556,7 @@ template <typename System, typename FixedSourcesTag,
           typename SourcesArgsTags =
               elliptic::get_sources_argument_tags<System, false>,
           typename ModifyBoundaryDataArgsTags =
-              elliptic::get_modify_boundary_data_args_tags<System>>
+              elliptic::get_modify_boundary_data_args_tags<System, false>>
 struct ImposeInhomogeneousBoundaryConditionsOnSource;
 
 /// \cond

src/Elliptic/DiscontinuousGalerkin/DgOperator.hpp

@@ -522,7 +522,8 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
             typename... DerivVars, typename... PrimalFluxesVars,
             typename... PrimalMortarVars, typename... PrimalMortarFluxes,
             typename TemporalId, typename... FluxesArgs,
-            typename... SourcesArgs, typename DataIsZero = NoDataIsZero,
+            typename... SourcesArgs, typename... ModifyBoundaryDataArgs,
+            typename DataIsZero = NoDataIsZero,
             typename DirectionsPredicate = AllDirections>
   static void apply_operator(
       const gsl::not_null<Variables<tmpl::list<OperatorTags...>>*>
@@ -561,6 +562,7 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
       const TemporalId& /*temporal_id*/,
       const DirectionMap<Dim, std::tuple<FluxesArgs...>>& fluxes_args_on_faces,
       const std::tuple<SourcesArgs...>& sources_args,
+      const std::tuple<ModifyBoundaryDataArgs...>& modify_boundary_data_args,
       const DataIsZero& data_is_zero = NoDataIsZero{},
       const DirectionsPredicate& directions_predicate = AllDirections{}) {
     static_assert(
@@ -702,6 +704,30 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
 
       const auto face_mesh = mesh.slice_away(direction.dimension());
       auto [local_data, remote_data] = mortar_data.extract();
+
+      if (is_internal) {
+        if constexpr (Linearized and
+                      not std::is_same_v<typename System::modify_boundary_data,
+                                         void>) {
+          // Apply a linearized modification to received boundary data.
+          // This allows modifications to depend linearly on the variables.
+          std::apply(
+              [&rem_data = remote_data, &loc_data = local_data,
+               m_id = mortar_id](const auto&... args) {
+                System::modify_boundary_data::apply_linearized(
+                    make_not_null(
+                        &get<PrimalMortarVars>(rem_data.field_data))...,
+                    make_not_null(&get<::Tags::NormalDotFlux<PrimalMortarVars>>(
+                        rem_data.field_data))...,
+                    get<PrimalMortarVars>(loc_data.field_data)...,
+                    get<::Tags::NormalDotFlux<PrimalMortarVars>>(
+                        loc_data.field_data)...,
+                    m_id, args...);
+              },
+              modify_boundary_data_args);
+        }
+      }
+
       const size_t face_num_points = face_mesh.number_of_grid_points();
       const auto& face_normal = face_normals.at(direction);
       const auto& face_normal_vector = face_normal_vectors.at(direction);
@@ -1098,7 +1124,7 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
         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);
+        fluxes_args_on_faces, sources_args, modify_boundary_data_args);
     // Impose the nonlinear (constant) boundary contribution as fixed sources on
     // the RHS of the equations
     *fixed_sources -= operator_applied_to_zero_vars;

src/Elliptic/DiscontinuousGalerkin/SubdomainOperator/SubdomainOperator.hpp

@@ -30,6 +30,7 @@
 #include "Elliptic/DiscontinuousGalerkin/SubdomainOperator/Tags.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/MortarHelpers.hpp"
@@ -182,18 +183,22 @@ struct SubdomainOperator
       elliptic::get_fluxes_argument_tags<System, linearized>;
   using sources_args_tags =
       elliptic::get_sources_argument_tags<System, linearized>;
-
+  using modify_boundary_data_args_tags =
+      elliptic::get_modify_boundary_data_args_tags<System, linearized>;
   // We need the fluxes args also on interfaces (internal and external). The
   // volume tags are the subset that don't have to be taken from interfaces.
   using fluxes_args_volume_tags =
       elliptic::get_fluxes_volume_tags<System, linearized>;
 
   // 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<
+  using args_tags_from_center = tmpl::remove_duplicates<tmpl::append<
       elliptic::get_fluxes_const_global_cache_tags<System, linearized>,
       elliptic::get_sources_const_global_cache_tags<System, linearized>,
-      elliptic::dg::Tags::Massive, elliptic::dg::Tags::Formulation>>;
+      elliptic::get_modify_boundary_data_const_global_cache_tags<System,
+                                                                 linearized>,
+      tmpl::list<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
@@ -210,6 +215,8 @@ struct SubdomainOperator
       tmpl::transform<fluxes_args_tags, make_overlap_tag>;
   using sources_args_tags_overlap =
       tmpl::transform<sources_args_tags, make_overlap_tag>;
+  using modify_boundary_data_args_tags_overlap =
+      tmpl::transform<modify_boundary_data_args_tags, make_overlap_tag>;
   using fluxes_args_tags_overlap_faces = tmpl::transform<
       domain::make_faces_tags<Dim, fluxes_args_tags, fluxes_args_volume_tags>,
       make_overlap_tag>;
@@ -279,6 +286,8 @@ struct SubdomainOperator
         db::apply<tags_to_retrieve>(get_items, box);
     const auto fluxes_args = db::apply<fluxes_args_tags>(get_items, box);
     const auto sources_args = db::apply<sources_args_tags>(get_items, box);
+    const auto modify_boundary_data_args =
+        db::apply<modify_boundary_data_args_tags>(get_items, box);
     using FluxesArgs = std::decay_t<decltype(fluxes_args)>;
     DirectionMap<Dim, FluxesArgs> fluxes_args_on_faces{};
     for (const auto& direction : Direction<Dim>::all_directions()) {
@@ -623,7 +632,8 @@ struct SubdomainOperator
               make_not_null(&central_mortar_data_), operand.element_data,
               central_deriv_vars_, central_primal_fluxes_, args...);
         },
-        box, temporal_id, fluxes_args_on_faces, sources_args, data_is_zero);
+        box, temporal_id, fluxes_args_on_faces, sources_args,
+        modify_boundary_data_args, data_is_zero);
     // Apply on neighbors
     for (const auto& [direction, neighbors] : central_element.neighbors()) {
       for (const auto& neighbor_id : neighbors) {
@@ -663,6 +673,9 @@ struct SubdomainOperator
             elliptic::util::apply_at<sources_args_tags_overlap,
                                      args_tags_from_center>(get_items, box,
                                                             overlap_id),
+            elliptic::util::apply_at<modify_boundary_data_args_tags_overlap,
+                                     args_tags_from_center>(get_items, box,
+                                                            overlap_id),
             data_is_zero);
 
         // Restrict the extended operator data back to the subdomain, assuming

src/Elliptic/Protocols/FirstOrderSystem.hpp

@@ -189,11 +189,19 @@ struct test_fields_and_fluxes<Dim, tmpl::list<PrimalFields...>,
  *   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`).
+ *   The `modify_boundary_data` may also be used for modifications that depend
+ *   linearly on the variables. In this case, the class must also provide an
+ *   `argument_tags_linearized` type alias and an `apply_linearized` function
+ *   that takes these arguments in this order:
+ *
+ *   1. The remote primal fields and the remote normal-dot-fluxes on the mortar
+ *      as not-null pointers. These hold the received data from the neighbor and
+ *      can be modified in-place.
+ *   2. The local primal fields and the local normal-dot-fluxes on the mortar
+ *      as const references. These hold the data from the local side of the
+ *      boundary.
+ *   3. The `const DirectionalId<Dim>& mortar_id` identifying the mortar.
+ *   4. The `argument_tags_linearized`.
  */
 struct FirstOrderSystem {
   template <typename ConformingType>

src/Elliptic/Systems/GetModifyBoundaryData.hpp

@@ -11,14 +11,26 @@ namespace elliptic {
 namespace detail {
 struct NoModifyBoundaryData {
   using argument_tags = tmpl::list<>;
+  using argument_tags_linearized = tmpl::list<>;
+  using const_global_cache_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<
+using get_modify_boundary_data = tmpl::conditional_t<
     std::is_same_v<typename System::modify_boundary_data, void>,
-    detail::NoModifyBoundaryData,
-    typename System::modify_boundary_data>::argument_tags;
+    detail::NoModifyBoundaryData, typename System::modify_boundary_data>;
+
+template <typename System, bool Linearized>
+using get_modify_boundary_data_args_tags = tmpl::conditional_t<
+    Linearized,
+    typename get_modify_boundary_data<System>::argument_tags_linearized,
+    typename get_modify_boundary_data<System>::argument_tags>;
+
+template <typename System, bool Linearized>
+using get_modify_boundary_data_const_global_cache_tags =
+    typename get_modify_boundary_data<System>::const_global_cache_tags;
+
 }  // namespace elliptic

src/Elliptic/Systems/SelfForce/GeneralRelativity/Actions/InitializeEffectiveSource.hpp

@@ -73,7 +73,8 @@ struct InitializeEffectiveSource : tt::ConformsTo<::amr::protocols::Projector> {
 
  public:  // Iterable action
   using const_global_cache_tags =
-      tmpl::list<elliptic::dg::Tags::Massive, BackgroundTag>;
+      tmpl::list<elliptic::dg::Tags::Massive, BackgroundTag,
+                 Tags::NullSlicingBlocks<Dim>>;
   using simple_tags =
       tmpl::list<fixed_sources_tag, singular_vars_tag,
                  ::Tags::Mortars<singular_vars_on_mortars_tag, Dim>,

src/Elliptic/Systems/SelfForce/GeneralRelativity/AnalyticData/CircularOrbit.cpp

@@ -43,6 +43,13 @@ tnsr::I<double, 2> CircularOrbit::puncture_position() const {
   return tnsr::I<double, 2>{{{r_star, M_PI_2}}};
 }
 
+double CircularOrbit::omega() const {
+  const double M = black_hole_mass_;
+  const double a = black_hole_spin_ * M;
+  const double r_0 = orbital_radius_;
+  return 1. / (a + sqrt(cube(r_0) / M));
+}
+
 // Background
 tuples::TaggedTuple<Tags::Alpha, Tags::Beta, Tags::GammaRstar, Tags::GammaTheta>
 CircularOrbit::variables(const tnsr::I<DataVector, 2>& x,
@@ -52,8 +59,7 @@ CircularOrbit::variables(const tnsr::I<DataVector, 2>& x,
   const double M = black_hole_mass_;
   const double r_plus = M * (1. + sqrt(1. - square(black_hole_spin_)));
   const double r_minus = M * (1. - sqrt(1. - square(black_hole_spin_)));
-  const double r_0 = orbital_radius_;
-  const double omega = 1. / (a + sqrt(cube(r_0) / M));
+  const double omega = this->omega();
   const auto& r_star = get<0>(x);
   const auto& theta = get<1>(x);
   const DataVector cos_theta = cos(theta);

src/Elliptic/Systems/SelfForce/GeneralRelativity/AnalyticData/CircularOrbit.hpp

@@ -89,6 +89,7 @@ class CircularOrbit : public elliptic::analytic_data::Background,
   double black_hole_spin() const { return black_hole_spin_; }
   double orbital_radius() const { return orbital_radius_; }
   int m_mode_number() const { return m_mode_number_; }
+  double omega() const;
 
   using background_tags =
       tmpl::list<Tags::Alpha, Tags::Beta, Tags::GammaRstar, Tags::GammaTheta>;

src/Elliptic/Systems/SelfForce/GeneralRelativity/CMakeLists.txt

@@ -26,6 +26,7 @@ target_link_libraries(
   DataStructures
   DomainStructure
   ErrorHandling
+  GrSelfForceAnalyticData
   Utilities
   )
 

src/Elliptic/Systems/SelfForce/GeneralRelativity/Equations.cpp

@@ -3,11 +3,14 @@
 
 #include "Elliptic/Systems/SelfForce/GeneralRelativity/Equations.hpp"
 
+#include <algorithm>
 #include <cstddef>
 
 #include "DataStructures/ComplexDataVector.hpp"
 #include "DataStructures/DataVector.hpp"
 #include "DataStructures/Tensor/Tensor.hpp"
+#include "Elliptic/Systems/SelfForce/GeneralRelativity/AnalyticData/CircularOrbit.hpp"
+#include "Utilities/Algorithm.hpp"
 
 namespace GrSelfForce {
 
@@ -104,4 +107,32 @@ void ModifyBoundaryData::apply(
   }
 }
 
+void ModifyBoundaryData::apply_linearized(
+    const gsl::not_null<tnsr::aa<ComplexDataVector, 3>*> field_remote,
+    const gsl::not_null<tnsr::aa<ComplexDataVector, 3>*>
+        n_dot_field_gradient_remote,
+    const tnsr::aa<ComplexDataVector, 3>& field_local,
+    const tnsr::aa<ComplexDataVector, 3>& /*n_dot_field_gradient_local*/,
+    const DirectionalId<Dim>& mortar_id, const Element<Dim>& element,
+    const std::vector<size_t>& null_slicing_blocks,
+    const elliptic::analytic_data::Background& background) {
+  if (alg::found(null_slicing_blocks, element.id().block_id()) ==
+      alg::found(null_slicing_blocks, mortar_id.id().block_id())) {
+    // Both elements use the same slicing. Nothing to do.
+    return;
+  }
+  // Apply the jump in the field gradient across the boundary to handle
+  // vtu-slicing. The signs are all the same (on both sides of the boundary and
+  // at both transition points).
+  const auto& circular_orbit =
+      dynamic_cast<const GrSelfForce::AnalyticData::CircularOrbit&>(background);
+  const double omega = circular_orbit.omega();
+  const double m_mode_number = circular_orbit.m_mode_number();
+  for (size_t j = 0; j < n_dot_field_gradient_remote->size(); ++j) {
+    (*n_dot_field_gradient_remote)[j] -=
+        std::complex<double>(0.0, m_mode_number * omega) *
+        ((field_local)[j] + (*field_remote)[j]) * 0.5;
+  }
+}
+
 }  // namespace GrSelfForce

src/Elliptic/Systems/SelfForce/GeneralRelativity/Equations.hpp

@@ -11,9 +11,12 @@
 #include "DataStructures/Tensor/Tensor.hpp"
 #include "DataStructures/Variables.hpp"
 #include "DataStructures/VariablesTag.hpp"
+#include "Domain/Tags.hpp"
 #include "Elliptic/Systems/SelfForce/GeneralRelativity/Tags.hpp"
+#include "Elliptic/Tags.hpp"
 #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
 #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
+#include "PointwiseFunctions/InitialDataUtilities/Background.hpp"
 #include "Utilities/Gsl.hpp"
 #include "Utilities/MakeWithValue.hpp"
 #include "Utilities/TMPL.hpp"
@@ -122,13 +125,29 @@ struct ModifyBoundaryData {
       tmpl::list<Tags::FieldIsRegularized,
                  ::Tags::Mortars<Tags::FieldIsRegularized, Dim>,
                  ::Tags::Mortars<singular_vars_on_mortars_tag, Dim>>;
+ public:
+  using argument_tags_linearized = tmpl::list<
+      domain::Tags::Element<Dim>, Tags::NullSlicingBlocks<Dim>,
+      elliptic::Tags::Background<elliptic::analytic_data::Background>>;
+  using const_global_cache_tags = tmpl::list<
+      Tags::NullSlicingBlocks<Dim>,
+      elliptic::Tags::Background<elliptic::analytic_data::Background>>;
   static void apply(
       gsl::not_null<tnsr::aa<ComplexDataVector, 3>*> field,
       gsl::not_null<tnsr::aa<ComplexDataVector, 3>*> n_dot_flux,
       const DirectionalId<Dim>& mortar_id, bool field_is_regularized,
       const DirectionalIdMap<Dim, bool>& neighbors_field_is_regularized,
       const DirectionalIdMap<Dim, typename singular_vars_on_mortars_tag::type>&
           singular_vars_on_mortars);
+  static void apply_linearized(
+      gsl::not_null<tnsr::aa<ComplexDataVector, 3>*> field_remote,
+      gsl::not_null<tnsr::aa<ComplexDataVector, 3>*>
+          n_dot_field_gradient_remote,
+      const tnsr::aa<ComplexDataVector, 3>& field_local,
+      const tnsr::aa<ComplexDataVector, 3>& n_dot_field_gradient_local,
+      const DirectionalId<Dim>& mortar_id, const Element<Dim>& element,
+      const std::vector<size_t>& null_slicing_blocks,
+      const elliptic::analytic_data::Background& background);
 };
 
 }  // namespace GrSelfForce