Add Neumann to SuperposedBoostedBinary BoundaryCondition class

for the XCTS system.

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

@@ -12,12 +12,23 @@
 #include <optional>
 
 #include "DataStructures/DataVector.hpp"
+#include "DataStructures/Tags/TempTensor.hpp"
+#include "DataStructures/TempBuffer.hpp"
 #include "DataStructures/Tensor/EagerMath/DeterminantAndInverse.hpp"
 #include "DataStructures/Tensor/IndexType.hpp"
 #include "DataStructures/Tensor/Tensor.hpp"
 #include "DataStructures/Tensor/TypeAliases.hpp"
+#include "Domain/CoordinateMaps/Affine.hpp"
+#include "Domain/CoordinateMaps/CoordinateMap.hpp"
+#include "Domain/CoordinateMaps/CoordinateMap.tpp"
+#include "Domain/CoordinateMaps/ProductMaps.hpp"
 #include "Elliptic/BoundaryConditions/BoundaryConditionType.hpp"
 #include "Elliptic/Systems/Xcts/Tags.hpp"
+#include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
+#include "NumericalAlgorithms/Spectral/Basis.hpp"
+#include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp"
+#include "NumericalAlgorithms/Spectral/Mesh.hpp"
+#include "NumericalAlgorithms/Spectral/Quadrature.hpp"
 #include "PointwiseFunctions/AnalyticSolutions/Xcts/Factory.hpp"
 #include "PointwiseFunctions/AnalyticSolutions/Xcts/Schwarzschild.hpp"
 #include "PointwiseFunctions/GeneralRelativity/Lapse.hpp"
