Merge pull request #5915 from nilsvu/xcts_robin

XCTS: add Robin outer boundary conditions

src/Elliptic/BoundaryConditions/AnalyticSolution.hpp

@@ -134,6 +134,9 @@ class AnalyticSolution<System, Dim, tmpl::list<FieldTags...>,
   template <typename Metavariables>
   void apply(const gsl::not_null<typename FieldTags::type*>... fields,
              const gsl::not_null<typename FieldTags::type*>... n_dot_fluxes,
+             const TensorMetafunctions::prepend_spatial_index<
+                 typename FieldTags::type, Dim, UpLo::Lo,
+                 Frame::Inertial>&... /*deriv_fields*/,
              const Metavariables& /*meta*/,
              const tnsr::I<DataVector, Dim>& face_inertial_coords,
              const tnsr::i<DataVector, Dim>& face_normal) const {
@@ -184,7 +187,10 @@ class AnalyticSolution<System, Dim, tmpl::list<FieldTags...>,
 
   void apply_linearized(
       const gsl::not_null<typename FieldTags::type*>... fields,
-      const gsl::not_null<typename FieldTags::type*>... n_dot_fluxes) const {
+      const gsl::not_null<typename FieldTags::type*>... n_dot_fluxes,
+      const TensorMetafunctions::prepend_spatial_index<
+          typename FieldTags::type, Dim, UpLo::Lo,
+          Frame::Inertial>&... /*deriv_fields*/) const {
     const auto impose_boundary_condition = [this](auto field_tag_v,
                                                   const auto field,
                                                   const auto n_dot_flux) {

src/Elliptic/BoundaryConditions/ApplyBoundaryCondition.hpp

@@ -68,7 +68,7 @@ void apply_boundary_condition(
     const elliptic::BoundaryConditions::BoundaryCondition<Dim>&
         boundary_condition,
     const db::DataBox<DbTagsList>& box, const MapKeys& map_keys_to_direction,
-    const gsl::not_null<FieldsAndFluxes*>... fields_and_fluxes) {
+    FieldsAndFluxes&&... fields_and_fluxes) {
   using boundary_condition_classes =
       typename detail::GetBoundaryConditionClasses<
           elliptic::BoundaryConditions::BoundaryCondition<Dim>,
@@ -99,9 +99,13 @@ void apply_boundary_condition(
                                  volume_tags_transformed>(
             [&derived, &fields_and_fluxes...](const auto&... args) {
               if constexpr (Linearized) {
-                derived->apply_linearized(fields_and_fluxes..., args...);
+                derived->apply_linearized(
+                    std::forward<FieldsAndFluxes>(fields_and_fluxes)...,
+                    args...);
               } else {
-                derived->apply(fields_and_fluxes..., args...);
+                derived->apply(
+                    std::forward<FieldsAndFluxes>(fields_and_fluxes)...,
+                    args...);
               }
             },
             box, map_keys_to_direction);

src/Elliptic/BoundaryConditions/BoundaryCondition.hpp

@@ -68,12 +68,11 @@ namespace BoundaryConditions {
  *   1. The dynamic fields as not-null pointers.
  *   2. The normal-dot-fluxes corresponding to the dynamic fields as not-null
  *      pointers. These have the same types as the dynamic fields.
- *   3. The types held by the argument tags.
+ *   3. The field derivatives as const-refs.
+ *   4. The types held by the argument tags.
  *
- *   For first-order systems that involve auxiliary variables, only the
- *   non-auxiliary ("primal") variables are included in the lists above. For
- *   example, boundary conditions for a first-order Poisson system might have an
- *   `apply` function signature that looks like this:
+ *   For example, boundary conditions for a first-order Poisson system might
+ *   have an `apply` function signature that looks like this:
  *
  *   \snippet Elliptic/BoundaryConditions/Test_BoundaryCondition.cpp example_poisson_fields
  *

src/Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp

@@ -181,7 +181,7 @@ struct PrepareAndSendMortarData<
             element_id.block_id());
     const auto apply_boundary_condition =
         [&box, &boundary_conditions, &element_id](
-            const Direction<Dim>& direction, const auto... fields_and_fluxes) {
+            const Direction<Dim>& direction, auto&&... fields_and_fluxes) {
           ASSERT(
               boundary_conditions.contains(direction),
               "No boundary condition is available in block "
@@ -198,7 +198,8 @@ struct PrepareAndSendMortarData<
               dynamic_cast<const BoundaryConditionsBase&>(
                   *boundary_conditions.at(direction));
           elliptic::apply_boundary_condition<Linearized>(
-              boundary_condition, box, direction, fields_and_fluxes...);
+              boundary_condition, box, direction,
+              std::forward<decltype(fields_and_fluxes)>(fields_and_fluxes)...);
         };
 
     // Can't `db::get` the arguments for the boundary conditions within
@@ -564,7 +565,7 @@ struct ImposeInhomogeneousBoundaryConditionsOnSource<
             element_id.block_id());
     const auto apply_boundary_condition =
         [&box, &boundary_conditions, &element_id](
-            const Direction<Dim>& direction, const auto... fields_and_fluxes) {
+            const Direction<Dim>& direction, auto&&... fields_and_fluxes) {
           ASSERT(
               boundary_conditions.contains(direction),
               "No boundary condition is available in block "
@@ -581,7 +582,8 @@ struct ImposeInhomogeneousBoundaryConditionsOnSource<
               dynamic_cast<const BoundaryConditionsBase&>(
                   *boundary_conditions.at(direction));
           elliptic::apply_boundary_condition<false>(
-              boundary_condition, box, direction, fields_and_fluxes...);
+              boundary_condition, box, direction,
+              std::forward<decltype(fields_and_fluxes)>(fields_and_fluxes)...);
         };
 
     // Can't `db::get` the arguments for the boundary conditions within

src/Elliptic/DiscontinuousGalerkin/DgOperator.hpp

@@ -455,11 +455,21 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
         // the boundary conditions is taken from the "interior" side of the
         // boundary, i.e. with a normal vector that points _out_ of the
         // computational domain.
+        Variables<tmpl::list<DerivTags...>> deriv_vars_on_boundary{};
+        if (AllDataIsZero or local_data_is_zero) {
+          deriv_vars_on_boundary.initialize(face_num_points, 0.);
+        } else {
+          deriv_vars_on_boundary.initialize(face_num_points);
+          ::dg::project_contiguous_data_to_boundary(
+              make_not_null(&deriv_vars_on_boundary), *deriv_vars, mesh,
+              direction);
+        }
         apply_boundary_condition(
             direction,
             make_not_null(&get<PrimalMortarVars>(boundary_data.field_data))...,
             make_not_null(&get<::Tags::NormalDotFlux<PrimalMortarVars>>(
-                boundary_data.field_data))...);
+                boundary_data.field_data))...,
+            get<DerivTags>(deriv_vars_on_boundary)...);
 
         // Invert the sign of the fluxes to account for the inverted normal on
         // exterior faces. Also multiply by 2 and add the interior fluxes to

src/Elliptic/DiscontinuousGalerkin/SubdomainOperator/SubdomainOperator.hpp

@@ -289,7 +289,7 @@ struct SubdomainOperator
          &override_boundary_conditions](const ElementId<Dim>& local_element_id,
                                         const Direction<Dim>& local_direction,
                                         auto is_overlap, const auto& map_keys,
-                                        const auto... fields_and_fluxes) {
+                                        auto&&... fields_and_fluxes) {
           constexpr bool is_overlap_v =
               std::decay_t<decltype(is_overlap)>::value;
           // Get boundary conditions from domain, or use overridden boundary
@@ -336,8 +336,9 @@ struct SubdomainOperator
           elliptic::apply_boundary_condition<
               linearized,
               tmpl::conditional_t<is_overlap_v, make_overlap_tag, void>,
-              BoundaryConditionClasses>(boundary_condition, box, map_keys,
-                                        fields_and_fluxes...);
+              BoundaryConditionClasses>(
+              boundary_condition, box, map_keys,
+              std::forward<decltype(fields_and_fluxes)>(fields_and_fluxes)...);
         };
 
     // Check if the subdomain data is sparse, i.e. if some elements have zero
