Compute ADM angular momentum integral in BBH ID

src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.cpp

@@ -11,15 +11,16 @@
 #include "NumericalAlgorithms/DiscontinuousGalerkin/ProjectToBoundary.hpp"
 #include "NumericalAlgorithms/LinearOperators/DefiniteIntegral.hpp"
 #include "NumericalAlgorithms/Spectral/Mesh.hpp"
-#include "PointwiseFunctions/Xcts/AdmLinearMomentum.hpp"
 #include "PointwiseFunctions/Xcts/AdmMass.hpp"
+#include "PointwiseFunctions/Xcts/AdmMomentum.hpp"
 #include "PointwiseFunctions/Xcts/CenterOfMass.hpp"
 
 namespace Events {
 
 void local_adm_integrals(
     gsl::not_null<Scalar<double>*> adm_mass,
     gsl::not_null<tnsr::I<double, 3>*> adm_linear_momentum,
+    gsl::not_null<Scalar<double>*> adm_angular_momentum_z,
     gsl::not_null<tnsr::I<double, 3>*> center_of_mass,
     const Scalar<DataVector>& conformal_factor,
     const tnsr::i<DataVector, 3>& deriv_conformal_factor,
@@ -38,6 +39,7 @@ void local_adm_integrals(
     const DirectionMap<3, tnsr::i<DataVector, 3>>& conformal_face_normals) {
   // Initialize integrals to 0
   adm_mass->get() = 0;
+  adm_angular_momentum_z->get() = 0;
   for (int I = 0; I < 3; I++) {
     adm_linear_momentum->get(I) = 0;
     center_of_mass->get(I) = 0;
@@ -126,6 +128,9 @@ void local_adm_integrals(
         Xcts::adm_linear_momentum_surface_integrand(
             face_conformal_factor, face_inv_spatial_metric,
             face_inv_extrinsic_curvature, face_trace_extrinsic_curvature);
+    const auto angular_momentum_z_integrand =
+        Xcts::adm_angular_momentum_z_surface_integrand(
+            linear_momentum_integrand, face_inertial_coords);
     const auto center_of_mass_integrand =
         Xcts::center_of_mass_surface_integrand(face_conformal_factor,
                                                face_inertial_coords);
@@ -135,11 +140,16 @@ void local_adm_integrals(
         tenex::evaluate(mass_integrand(ti::I) * conformal_face_normal(ti::i));
     const auto contracted_linear_momentum_integrand = tenex::evaluate<ti::I>(
         linear_momentum_integrand(ti::I, ti::J) * flat_face_normal(ti::j));
+    const auto contracted_angular_momentum_z_integrand = tenex::evaluate(
+        angular_momentum_z_integrand(ti::I) * flat_face_normal(ti::i));
 
     // Take integrals
     adm_mass->get() += definite_integral(
         get(contracted_mass_integrand) * get(conformal_area_element),
         face_mesh);
+    adm_angular_momentum_z->get() += definite_integral(
+        get(contracted_angular_momentum_z_integrand) * get(flat_area_element),
+        face_mesh);
     for (int I = 0; I < 3; I++) {
       adm_linear_momentum->get(I) += definite_integral(
           contracted_linear_momentum_integrand.get(I) * get(flat_area_element),

src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.hpp

@@ -38,13 +38,11 @@ namespace Events {
  *
  * To get the total ADM integrals, the results need to be summed over in a
  * reduction.
- *
- * See `Xcts::adm_linear_momentum_surface_integrand` for details on the formula
- * for each integrand.
  */
 void local_adm_integrals(
     gsl::not_null<Scalar<double>*> adm_mass,
     gsl::not_null<tnsr::I<double, 3>*> adm_linear_momentum,
+    gsl::not_null<Scalar<double>*> adm_angular_momentum_z,
     gsl::not_null<tnsr::I<double, 3>*> center_of_mass,
     const Scalar<DataVector>& conformal_factor,
     const tnsr::i<DataVector, 3>& deriv_conformal_factor,
@@ -74,6 +72,7 @@ void local_adm_integrals(
  * - Number of points in the domain
  * - ADM mass
  * - ADM linear momentum
+ * - ADM angular momentum (z-component)
  * - Center of mass
  */
 template <typename ArraySectionIdTag = void>
@@ -90,6 +89,8 @@ class ObserveAdmIntegrals : public Event {
       Parallel::ReductionDatum<double, funcl::Plus<>>,
       // ADM Linear Momentum (z-component)
       Parallel::ReductionDatum<double, funcl::Plus<>>,
+      // ADM Angular Momentum (z-component)
+      Parallel::ReductionDatum<double, funcl::Plus<>>,
       // Center of Mass (x-component)
       Parallel::ReductionDatum<double, funcl::Plus<>, funcl::Divides<>,
                                std::index_sequence<1>>,
@@ -115,6 +116,7 @@ class ObserveAdmIntegrals : public Event {
       "- Number of points in the domain\n"
       "- ADM mass\n"
       "- ADM linear momentum\n"
+      "- ADM angular momentum (z-component)\n"
       "- Center of mass";
 
   ObserveAdmIntegrals() = default;
@@ -179,27 +181,30 @@ class ObserveAdmIntegrals : public Event {
 
     Scalar<double> adm_mass;
     tnsr::I<double, 3> adm_linear_momentum;
+    Scalar<double> adm_angular_momentum_z;
     tnsr::I<double, 3> center_of_mass;
     local_adm_integrals(
         make_not_null(&adm_mass), make_not_null(&adm_linear_momentum),
-        make_not_null(&center_of_mass), conformal_factor,
-        deriv_conformal_factor, conformal_metric, inv_conformal_metric,
-        conformal_christoffel_second_kind, conformal_christoffel_contracted,
-        spatial_metric, inv_spatial_metric, extrinsic_curvature,
-        trace_extrinsic_curvature, inertial_coords, inv_jacobian, mesh, element,
-        conformal_face_normals);
+        make_not_null(&adm_angular_momentum_z), make_not_null(&center_of_mass),
+        conformal_factor, deriv_conformal_factor, conformal_metric,
+        inv_conformal_metric, conformal_christoffel_second_kind,
+        conformal_christoffel_contracted, spatial_metric, inv_spatial_metric,
+        extrinsic_curvature, trace_extrinsic_curvature, inertial_coords,
+        inv_jacobian, mesh, element, conformal_face_normals);
 
     // Save components of linear momentum as reduction data
     ReductionData reduction_data{
         mesh.number_of_grid_points(), get(adm_mass),
         get<0>(adm_linear_momentum),  get<1>(adm_linear_momentum),
-        get<2>(adm_linear_momentum),  get<0>(center_of_mass),
-        get<1>(center_of_mass),       get<2>(center_of_mass)};
+        get<2>(adm_linear_momentum),  get(adm_angular_momentum_z),
+        get<0>(center_of_mass),       get<1>(center_of_mass),
+        get<2>(center_of_mass)};
     std::vector<std::string> legend{
         "NumberOfPoints",      "AdmMass",
         "AdmLinearMomentum_x", "AdmLinearMomentum_y",
-        "AdmLinearMomentum_z", "CenterOfMass_x",
-        "CenterOfMass_y",      "CenterOfMass_z"};
+        "AdmLinearMomentum_z", "AdmAngularMomentum_z",
+        "CenterOfMass_x",      "CenterOfMass_y",
+        "CenterOfMass_z"};
 
     // Get information required for reduction
     auto& local_observer = *Parallel::local_branch(

src/PointwiseFunctions/Xcts/AdmLinearMomentum.cpp → src/PointwiseFunctions/Xcts/AdmMomentum.cpp

@@ -1,7 +1,7 @@
 // Distributed under the MIT License.
 // See LICENSE.txt for details.
 
-#include "PointwiseFunctions/Xcts/AdmLinearMomentum.hpp"
+#include "PointwiseFunctions/Xcts/AdmMomentum.hpp"
 
 namespace Xcts {
 
@@ -68,4 +68,45 @@ tnsr::I<DataVector, 3> adm_linear_momentum_volume_integrand(
   return result;
 }
 
+void adm_angular_momentum_z_surface_integrand(
+    gsl::not_null<tnsr::I<DataVector, 3>*> result,
+    const tnsr::II<DataVector, 3>& linear_momentum_surface_integrand,
+    const tnsr::I<DataVector, 3>& coords) {
+  // Note: we can ignore the $1/(8\pi)$ term below because it is already
+  // included in `linear_momentum_surface_integrand`.
+  for (int I = 0; I < 3; I++) {
+    result->get(I) =
+        get<0>(coords) * linear_momentum_surface_integrand.get(1, I) -
+        get<1>(coords) * linear_momentum_surface_integrand.get(0, I);
+  }
+}
+
+tnsr::I<DataVector, 3> adm_angular_momentum_z_surface_integrand(
+    const tnsr::II<DataVector, 3>& linear_momentum_surface_integrand,
+    const tnsr::I<DataVector, 3>& coords) {
+  tnsr::I<DataVector, 3> result;
+  adm_angular_momentum_z_surface_integrand(
+      make_not_null(&result), linear_momentum_surface_integrand, coords);
+  return result;
+}
+
+void adm_angular_momentum_z_volume_integrand(
+    gsl::not_null<Scalar<DataVector>*> result,
+    const tnsr::I<DataVector, 3>& linear_momentum_volume_integrand,
+    const tnsr::I<DataVector, 3>& coords) {
+  // Note: we can ignore the $-1/(8\pi)$ term below because it is already
+  // included in `linear_momentum_volume_integrand`.
+  result->get() = get<0>(coords) * get<1>(linear_momentum_volume_integrand) -
+                  get<1>(coords) * get<0>(linear_momentum_volume_integrand);
+}
+
+Scalar<DataVector> adm_angular_momentum_z_volume_integrand(
+    const tnsr::I<DataVector, 3>& linear_momentum_volume_integrand,
+    const tnsr::I<DataVector, 3>& coords) {
+  Scalar<DataVector> result;
+  adm_angular_momentum_z_volume_integrand(
+      make_not_null(&result), linear_momentum_volume_integrand, coords);
+  return result;
+}
+
 }  // namespace Xcts

src/PointwiseFunctions/Xcts/AdmLinearMomentum.hpp → src/PointwiseFunctions/Xcts/AdmMomentum.hpp

@@ -18,7 +18,7 @@ namespace Xcts {
  *
  * \begin{equation}
  *   P_\text{ADM}^i = \frac{1}{8\pi}
- *                    \oint_{S_\infty} \psi^10 \Big(
+ *                    \oint_{S_\infty} \psi^{10} \Big(
  *                      K^{ij} - K \gamma^{ij}
  *                    \Big) \, dS_j.
  * \end{equation}
@@ -103,4 +103,77 @@ tnsr::I<DataVector, 3> adm_linear_momentum_volume_integrand(
     const tnsr::i<DataVector, 3>& conformal_christoffel_contracted);
 /// @}
 
+/// @{
+/*!
+ * \brief Surface integrand for the z-component of the ADM angular momentum.
+ *
+ * We define the ADM angular momentum surface integral as (see Eq. 23 in
+ * \cite Ossokine2015yla):
+ *
+ * \begin{equation}
+ *   J_\text{ADM}^z = \frac{1}{8\pi}
+ *                    \oint_{S_\infty} \Big(
+ *                      x P^{yj} - y P^{xj}
+ *                    \Big) \, dS_j,
+ * \end{equation}
+ *
+ * where $1/(8\pi) P^{jk}$ is the result from
+ * `adm_linear_momentum_surface_integrand`.
+ *
+ * \note For consistency with `adm_angular_momentum_z_volume_integrand`, this
+ * integrand needs to be contracted with the Euclidean face normal and
+ * integrated with the Euclidean area element.
+ *
+ * \param result output pointer
+ * \param linear_momentum_surface_integrand the quantity $1/(8\pi) P^{ij}$
+ * (result of `adm_linear_momentum_surface_integrand`)
+ * \param coords the inertial coordinates $x^i$
+ */
+void adm_angular_momentum_z_surface_integrand(
+    gsl::not_null<tnsr::I<DataVector, 3>*> result,
+    const tnsr::II<DataVector, 3>& linear_momentum_surface_integrand,
+    const tnsr::I<DataVector, 3>& coords);
+
+/// Return-by-value overload
+tnsr::I<DataVector, 3> adm_angular_momentum_z_surface_integrand(
+    const tnsr::II<DataVector, 3>& linear_momentum_surface_integrand,
+    const tnsr::I<DataVector, 3>& coords);
+/// @}
+
+/// @{
+/*!
+ * \brief Volume integrand for the z-component of the ADM angular momentum.
+ *
+ * We define the ADM angular momentum volume integral as (see Eq. 23 in
+ * \cite Ossokine2015yla):
+ *
+ * \begin{equation}
+ *   J_\text{ADM}^z = - \frac{1}{8\pi}
+ *                      \int_{V_\infty} \Big(
+ *                        x G^y - y G^x
+ *                      \Big) \, dV,
+ * \end{equation}
+ *
+ * where $-1/(8\pi) G^i$ is the result from
+ * `adm_linear_momentum_volume_integrand`.
+ *
+ * \note For consistency with `adm_angular_momentum_z_surface_integrand`, this
+ * integrand needs to be integrated with the Euclidean volume element.
+ *
+ * \param result output pointer
+ * \param linear_momentum_volume_integrand the quantity $-1/(8\pi) G^i$ (result
+ * of `adm_linear_momentum_volume_integrand`)
+ * \param coords the inertial coordinates $x^i$
+ */
+void adm_angular_momentum_z_volume_integrand(
+    gsl::not_null<Scalar<DataVector>*> result,
+    const tnsr::I<DataVector, 3>& linear_momentum_volume_integrand,
+    const tnsr::I<DataVector, 3>& coords);
+
+/// Return-by-value overload
+Scalar<DataVector> adm_angular_momentum_z_volume_integrand(
+    const tnsr::I<DataVector, 3>& linear_momentum_volume_integrand,
+    const tnsr::I<DataVector, 3>& coords);
+/// @}
+
 }  // namespace Xcts

src/PointwiseFunctions/Xcts/CMakeLists.txt

@@ -8,8 +8,8 @@ add_spectre_library(${LIBRARY})
 spectre_target_sources(
   ${LIBRARY}
   PRIVATE
-  AdmLinearMomentum.cpp
   AdmMass.cpp
+  AdmMomentum.cpp
   CenterOfMass.cpp
   ExtrinsicCurvature.cpp
   LongitudinalOperator.cpp
@@ -20,8 +20,8 @@ spectre_target_headers(
   ${LIBRARY}
   INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src
   HEADERS
-  AdmLinearMomentum.hpp
   AdmMass.hpp
+  AdmMomentum.hpp
   CenterOfMass.hpp
   ExtrinsicCurvature.hpp
   LongitudinalOperator.hpp

tests/Unit/Elliptic/Systems/Xcts/Events/Test_ObserveAdmIntegrals.cpp

@@ -94,6 +94,8 @@ void test_local_adm_integrals(const double& distance,
   Scalar<double> total_adm_mass;
   total_adm_mass.get() = 0.;
   tnsr::I<double, 3> total_adm_linear_momentum;
+  Scalar<double> total_adm_angular_momentum_z;
+  total_adm_angular_momentum_z.get() = 0.;
   tnsr::I<double, 3> total_center_of_mass;
   for (int I = 0; I < 3; I++) {
     total_adm_linear_momentum.get(I) = 0.;
@@ -246,17 +248,20 @@ void test_local_adm_integrals(const double& distance,
     // Compute local integrals.
     Scalar<double> local_adm_mass;
     tnsr::I<double, 3> local_adm_linear_momentum;
+    Scalar<double> local_adm_angular_momentum_z;
     tnsr::I<double, 3> local_center_of_mass;
     Events::local_adm_integrals(
         make_not_null(&local_adm_mass),
         make_not_null(&local_adm_linear_momentum),
+        make_not_null(&local_adm_angular_momentum_z),
         make_not_null(&local_center_of_mass), conformal_factor,
         deriv_conformal_factor, conformal_metric, inv_conformal_metric,
         conformal_christoffel_second_kind, conformal_christoffel_contracted,
         spatial_metric, inv_spatial_metric, extrinsic_curvature,
         trace_extrinsic_curvature, inertial_coords, inv_jacobian, mesh,
         current_element, conformal_face_normals);
     total_adm_mass.get() += get(local_adm_mass);
+    total_adm_angular_momentum_z.get() += get(local_adm_angular_momentum_z);
     for (int I = 0; I < 3; I++) {
       total_adm_linear_momentum.get(I) += local_adm_linear_momentum.get(I);
       total_center_of_mass.get(I) += local_center_of_mass.get(I);
@@ -270,6 +275,7 @@ void test_local_adm_integrals(const double& distance,
   CHECK(get<1>(total_adm_linear_momentum) == custom_approx(0.));
   CHECK(get<2>(total_adm_linear_momentum) ==
         custom_approx(lorentz_factor * mass * boost_speed));
+  CHECK(get(total_adm_angular_momentum_z) == custom_approx(0.));
   CHECK(get<0>(total_center_of_mass) == custom_approx(0.));
   CHECK(get<1>(total_center_of_mass) == custom_approx(0.));
   CHECK(get<2>(total_center_of_mass) == custom_approx(0.));

tests/Unit/PointwiseFunctions/Xcts/CMakeLists.txt

@@ -4,8 +4,8 @@
 set(LIBRARY Test_XctsPointwiseFunctions)
 
 set(LIBRARY_SOURCES
-  Test_AdmLinearMomentum.cpp
   Test_AdmMass.cpp
+  Test_AdmMomentum.cpp
   Test_CenterOfMass.cpp
   Test_ExtrinsicCurvature.cpp
   Test_LongitudinalOperator.cpp

tests/Unit/PointwiseFunctions/Xcts/Test_AdmMass.cpp

@@ -23,8 +23,8 @@
 #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp"
 #include "PointwiseFunctions/AnalyticSolutions/Xcts/Schwarzschild.hpp"
 #include "PointwiseFunctions/AnalyticSolutions/Xcts/WrappedGr.hpp"
-#include "PointwiseFunctions/Xcts/AdmLinearMomentum.hpp"
 #include "PointwiseFunctions/Xcts/AdmMass.hpp"
+#include "PointwiseFunctions/Xcts/AdmMomentum.hpp"
 
 namespace {