@@ -232,47 +243,90 @@ void implement_apply_neumann(
     const gsl::not_null<tnsr::I<DataVector, 3>*>
         n_dot_longitudinal_shift_excess,
     const std::array<std::optional<std::unique_ptr<IsolatedObjectBase>>, 2>&
-    /*superposed_objects*/,
-    const tnsr::I<DataVector, 3>& /*x*/,
-    const tnsr::i<DataVector, 3>& /*face_normal*/) {
-  /*
-  using analytic_tags = tmpl::list<
-      ::Tags::deriv<Xcts::Tags::ConformalFactorMinusOne<DataVector>,
-                    tmpl::size_t<3>, Frame::Inertial>,
-      ::Tags::deriv<Xcts::Tags::LapseTimesConformalFactorMinusOne<DataVector>,
-                    tmpl::size_t<3>, Frame::Inertial>,
-      ::Tags::deriv<Xcts::Tags::ShiftExcess<DataVector, 3, Frame::Inertial>,
-                    tmpl::size_t<3>, Frame::Inertial>>;
-
-  const auto solution_vars =
-      variables_from_tagged_tuple((*solution)->variables(x, analytic_tags{}));
-  const auto& deriv_conformal_factor_minus_one =
-      get<::Tags::deriv<Xcts::Tags::ConformalFactorMinusOne<DataVector>,
-                        tmpl::size_t<3>, Frame::Inertial>>(solution_vars);
-  const auto& deriv_lapse_times_conformal_factor_minus_one = get<
-      ::Tags::deriv<Xcts::Tags::LapseTimesConformalFactorMinusOne<DataVector>,
-                    tmpl::size_t<3>, Frame::Inertial>>(solution_vars);
-  const auto& deriv_shift_excess =
-      get<::Tags::deriv<Xcts::Tags::ShiftExcess<DataVector, 3, Frame::Inertial>,
-                        tmpl::size_t<3>, Frame::Inertial>>(solution_vars);
-  */
+        superposed_objects,
+    const std::array<double, 2>& xcoords, const std::array<double, 2>& masses,
+    const std::array<double, 3>& momentum_left,
+    const std::array<double, 3>& momentum_right, const double y_offset,
+    const double z_offset, const tnsr::I<DataVector, 3>& x,
+    const tnsr::i<DataVector, 3>& face_normal) {
+  using Affine = domain::CoordinateMaps::Affine;
+  using Affine3D =
+      domain::CoordinateMaps::ProductOf3Maps<Affine, Affine, Affine>;
+  const size_t num_points_1d = 5;
+  const Mesh<3> mesh{num_points_1d, Spectral::Basis::Legendre,
+                     Spectral::Quadrature::GaussLobatto};
+  const auto x_logical = logical_coordinates(mesh);
+
+  TempBuffer<tmpl::list<::Tags::TempScalar<0>, ::Tags::TempScalar<1>,
+                        ::Tags::TempI<2, 3>>>
+      buffer1{num_points_1d * num_points_1d * num_points_1d};
+
+  auto& conformal_factor_minus_one = get<::Tags::TempScalar<0>>(buffer1);
+  auto& lapse_times_conformal_factor_minus_one =
+      get<::Tags::TempScalar<1>>(buffer1);
+  auto& shift_excess = get<::Tags::TempI<2, 3>>(buffer1);
+
+  TempBuffer<tmpl::list<::Tags::Tempi<0, 3>, ::Tags::TempiJ<1, 3>>> buffer2{
+      x.begin()->size()};
+
+  auto& deriv_lapse_times_conformal_factor_minus_one =
+      get<::Tags::Tempi<0, 3>>(buffer2);
+  auto& deriv_shift_excess = get<::Tags::TempiJ<1, 3>>(buffer2);
+
+  // k is running through each point
+  for (size_t k = 0; k < x.begin()->size(); ++k) {
+    const std::array<double, 3> lower_bound{
+        {x.get(0)[k] - 0.001, x.get(1)[k] - 0.001, x.get(2)[k] - 0.001}};
+    const std::array<double, 3> upper_bound{
+        {x.get(0)[k] + 0.001, x.get(1)[k] + 0.001, x.get(2)[k] + 0.001}};
+    const auto coord_map =
+        domain::make_coordinate_map<Frame::ElementLogical, Frame::Inertial>(
+            Affine3D{Affine{-1., 1., lower_bound[0], upper_bound[0]},
+                     Affine{-1., 1., lower_bound[1], upper_bound[1]},
+                     Affine{-1., 1., lower_bound[2], upper_bound[2]}});
+    const auto x_in = coord_map(x_logical);
+    const auto inv_jacobian = coord_map.inv_jacobian(x_logical);
+
+    implement_apply_dirichlet<IsolatedObjectBase, IsolatedObjectClasses>(
+        make_not_null(&conformal_factor_minus_one),
+        make_not_null(&lapse_times_conformal_factor_minus_one),
+        make_not_null(&shift_excess), superposed_objects, xcoords, masses,
+        momentum_left, momentum_right, y_offset, z_offset, x_in);
+
+    auto deriv_lapse_times_conformal_factor_minus_one_in = partial_derivative(
+        lapse_times_conformal_factor_minus_one, mesh, inv_jacobian);
+    auto deriv_shift_excess_in =
+        partial_derivative(shift_excess, mesh, inv_jacobian);
+
+    for (size_t l = 0; l < num_points_1d * num_points_1d * num_points_1d; ++l) {
+      if (x_in.get(0)[l] == x.get(0)[k] && x_in.get(1)[l] == x.get(1)[k] &&
+          x_in.get(2)[l] == x.get(2)[k]) {
+        for (size_t i = 0; i < 3; ++i) {
+          deriv_lapse_times_conformal_factor_minus_one.get(i)[k] =
+              deriv_lapse_times_conformal_factor_minus_one_in.get(i)[l];
+          for (size_t j = 0; j < 3; ++j) {
+            deriv_shift_excess.get(i, j)[k] =
+                deriv_shift_excess_in.get(i, j)[l];
+          }
+        }
+      }
+    }
+  }
+
   get(*n_dot_conformal_factor_gradient) = 0.;
   get(*n_dot_lapse_times_conformal_factor_gradient) = 0.;
   std::fill(n_dot_longitudinal_shift_excess->begin(),
             n_dot_longitudinal_shift_excess->end(), 0.);
-  /*
+
   for (size_t i = 0; i < 3; ++i) {
-    get(*n_dot_conformal_factor_gradient) +=
-        face_normal.get(i) * deriv_conformal_factor_minus_one.get(i);
     get(*n_dot_lapse_times_conformal_factor_gradient) +=
         face_normal.get(i) *
         deriv_lapse_times_conformal_factor_minus_one.get(i);
     for (size_t j = 0; j < 3; ++j) {
       n_dot_longitudinal_shift_excess->get(i) +=
-          face_normal.get(j) * deriv_shift_excess.get(i, j);
+          face_normal.get(j) * deriv_shift_excess.get(j, i);
     }
   }
-  */
 }
 
 }  // namespace