@@ -393,10 +394,11 @@ struct SubdomainOperator
     const auto apply_boundary_condition_center =
         [&apply_boundary_condition, &local_central_element = central_element](
             const Direction<Dim>& local_direction,
-            const auto... fields_and_fluxes) {
-          apply_boundary_condition(local_central_element.id(), local_direction,
-                                   std::false_type{}, local_direction,
-                                   fields_and_fluxes...);
+            auto&&... fields_and_fluxes) {
+          apply_boundary_condition(
+              local_central_element.id(), local_direction, std::false_type{},
+              local_direction,
+              std::forward<decltype(fields_and_fluxes)>(fields_and_fluxes)...);
         };
     db::apply<prepare_args_tags>(
         [this, &operand](const auto&... args) {
@@ -467,11 +469,12 @@ struct SubdomainOperator
         const auto apply_boundary_condition_neighbor =
             [&apply_boundary_condition, &local_neighbor_id = neighbor_id,
              &overlap_id](const Direction<Dim>& local_direction,
-                          const auto... fields_and_fluxes) {
+                          auto&&... fields_and_fluxes) {
               apply_boundary_condition(
                   local_neighbor_id, local_direction, std::true_type{},
                   std::forward_as_tuple(overlap_id, local_direction),
-                  fields_and_fluxes...);
+                  std::forward<decltype(fields_and_fluxes)>(
+                      fields_and_fluxes)...);
             };
 
         const auto fluxes_args_on_overlap =

src/Elliptic/Systems/Elasticity/BoundaryConditions/LaserBeam.cpp

@@ -15,6 +15,7 @@ namespace Elasticity::BoundaryConditions {
 void LaserBeam::apply(
     const gsl::not_null<tnsr::I<DataVector, 3>*> /*displacement*/,
     const gsl::not_null<tnsr::I<DataVector, 3>*> n_dot_minus_stress,
+    const tnsr::iJ<DataVector, 3>& /*deriv_displacement*/,
     const tnsr::I<DataVector, 3>& x,
     const tnsr::i<DataVector, 3>& face_normal) const {
   const auto n_dot_x = get<0>(face_normal) * get<0>(x) +
@@ -31,7 +32,8 @@ void LaserBeam::apply(
 
 void LaserBeam::apply_linearized(
     const gsl::not_null<tnsr::I<DataVector, 3>*> /*displacement*/,
-    const gsl::not_null<tnsr::I<DataVector, 3>*> n_dot_minus_stress) {
+    const gsl::not_null<tnsr::I<DataVector, 3>*> n_dot_minus_stress,
+    const tnsr::iJ<DataVector, 3>& /*deriv_displacement*/) {
   get<0>(*n_dot_minus_stress) = 0.;
   get<1>(*n_dot_minus_stress) = 0.;
   get<2>(*n_dot_minus_stress) = 0.;

src/Elliptic/Systems/Elasticity/BoundaryConditions/LaserBeam.hpp

@@ -102,6 +102,7 @@ class LaserBeam : public elliptic::BoundaryConditions::BoundaryCondition<3> {
 
   void apply(gsl::not_null<tnsr::I<DataVector, 3>*> displacement,
              gsl::not_null<tnsr::I<DataVector, 3>*> n_dot_minus_stress,
+             const tnsr::iJ<DataVector, 3>& deriv_displacement,
              const tnsr::I<DataVector, 3>& x,
              const tnsr::i<DataVector, 3>& face_normal) const;
 
@@ -110,7 +111,8 @@ class LaserBeam : public elliptic::BoundaryConditions::BoundaryCondition<3> {
 
   static void apply_linearized(
       gsl::not_null<tnsr::I<DataVector, 3>*> displacement,
-      gsl::not_null<tnsr::I<DataVector, 3>*> n_dot_minus_stress);
+      gsl::not_null<tnsr::I<DataVector, 3>*> n_dot_minus_stress,
+      const tnsr::iJ<DataVector, 3>& deriv_displacement);
 
   // NOLINTNEXTLINE(google-runtime-references)
   void pup(PUP::er& p) override { p | beam_width_; }

src/Elliptic/Systems/Elasticity/BoundaryConditions/Zero.cpp

@@ -28,7 +28,8 @@ std::string Zero<Dim, BoundaryConditionType>::name() {
 template <size_t Dim, elliptic::BoundaryConditionType BoundaryConditionType>
 void Zero<Dim, BoundaryConditionType>::apply(
     const gsl::not_null<tnsr::I<DataVector, Dim>*> displacement,
-    const gsl::not_null<tnsr::I<DataVector, Dim>*> n_dot_minus_stress) {
+    const gsl::not_null<tnsr::I<DataVector, Dim>*> n_dot_minus_stress,
+    const tnsr::iJ<DataVector, Dim>& /*deriv_displacement*/) {
   if constexpr (BoundaryConditionType ==
                 elliptic::BoundaryConditionType::Dirichlet) {
     (void)n_dot_minus_stress;
@@ -46,8 +47,9 @@ void Zero<Dim, BoundaryConditionType>::apply(
 template <size_t Dim, elliptic::BoundaryConditionType BoundaryConditionType>
 void Zero<Dim, BoundaryConditionType>::apply_linearized(
     const gsl::not_null<tnsr::I<DataVector, Dim>*> displacement,
-    const gsl::not_null<tnsr::I<DataVector, Dim>*> n_dot_minus_stress) {
-  apply(displacement, n_dot_minus_stress);
+    const gsl::not_null<tnsr::I<DataVector, Dim>*> n_dot_minus_stress,
+    const tnsr::iJ<DataVector, Dim>& deriv_displacement) {
+  apply(displacement, n_dot_minus_stress, deriv_displacement);
 }
 
 #define DIM(data) BOOST_PP_TUPLE_ELEM(0, data)

src/Elliptic/Systems/Elasticity/BoundaryConditions/Zero.hpp

@@ -86,16 +86,17 @@ class Zero : public elliptic::BoundaryConditions::BoundaryCondition<Dim> {
   using argument_tags = tmpl::list<>;
   using volume_tags = tmpl::list<>;
 
-  static void apply(
-      gsl::not_null<tnsr::I<DataVector, Dim>*> displacement,
-      gsl::not_null<tnsr::I<DataVector, Dim>*> n_dot_minus_stress);
+  static void apply(gsl::not_null<tnsr::I<DataVector, Dim>*> displacement,
+                    gsl::not_null<tnsr::I<DataVector, Dim>*> n_dot_minus_stress,
+                    const tnsr::iJ<DataVector, Dim>& deriv_displacement);
 
   using argument_tags_linearized = tmpl::list<>;
   using volume_tags_linearized = tmpl::list<>;
 
   static void apply_linearized(
       gsl::not_null<tnsr::I<DataVector, Dim>*> displacement,
-      gsl::not_null<tnsr::I<DataVector, Dim>*> n_dot_minus_stress);
+      gsl::not_null<tnsr::I<DataVector, Dim>*> n_dot_minus_stress,
+      const tnsr::iJ<DataVector, Dim>& deriv_displacement);
 };
 
 template <size_t Dim, elliptic::BoundaryConditionType BoundaryConditionType>

src/Elliptic/Systems/Poisson/BoundaryConditions/Robin.cpp

@@ -30,7 +30,8 @@ Robin<Dim>::Robin(const double dirichlet_weight, const double neumann_weight,
 template <size_t Dim>
 void Robin<Dim>::apply(
     const gsl::not_null<Scalar<DataVector>*> field,
-    const gsl::not_null<Scalar<DataVector>*> n_dot_field_gradient) const {
+    const gsl::not_null<Scalar<DataVector>*> n_dot_field_gradient,
+    const tnsr::i<DataVector, Dim>& /*deriv_field*/) const {
   if (neumann_weight_ == 0.) {
     ASSERT(
         not equal_within_roundoff(dirichlet_weight_, 0.),
@@ -49,8 +50,8 @@ void Robin<Dim>::apply(
 template <size_t Dim>
 void Robin<Dim>::apply_linearized(
     const gsl::not_null<Scalar<DataVector>*> field_correction,
-    const gsl::not_null<Scalar<DataVector>*> n_dot_field_gradient_correction)
-    const {
+    const gsl::not_null<Scalar<DataVector>*> n_dot_field_gradient_correction,
+    const tnsr::i<DataVector, Dim>& /*deriv_field_correction*/) const {
   if (neumann_weight_ == 0.) {
     get(*field_correction) = 0.;
   } else {

src/Elliptic/Systems/Poisson/BoundaryConditions/Robin.hpp

@@ -90,14 +90,16 @@ class Robin : public elliptic::BoundaryConditions::BoundaryCondition<Dim> {
   using volume_tags = tmpl::list<>;
 
   void apply(gsl::not_null<Scalar<DataVector>*> field,
-             gsl::not_null<Scalar<DataVector>*> n_dot_field_gradient) const;
+             gsl::not_null<Scalar<DataVector>*> n_dot_field_gradient,
+             const tnsr::i<DataVector, Dim>& deriv_field) const;
 
   using argument_tags_linearized = tmpl::list<>;
   using volume_tags_linearized = tmpl::list<>;
 
   void apply_linearized(
       gsl::not_null<Scalar<DataVector>*> field_correction,
-      gsl::not_null<Scalar<DataVector>*> n_dot_field_gradient_correction) const;
+      gsl::not_null<Scalar<DataVector>*> n_dot_field_gradient_correction,
+      const tnsr::i<DataVector, Dim>& deriv_field_correction) const;
 
   void pup(PUP::er& p) override;
 

src/Elliptic/Systems/Punctures/BoundaryConditions/Flatness.cpp

@@ -13,13 +13,15 @@ namespace Punctures::BoundaryConditions {
 void Flatness::apply(
     const gsl::not_null<Scalar<DataVector>*> field,
     const gsl::not_null<Scalar<DataVector>*> n_dot_field_gradient,
+    const tnsr::i<DataVector, 3>& /*field_gradient*/,
     const tnsr::I<DataVector, 3>& x) {
   get(*n_dot_field_gradient) = -get(*field) / get(magnitude(x));
 }
 
 void Flatness::apply_linearized(
     const gsl::not_null<Scalar<DataVector>*> field_correction,
     const gsl::not_null<Scalar<DataVector>*> n_dot_field_gradient_correction,
+    const tnsr::i<DataVector, 3>& /*field_gradient*/,
     const tnsr::I<DataVector, 3>& x) {
   get(*n_dot_field_gradient_correction) =
       -get(*field_correction) / get(magnitude(x));

src/Elliptic/Systems/Punctures/BoundaryConditions/Flatness.hpp

@@ -59,6 +59,7 @@ class Flatness : public elliptic::BoundaryConditions::BoundaryCondition<3> {
 
   static void apply(gsl::not_null<Scalar<DataVector>*> field,
                     gsl::not_null<Scalar<DataVector>*> n_dot_field_gradient,
+                    const tnsr::i<DataVector, 3>& field_gradient,
                     const tnsr::I<DataVector, 3>& x);
 
   using argument_tags_linearized =
@@ -68,6 +69,7 @@ class Flatness : public elliptic::BoundaryConditions::BoundaryCondition<3> {
   static void apply_linearized(
       gsl::not_null<Scalar<DataVector>*> field_correction,
       gsl::not_null<Scalar<DataVector>*> n_dot_field_gradient_correction,
+      const tnsr::i<DataVector, 3>& field_gradient,
       const tnsr::I<DataVector, 3>& x);
 };
 

src/Elliptic/Systems/Xcts/BoundaryConditions/ApparentHorizon.cpp

@@ -506,6 +506,9 @@ void ApparentHorizon<ConformalGeometry>::apply(
         n_dot_lapse_times_conformal_factor_gradient,
     const gsl::not_null<tnsr::I<DataVector, 3>*>
         n_dot_longitudinal_shift_excess,
+    const tnsr::i<DataVector, 3>& /*deriv_conformal_factor*/,
+    const tnsr::i<DataVector, 3>& /*deriv_lapse_times_conformal_factor*/,
+    const tnsr::iJ<DataVector, 3>& /*deriv_shift_excess*/,
     const tnsr::i<DataVector, 3>& face_normal,
     const tnsr::ij<DataVector, 3>& deriv_unnormalized_face_normal,
     const Scalar<DataVector>& face_normal_magnitude,
@@ -535,6 +538,9 @@ void ApparentHorizon<ConformalGeometry>::apply(
         n_dot_lapse_times_conformal_factor_gradient,
     const gsl::not_null<tnsr::I<DataVector, 3>*>
         n_dot_longitudinal_shift_excess,
+    const tnsr::i<DataVector, 3>& /*deriv_conformal_factor*/,
+    const tnsr::i<DataVector, 3>& /*deriv_lapse_times_conformal_factor*/,
+    const tnsr::iJ<DataVector, 3>& /*deriv_shift_excess*/,
     const tnsr::i<DataVector, 3>& face_normal,
     const tnsr::ij<DataVector, 3>& deriv_unnormalized_face_normal,
     const Scalar<DataVector>& face_normal_magnitude,
@@ -568,6 +574,10 @@ void ApparentHorizon<ConformalGeometry>::apply_linearized(
         n_dot_lapse_times_conformal_factor_gradient_correction,
     const gsl::not_null<tnsr::I<DataVector, 3>*>
         n_dot_longitudinal_shift_excess_correction,
+    const tnsr::i<DataVector, 3>& /*deriv_conformal_factor_correction*/,
+    const tnsr::i<DataVector,
+                  3>& /*deriv_lapse_times_conformal_factor_correction*/,
+    const tnsr::iJ<DataVector, 3>& /*deriv_shift_excess_correction*/,
     const tnsr::i<DataVector, 3>& face_normal,
     const tnsr::ij<DataVector, 3>& deriv_unnormalized_face_normal,
     const Scalar<DataVector>& face_normal_magnitude,
@@ -601,6 +611,10 @@ void ApparentHorizon<ConformalGeometry>::apply_linearized(
         n_dot_lapse_times_conformal_factor_gradient_correction,
     const gsl::not_null<tnsr::I<DataVector, 3>*>
         n_dot_longitudinal_shift_excess_correction,
+    const tnsr::i<DataVector, 3>& /*deriv_conformal_factor_correction*/,
+    const tnsr::i<DataVector,
+                  3>& /*deriv_lapse_times_conformal_factor_correction*/,
+    const tnsr::iJ<DataVector, 3>& /*deriv_shift_excess_correction*/,
     const tnsr::i<DataVector, 3>& face_normal,
     const tnsr::ij<DataVector, 3>& deriv_unnormalized_face_normal,
     const Scalar<DataVector>& face_normal_magnitude,