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
|
