fixup: change CoM parameters.

sxs-collaboration · Sep 16, 2024 · 14d7a3a · 14d7a3a
1 parent a487791
commit 14d7a3a
Show file tree

Hide file tree

Showing 3 changed files with 27 additions and 12 deletions.
diff --git a/src/PointwiseFunctions/Xcts/CenterOfMass.cpp b/src/PointwiseFunctions/Xcts/CenterOfMass.cpp
@@ -3,6 +3,8 @@
 
 #include "PointwiseFunctions/Xcts/CenterOfMass.hpp"
 
+#include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
+
 namespace Xcts {
 
 void center_of_mass_surface_integrand(
@@ -27,20 +29,22 @@ void center_of_mass_volume_integrand(
     gsl::not_null<tnsr::I<DataVector, 3>*> result,
     const Scalar<DataVector>& conformal_factor,
     const tnsr::i<DataVector, 3, Frame::Inertial>& deriv_conformal_factor,
-    const tnsr::I<DataVector, 3>& unit_normal) {
-  tenex::evaluate<ti::I>(result, 3. / (8. * M_PI) * 4. *
-                                     pow<3>(conformal_factor()) *
-                                     deriv_conformal_factor(ti::j) *
-                                     unit_normal(ti::I) * unit_normal(ti::J));
+    const tnsr::I<DataVector, 3>& coords) {
+  tenex::evaluate<ti::I>(
+      result,
+      3. / (4. * M_PI * pow<2>(magnitude(coords)())) *
+          (2. * pow<3>(conformal_factor()) * deriv_conformal_factor(ti::j) *
+               coords(ti::I) * coords(ti::J) +
+           pow<4>(conformal_factor()) * coords(ti::I)));
 }
 
 tnsr::I<DataVector, 3> center_of_mass_volume_integrand(
     const Scalar<DataVector>& conformal_factor,
     const tnsr::i<DataVector, 3, Frame::Inertial>& deriv_conformal_factor,
-    const tnsr::I<DataVector, 3>& unit_normal) {
+    const tnsr::I<DataVector, 3>& coords) {
   tnsr::I<DataVector, 3> result;
   center_of_mass_volume_integrand(make_not_null(&result), conformal_factor,
-                                  deriv_conformal_factor, unit_normal);
+                                  deriv_conformal_factor, coords);
   return result;
 }
 

diff --git a/src/PointwiseFunctions/Xcts/CenterOfMass.hpp b/src/PointwiseFunctions/Xcts/CenterOfMass.hpp
@@ -56,28 +56,38 @@ tnsr::II<DataVector, 3> center_of_mass_surface_integrand(
  *
  * \begin{equation}
  *   C_\text{CoM}^i = \frac{3}{8 \pi M_\text{ADM}}
- *               \int_{V_\infty} 4 \psi^3 \partial_j \psi n^i n^j \, dV.
+ *                    \int_{V_\infty} \Big(
+ *                      4 \psi^3 \partial_j \psi n^i n^j
+ *                      + \frac{2}{r} \psi^4 n^i
+ *                    \Big) dV
+ *                  = \frac{3}{4 \pi M_\text{ADM}}
+ *                    \int_{V_\infty} \frac{1}{r^2} \Big(
+ *                      2 \psi^3 \partial_j \psi x^i x^j
+ *                      + \psi^4 x^i
+ *                    \Big) dV,
  * \end{equation}
  *
+ * where $n^i = x^i / r$ and $r = \sqrt{x^2 + y^2 + z^2}$.
+ *
  * \see `center_of_mass_surface_integrand`
  *
  * \param result output pointer
  * \param conformal_factor the conformal factor $\psi$
  * \param deriv_conformal_factor the gradient of the conformal factor
  * $\partial_i \psi$
- * \param unit_normal the outward-pointing unit normal $n^i = x^i / r$
+ * \param coords the inertial coordinates $x^i$
  */
 void center_of_mass_volume_integrand(
     gsl::not_null<tnsr::I<DataVector, 3>*> result,
     const Scalar<DataVector>& conformal_factor,
     const tnsr::i<DataVector, 3, Frame::Inertial>& deriv_conformal_factor,
-    const tnsr::I<DataVector, 3>& unit_normal);
+    const tnsr::I<DataVector, 3>& coords);
 
 /// Return-by-value overload
 tnsr::I<DataVector, 3> center_of_mass_volume_integrand(
     const Scalar<DataVector>& conformal_factor,
     const tnsr::i<DataVector, 3, Frame::Inertial>& deriv_conformal_factor,
-    const tnsr::I<DataVector, 3>& unit_normal);
+    const tnsr::I<DataVector, 3>& coords);
 /// @}
 
 }  // namespace Xcts
diff --git a/tests/Unit/PointwiseFunctions/Xcts/Test_CenterOfMass.cpp b/tests/Unit/PointwiseFunctions/Xcts/Test_CenterOfMass.cpp
@@ -235,7 +235,7 @@ void test_infinite_volume_integral(const double distance,
 
     // Evaluate volume integral.
     const auto volume_integrand = Xcts::center_of_mass_volume_integrand(
-        conformal_factor, deriv_conformal_factor, unit_normal);
+        conformal_factor, deriv_conformal_factor, inertial_coords);
     for (int I = 0; I < 3; I++) {
       total_integral.get(I) +=
           definite_integral(volume_integrand.get(I) * get(det_jacobian), mesh);
@@ -302,6 +302,7 @@ SPECTRE_TEST_CASE("Unit.PointwiseFunctions.Xcts.CenterOfMass",
   for (const double distance : std::array<double, 3>({1.e4, 1.e5, 1.e6})) {
     test_infinite_surface_integral(distance, 0.1);
     // Note: the infinite volume integral blows up for non-zero shifts.
+    // TODO: figure out why the test below fails with the recent changes.
     test_infinite_volume_integral(distance, 0.);
   }
 }