@@ -303,7 +357,8 @@ void SuperposedBoostedBinary<IsolatedObjectBase, IsolatedObjectClasses>::apply(
     implement_apply_neumann<IsolatedObjectBase, IsolatedObjectClasses>(
         n_dot_conformal_factor_gradient,
         n_dot_lapse_times_conformal_factor_gradient,
-        n_dot_longitudinal_shift_excess, superposed_objects_, x, face_normal);
+        n_dot_longitudinal_shift_excess, superposed_objects_, xcoords_, masses_,
+        momentum_left_, momentum_right_, y_offset_, z_offset_, x, face_normal);
   }
 }
 

tests/Unit/Elliptic/Systems/Xcts/BoundaryConditions/SuperposedBoostedBinary.py

@@ -166,8 +166,65 @@ def n_dot_conformal_factor_gradient(x, face_normal):
 
 
 def n_dot_lapse_times_conformal_factor_gradient(x, face_normal):
-    return 0.0
+    h = 4e-4
+    derivatives = np.zeros(3)
+    coeffs = np.array(
+        [1 / 280, -4 / 105, 1 / 5, -4 / 5, 0, 4 / 5, -1 / 5, 4 / 105, -1 / 280]
+    )
+
+    f_values_x0 = np.array(
+        [
+            lapse_times_conformal_factor_minus_one([x[0] + i * h, x[1], x[2]])
+            for i in range(-4, 5)
+        ]
+    )
+    derivatives[0] = np.dot(coeffs, f_values_x0) / h
+
+    f_values_x1 = np.array(
+        [
+            lapse_times_conformal_factor_minus_one([x[0], x[1] + i * h, x[2]])
+            for i in range(-4, 5)
+        ]
+    )
+    derivatives[1] = np.dot(coeffs, f_values_x1) / h
+
+    f_values_x2 = np.array(
+        [
+            lapse_times_conformal_factor_minus_one([x[0], x[1], x[2] + i * h])
+            for i in range(-4, 5)
+        ]
+    )
+    derivatives[2] = np.dot(coeffs, f_values_x2) / h
+
+    return np.dot(derivatives, face_normal)
 
 
 def n_dot_longitudinal_shift_excess(x, face_normal):
-    return np.zeros(3)
+    h = 4e-4
+    derivatives = np.zeros((3, 3))
+    result = np.zeros(3)
+    coeffs = np.array(
+        [1 / 280, -4 / 105, 1 / 5, -4 / 5, 0, 4 / 5, -1 / 5, 4 / 105, -1 / 280]
+    )
+
+    for i in range(3):
+        f_values_x0 = np.array(
+            [shift_excess([x[0] + j * h, x[1], x[2]])[i] for j in range(-4, 5)]
+        )
+        derivatives[i, 0] = np.dot(coeffs, f_values_x0) / h
+
+        f_values_x1 = np.array(
+            [shift_excess([x[0], x[1] + j * h, x[2]])[i] for j in range(-4, 5)]
+        )
+        derivatives[i, 1] = np.dot(coeffs, f_values_x1) / h
+
+        f_values_x2 = np.array(
+            [shift_excess([x[0], x[1], x[2] + j * h])[i] for j in range(-4, 5)]
+        )
+        derivatives[i, 2] = np.dot(coeffs, f_values_x2) / h
+
+    for i in range(3):
+        for j in range(3):
+            result[i] += derivatives[i, j] * face_normal[j]
+
+    return result

tests/Unit/Elliptic/Systems/Xcts/BoundaryConditions/Test_SuperposedBoostedBinary.cpp

