Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include "DataStructures/DataVector.hpp"
       7             : #include "DataStructures/Tensor/Tensor.hpp"
       8             : #include "Utilities/Gsl.hpp"
       9             : 
      10             : namespace Xcts {
      11             : 
      12             : /// @{
      13             : /*!
      14             :  * \brief Surface integrand for the center of mass calculation.
      15             :  *
      16             :  * We define the center of mass integral as
      17             :  *
      18             :  * \begin{equation}
      19             :  *   C_\text{CoM}^i = \frac{3}{8 \pi M_\text{ADM}} \oint_{S_\infty} \Big[
      20             :  *                      (\psi - 1)^4 + 4 (\psi - 1)^3
      21             :  *                      + 6 (\psi - 1)^2 + 4 (\psi - 1)
      22             :  *                    \Big] n^i \, dA,
      23             :  * \end{equation}
      24             :  *
      25             :  * where $n^i = x^i / r$ and $r = \sqrt{x^2 + y^2 + z^2}$. We use $\psi-1$
      26             :  * instead of $\psi$ because it is the variable solved for in the XCTS system.
      27             :  * Expanding the integrand, we see that this is identical to
      28             :  *
      29             :  * \begin{equation}
      30             :  *   C_\text{CoM}^i = \frac{3}{8 \pi M_\text{ADM}}
      31             :  *               \oint_{S_\infty} (\psi^4 - 1) n^i \, dA,
      32             :  * \end{equation}
      33             :  *
      34             :  * Analytically, this is equivalent to the definition in Eq. (25) of
      35             :  * \cite Ossokine2015yla because
      36             :  * \begin{equation}
      37             :  *   \oint_{S_\infty} n^i \, dA = 0.
      38             :  * \end{equation}
      39             :  * Numerically, we have found that subtracting $1$ from $\psi^4$ results in less
      40             :  * round-off errors, leading to a more accurate center of mass.
      41             :  *
      42             :  * One way to interpret this integral is that we are summing over the unit
      43             :  * vectors $n^i$, rescaled by $\psi^4$, in all directions. If $\psi(\vec r)$ is
      44             :  * constant, no rescaling happens and all the unit vectors cancel out. If
      45             :  * $\psi(\vec r)$ is not constant, then $\vec C_\text{CoM}$ will emerge from the
      46             :  * difference of large numbers (sum of rescaled $n^i$ in each subdomain). With
      47             :  * larger and larger numbers being involved in this cancellation (i.e. with
      48             :  * increasing radius of $S_\infty$), we loose numerical accuracy. In other
      49             :  * words, we are seeking the subdominant terms. Since $\psi^4 \to 1$ in
      50             :  * conformal flatness, subtracting $1$ from it in the integrand makes the
      51             :  * numbers involved in this cancellation smaller, reducing this issue.
      52             :  *
      53             :  * \note We don't include the ADM mass $M_{ADM}$ in this integrand. After
      54             :  * integrating the result of this function, you have to divide by $M_{ADM}$.
      55             :  *
      56             :  * \note For consistency with `center_of_mass_volume_integrand`, this
      57             :  * integrand needs to be integrated with the Euclidean area element.
      58             :  *
      59             :  * \see `Xcts::adm_mass_surface_integrand`
      60             :  *
      61             :  * \warning This integral assumes that the conformal metric falls off to
      62             :  * flatness faster than $1/r^2$. That means that it cannot be directly used
      63             :  * with the Kerr-Schild metric, which falls off as $1/r$. This is not a problem
      64             :  * for XCTS with Superposed Kerr-Schild (SKS) because of the exponential
      65             :  * fall-off terms.
      66             :  *
      67             :  * \param result output pointer
      68             :  * \param conformal_factor_minus_one the conformal factor $\psi - 1$
      69             :  * \param coords the inertial coordinates $x^i$
      70             :  */
      71           1 : void center_of_mass_surface_integrand(
      72             :     gsl::not_null<tnsr::I<DataVector, 3>*> result,
      73             :     const Scalar<DataVector>& conformal_factor_minus_one,
      74             :     const tnsr::I<DataVector, 3>& coords);
      75             : 
      76             : /// Return-by-value overload
      77           1 : tnsr::I<DataVector, 3> center_of_mass_surface_integrand(
      78             :     const Scalar<DataVector>& conformal_factor_minus_one,
      79             :     const tnsr::I<DataVector, 3>& coords);
      80             : /// @}
      81             : 
      82             : /// @{
      83             : /*!
      84             :  * \brief Volume integrand for the center of mass calculation.
      85             :  *
      86             :  * We cast the center of mass as an infinite volume integral by applying Gauss'
      87             :  * law on the surface integral defined in `center_of_mass_surface_integrand`:
      88             :  *
      89             :  * \begin{equation}
      90             :  *   C_\text{CoM}^i = \frac{3}{8 \pi M_\text{ADM}}
      91             :  *                    \int_{V_\infty} \Big(
      92             :  *                      4 \psi^3 \partial_j \psi n^i n^j
      93             :  *                      + \frac{2}{r} \psi^4 n^i
      94             :  *                    \Big) dV
      95             :  *                  = \frac{3}{4 \pi M_\text{ADM}}
      96             :  *                    \int_{V_\infty} \frac{1}{r^2} \Big(
      97             :  *                      2 \psi^3 \partial_j \psi x^i x^j
      98             :  *                      + \psi^4 x^i
      99             :  *                    \Big) dV,
     100             :  * \end{equation}
     101             :  *
     102             :  * where $n^i = x^i / r$ and $r = \sqrt{x^2 + y^2 + z^2}$.
     103             :  *
     104             :  * \note For consistency with `center_of_mass_surface_integrand`, this
     105             :  * integrand needs to be integrated with the Euclidean volume element.
     106             :  *
     107             :  * \see `center_of_mass_surface_integrand`
     108             :  *
     109             :  * \warning This integral currently suffers from significant round-off errors.
     110             :  * It is recommended to use only the surface integral if possible.
     111             :  *
     112             :  * \param result output pointer
     113             :  * \param conformal_factor the conformal factor $\psi$
     114             :  * \param deriv_conformal_factor the gradient of the conformal factor
     115             :  * $\partial_i \psi$
     116             :  * \param coords the inertial coordinates $x^i$
     117             :  */
     118           1 : void center_of_mass_volume_integrand(
     119             :     gsl::not_null<tnsr::I<DataVector, 3>*> result,
     120             :     const Scalar<DataVector>& conformal_factor,
     121             :     const tnsr::i<DataVector, 3, Frame::Inertial>& deriv_conformal_factor,
     122             :     const tnsr::I<DataVector, 3>& coords);
     123             : 
     124             : /// Return-by-value overload
     125           1 : tnsr::I<DataVector, 3> center_of_mass_volume_integrand(
     126             :     const Scalar<DataVector>& conformal_factor,
     127             :     const tnsr::i<DataVector, 3, Frame::Inertial>& deriv_conformal_factor,
     128             :     const tnsr::I<DataVector, 3>& coords);
     129             : /// @}
     130             : 
     131             : }  // namespace Xcts