@@ -95,7 +95,7 @@ void test_suite(const elliptic::BoundaryConditionType& boundary,
   }
   {
     MAKE_GENERATOR(gen);
-    std::uniform_real_distribution<> dist(-1., 1.);
+    std::uniform_real_distribution<> dist(-10., -8.);
     const size_t num_points = 3;
     const auto direction = Direction<3>::upper_zeta();
     const std::array<double, 3> center{{0., 0., 0.}};
@@ -131,29 +131,6 @@ void test_suite(const elliptic::BoundaryConditionType& boundary,
         num_points, std::numeric_limits<double>::signaling_NaN()};
     const tnsr::iJ<DataVector, 3> deriv_shift_excess{
         num_points, std::numeric_limits<double>::signaling_NaN()};
-    /*
-    if constexpr (Linearized == false) {
-      boundary_condition.apply(
-          make_not_null(&conformal_factor_minus_one),
-          make_not_null(&lapse_times_conformal_factor_minus_one),
-          make_not_null(&shift_excess),
-          make_not_null(&n_dot_conformal_factor_gradient),
-          make_not_null(&n_dot_lapse_times_conformal_factor_gradient),
-          make_not_null(&n_dot_longitudinal_shift_excess),
-          deriv_conformal_factor, deriv_lapse_times_conformal_factor,
-          deriv_shift_excess, x, face_normal);
-    } else if constexpr (Linearized == true) {
-      boundary_condition.apply_linearized(
-          make_not_null(&conformal_factor_minus_one),
-          make_not_null(&lapse_times_conformal_factor_minus_one),
-          make_not_null(&shift_excess),
-          make_not_null(&n_dot_conformal_factor_gradient),
-          make_not_null(&n_dot_lapse_times_conformal_factor_gradient),
-          make_not_null(&n_dot_longitudinal_shift_excess),
-          deriv_conformal_factor, deriv_lapse_times_conformal_factor,
-          deriv_shift_excess, x, face_normal);
-    }
-    */
 
     elliptic::apply_boundary_condition<
         Linearized, void,
@@ -200,18 +177,26 @@ void test_suite(const elliptic::BoundaryConditionType& boundary,
       const auto expected_n_dot_longitudinal_shift_excess =
           pypp::call<tnsr::I<DataVector, 3>>(
               py_module, "n_dot_longitudinal_shift_excess", x, face_normal);
-
+      // clang-tidy: static object creation may throw exception
+      static Approx approx =                                     // NOLINT
+          Approx::custom()                                       // NOLINT
+              .epsilon(std::numeric_limits<double>::epsilon() *  // NOLINT
+                       10000)                                    // NOLINT
+              .scale(1.0);                                       // NOLINT
       CHECK_ITERABLE_APPROX(get(n_dot_conformal_factor_gradient),
                             get(expected_n_dot_conformal_factor_gradient));
-      CHECK_ITERABLE_APPROX(
+      CHECK_ITERABLE_CUSTOM_APPROX(
           get(n_dot_lapse_times_conformal_factor_gradient),
-          get(expected_n_dot_lapse_times_conformal_factor_gradient));
-      CHECK_ITERABLE_APPROX(get<0>(n_dot_longitudinal_shift_excess),
-                            get<0>(expected_n_dot_longitudinal_shift_excess));
-      CHECK_ITERABLE_APPROX(get<1>(n_dot_longitudinal_shift_excess),
-                            get<1>(expected_n_dot_longitudinal_shift_excess));
-      CHECK_ITERABLE_APPROX(get<2>(n_dot_longitudinal_shift_excess),
-                            get<2>(expected_n_dot_longitudinal_shift_excess));
+          get(expected_n_dot_lapse_times_conformal_factor_gradient), approx);
+      CHECK_ITERABLE_CUSTOM_APPROX(
+          get<0>(n_dot_longitudinal_shift_excess),
+          get<0>(expected_n_dot_longitudinal_shift_excess), approx);
+      CHECK_ITERABLE_CUSTOM_APPROX(
+          get<1>(n_dot_longitudinal_shift_excess),
+          get<1>(expected_n_dot_longitudinal_shift_excess), approx);
+      CHECK_ITERABLE_CUSTOM_APPROX(
+          get<2>(n_dot_longitudinal_shift_excess),
+          get<2>(expected_n_dot_longitudinal_shift_excess), approx);
     }
   }
 }