From 07af5523c36e057c55fd73106c0c92e057f6deb8 Mon Sep 17 00:00:00 2001 From: Iago Mendes Date: Thu, 5 Sep 2024 21:26:43 -0700 Subject: [PATCH] Add volume integrands for XCTS asymptotic quantities. --- .../Xcts/Events/ObserveAdmIntegrals.cpp | 9 +- .../Xcts/AdmLinearMomentum.cpp | 23 +- .../Xcts/AdmLinearMomentum.hpp | 50 ++-- src/PointwiseFunctions/Xcts/AdmMass.cpp | 52 ++++ src/PointwiseFunctions/Xcts/AdmMass.hpp | 92 ++++++- src/PointwiseFunctions/Xcts/CenterOfMass.cpp | 32 ++- src/PointwiseFunctions/Xcts/CenterOfMass.hpp | 47 +++- .../Xcts/Events/Test_ObserveAdmIntegrals.cpp | 27 +- .../Xcts/Test_AdmLinearMomentum.cpp | 250 ++++++++++++++--- .../PointwiseFunctions/Xcts/Test_AdmMass.cpp | 251 ++++++++++++++++-- .../Xcts/Test_CenterOfMass.cpp | 184 +++++++++++-- 11 files changed, 872 insertions(+), 145 deletions(-) diff --git a/src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.cpp b/src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.cpp index 37309ad15105..7537cc69c83f 100644 --- a/src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.cpp +++ b/src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.cpp @@ -131,6 +131,8 @@ void local_adm_integrals( tenex::evaluate(mass_integrand(ti::I) * conformal_face_normal(ti::i)); const auto contracted_linear_momentum_integrand = tenex::evaluate( linear_momentum_integrand(ti::I, ti::J) * conformal_face_normal(ti::j)); + const auto contracted_center_of_mass_integrand = tenex::evaluate( + center_of_mass_integrand(ti::I, ti::J) * conformal_face_normal(ti::j)); // Take integrals adm_mass->get() += definite_integral( @@ -141,9 +143,10 @@ void local_adm_integrals( definite_integral(contracted_linear_momentum_integrand.get(I) * get(conformal_area_element), face_mesh); - center_of_mass->get(I) += definite_integral( - center_of_mass_integrand.get(I) * get(conformal_area_element), - face_mesh); + center_of_mass->get(I) += + definite_integral(contracted_center_of_mass_integrand.get(I) * + get(conformal_area_element), + face_mesh); } } } diff --git a/src/PointwiseFunctions/Xcts/AdmLinearMomentum.cpp b/src/PointwiseFunctions/Xcts/AdmLinearMomentum.cpp index 1d06c32f18d5..8e47fd4cd569 100644 --- a/src/PointwiseFunctions/Xcts/AdmLinearMomentum.cpp +++ b/src/PointwiseFunctions/Xcts/AdmLinearMomentum.cpp @@ -3,19 +3,6 @@ #include "PointwiseFunctions/Xcts/AdmLinearMomentum.hpp" -#include - -#include "DataStructures/DataVector.hpp" -#include "DataStructures/Tensor/Slice.hpp" -#include "DataStructures/Tensor/Tensor.hpp" -#include "Domain/Structure/Direction.hpp" -#include "Domain/Structure/IndexToSliceAt.hpp" -#include "NumericalAlgorithms/LinearOperators/Divergence.hpp" -#include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp" -#include "NumericalAlgorithms/Spectral/Mesh.hpp" -#include "Utilities/GenerateInstantiations.hpp" -#include "Utilities/Gsl.hpp" - namespace Xcts { void adm_linear_momentum_surface_integrand( @@ -47,11 +34,13 @@ void adm_linear_momentum_volume_integrand( gsl::not_null*> result, const tnsr::II& surface_integrand, const Scalar& conformal_factor, - const tnsr::i& conformal_factor_deriv, + const tnsr::i& deriv_conformal_factor, const tnsr::ii& conformal_metric, const tnsr::II& inv_conformal_metric, const tnsr::Ijj& conformal_christoffel_second_kind, const tnsr::i& conformal_christoffel_contracted) { + // Note: we can ignore the $1/(8\pi)$ term below because it is already + // included in `surface_integrand`. tenex::evaluate( result, -(conformal_christoffel_second_kind(ti::I, ti::j, ti::k) * @@ -59,14 +48,14 @@ void adm_linear_momentum_volume_integrand( conformal_christoffel_contracted(ti::k) * surface_integrand(ti::I, ti::K) - 2. * conformal_metric(ti::j, ti::k) * surface_integrand(ti::J, ti::K) * - inv_conformal_metric(ti::I, ti::L) * conformal_factor_deriv(ti::l) / + inv_conformal_metric(ti::I, ti::L) * deriv_conformal_factor(ti::l) / conformal_factor())); } tnsr::I adm_linear_momentum_volume_integrand( const tnsr::II& surface_integrand, const Scalar& conformal_factor, - const tnsr::i& conformal_factor_deriv, + const tnsr::i& deriv_conformal_factor, const tnsr::ii& conformal_metric, const tnsr::II& inv_conformal_metric, const tnsr::Ijj& conformal_christoffel_second_kind, @@ -74,7 +63,7 @@ tnsr::I adm_linear_momentum_volume_integrand( tnsr::I result; adm_linear_momentum_volume_integrand( make_not_null(&result), surface_integrand, conformal_factor, - conformal_factor_deriv, conformal_metric, inv_conformal_metric, + deriv_conformal_factor, conformal_metric, inv_conformal_metric, conformal_christoffel_second_kind, conformal_christoffel_contracted); return result; } diff --git a/src/PointwiseFunctions/Xcts/AdmLinearMomentum.hpp b/src/PointwiseFunctions/Xcts/AdmLinearMomentum.hpp index 4d25dbe8981a..7e3a32737085 100644 --- a/src/PointwiseFunctions/Xcts/AdmLinearMomentum.hpp +++ b/src/PointwiseFunctions/Xcts/AdmLinearMomentum.hpp @@ -4,30 +4,26 @@ #pragma once #include "DataStructures/DataVector.hpp" -#include "DataStructures/Tensor/Slice.hpp" #include "DataStructures/Tensor/Tensor.hpp" -#include "Domain/Structure/Direction.hpp" -#include "Domain/Structure/IndexToSliceAt.hpp" -#include "NumericalAlgorithms/Spectral/Mesh.hpp" #include "Utilities/Gsl.hpp" -#include - -#include "NumericalAlgorithms/LinearOperators/Divergence.hpp" -#include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp" - namespace Xcts { /// @{ /*! - * \brief Surface integrand for ADM linear momentum calculation defined as (see - * Eq. 20 in \cite Ossokine2015yla): + * \brief Surface integrand for the ADM linear momentum calculation. + * + * We define the ADM linear momentum integral as (see Eqs. 19-20 in + * \cite Ossokine2015yla): * * \begin{equation} - * \frac{1}{8\pi} \psi^{10} (K^{ij} - K \gamma^{ij}) + * P_\text{ADM}^i = \frac{1}{8\pi} + * \oint_{S_\infty} \psi^10 \Big( + * K^{ij} - K \gamma^{ij} + * \Big) \, dS_j. * \end{equation} * - * \param result output buffer for the surface integrand + * \param result output pointer * \param conformal_factor the conformal factor $\psi$ * \param inv_spatial_metric the inverse spatial metric $\gamma^{ij}$ * \param inv_extrinsic_curvature the inverse extrinsic curvature $K^{ij}$ @@ -54,23 +50,23 @@ tnsr::II adm_linear_momentum_surface_integrand( * Eq. 20 in \cite Ossokine2015yla): * * \begin{equation} - * - \frac{1}{8\pi} ( - * \bar\Gamma^i_{jk} P^{jk} - * + \bar\Gamma^j_{jk} P^{jk} - * - 2 \bar\gamma_{jk} P^{jk} \bar\gamma^{il} \partial_l(\ln\psi) - * ), + * P_\text{ADM}^i = - \frac{1}{8\pi} + * \int_{V_\infty} \Big( + * \bar\Gamma^i_{jk} P^{jk} + * + \bar\Gamma^j_{jk} P^{jk} + * - 2 \bar\gamma_{jk} P^{jk} \bar\gamma^{il} + * \partial_l(\ln\psi) + * \Big) \, dV, * \end{equation} * - * where $\frac{1}{8\pi} P^{jk}$ is the result from + * where $1/(8\pi) P^{jk}$ is the result from * `adm_linear_momentum_surface_integrand`. * - * Note that we are including the negative sign in the integrand. - * - * \param result output buffer for the surface integrand - * \param surface_integrand the result of - * `adm_linear_momentum_surface_integrand` + * \param result output pointer + * \param surface_integrand the quantity $1/(8\pi) P^{ij}$ (result of + * `adm_linear_momentum_surface_integrand`) * \param conformal_factor the conformal factor $\psi$ - * \param conformal_factor_deriv the derivative of the conformal factor + * \param deriv_conformal_factor the gradient of the conformal factor * $\partial_i\psi$ * \param conformal_metric the conformal metric $\bar\gamma_{ij}$ * \param inv_conformal_metric the inverse conformal metric $\bar\gamma^{ij}$ @@ -83,7 +79,7 @@ void adm_linear_momentum_volume_integrand( gsl::not_null*> result, const tnsr::II& surface_integrand, const Scalar& conformal_factor, - const tnsr::i& conformal_factor_deriv, + const tnsr::i& deriv_conformal_factor, const tnsr::ii& conformal_metric, const tnsr::II& inv_conformal_metric, const tnsr::Ijj& conformal_christoffel_second_kind, @@ -93,7 +89,7 @@ void adm_linear_momentum_volume_integrand( tnsr::I adm_linear_momentum_volume_integrand( const tnsr::II& surface_integrand, const Scalar& conformal_factor, - const tnsr::i& conformal_factor_deriv, + const tnsr::i& deriv_conformal_factor, const tnsr::ii& conformal_metric, const tnsr::II& inv_conformal_metric, const tnsr::Ijj& conformal_christoffel_second_kind, diff --git a/src/PointwiseFunctions/Xcts/AdmMass.cpp b/src/PointwiseFunctions/Xcts/AdmMass.cpp index e635f064c67c..918f571526e2 100644 --- a/src/PointwiseFunctions/Xcts/AdmMass.cpp +++ b/src/PointwiseFunctions/Xcts/AdmMass.cpp @@ -33,4 +33,56 @@ tnsr::I adm_mass_surface_integrand( return result; } +void adm_mass_volume_integrand( + gsl::not_null*> result, + const Scalar& conformal_factor, + const Scalar& conformal_ricci_scalar, + const Scalar& trace_extrinsic_curvature, + const Scalar& + longitudinal_shift_minus_dt_conformal_metric_over_lapse_square, + const tnsr::II& inv_conformal_metric, + const tnsr::iJK& deriv_inv_conformal_metric, + const tnsr::Ijj& conformal_christoffel_second_kind, + const tnsr::iJkk& deriv_conformal_christoffel_second_kind) { + tenex::evaluate( + result, + 1. / (16. * M_PI) * + (inv_conformal_metric(ti::J, ti::K) * + deriv_conformal_christoffel_second_kind(ti::i, ti::I, ti::j, + ti::k) + + conformal_christoffel_second_kind(ti::I, ti::j, ti::k) * + deriv_inv_conformal_metric(ti::i, ti::J, ti::K) - + inv_conformal_metric(ti::I, ti::J) * + deriv_conformal_christoffel_second_kind(ti::i, ti::K, ti::k, + ti::j) - + conformal_christoffel_second_kind(ti::K, ti::k, ti::j) * + deriv_inv_conformal_metric(ti::i, ti::I, ti::J) - + conformal_factor() * conformal_ricci_scalar() - + 2. / 3. * pow<5>(conformal_factor()) * + square(trace_extrinsic_curvature()) + + pow<5>(conformal_factor()) / 4. * + longitudinal_shift_minus_dt_conformal_metric_over_lapse_square())); +} + +Scalar adm_mass_volume_integrand( + const Scalar& conformal_factor, + const Scalar& conformal_ricci_scalar, + const Scalar& trace_extrinsic_curvature, + const Scalar& + longitudinal_shift_minus_dt_conformal_metric_over_lapse_square, + const tnsr::II& inv_conformal_metric, + const tnsr::iJK& deriv_inv_conformal_metric, + const tnsr::Ijj& conformal_christoffel_second_kind, + const tnsr::iJkk& deriv_conformal_christoffel_second_kind) { + Scalar result; + adm_mass_volume_integrand( + make_not_null(&result), conformal_factor, conformal_ricci_scalar, + trace_extrinsic_curvature, + longitudinal_shift_minus_dt_conformal_metric_over_lapse_square, + inv_conformal_metric, deriv_inv_conformal_metric, + conformal_christoffel_second_kind, + deriv_conformal_christoffel_second_kind); + return result; +} + } // namespace Xcts diff --git a/src/PointwiseFunctions/Xcts/AdmMass.hpp b/src/PointwiseFunctions/Xcts/AdmMass.hpp index 7285df116770..6053d318c2b1 100644 --- a/src/PointwiseFunctions/Xcts/AdmMass.hpp +++ b/src/PointwiseFunctions/Xcts/AdmMass.hpp @@ -16,20 +16,20 @@ namespace Xcts { * We define the ADM mass integral as (see Eq. 3.139 in \cite BaumgarteShapiro): * * \begin{equation} - * M_{ADM} - * = \int_{S_\infty} \frac{1}{16\pi} ( - * \bar\gamma^{jk} \bar\Gamma^i_{jk} - * - \bar\gamma^{ij} \bar\Gamma_{j} - * - 8 \bar\gamma^{ij} \partial_j \psi - * ) d\bar{S}_i. + * M_\text{ADM} = \frac{1}{16\pi} + * \oint_{S_\infty} \Big( + * \bar\gamma^{jk} \bar\Gamma^i_{jk} + * - \bar\gamma^{ij} \bar\Gamma_{j} + * - 8 \bar\gamma^{ij} \partial_j \psi + * \Big) dS_i. * \end{equation} * - * Note that we don't use the other versions presented in \cite BaumgarteShapiro - * of this integral because they make assumptions like $\bar\gamma = 1$, - * $\bar\Gamma^i_{ij} = 0$ and fast fall-off of the conformal metric. + * \note We don't use the other versions presented in \cite BaumgarteShapiro of + * this integral because they make assumptions like $\bar\gamma = 1$, + * $\bar\Gamma^i_{ij} = 0$ and faster fall-off of the conformal metric. * - * \param result output buffer for the surface integrand - * \param deriv_conformal_factor the partial derivatives of the conformal factor + * \param result output pointer + * \param deriv_conformal_factor the gradient of the conformal factor * $\partial_i \psi$ * \param inv_conformal_metric the inverse conformal metric $\bar\gamma^{ij}$ * \param conformal_christoffel_second_kind the conformal christoffel symbol @@ -52,4 +52,74 @@ tnsr::I adm_mass_surface_integrand( const tnsr::i& conformal_christoffel_contracted); /// @} +/// @{ +/*! + * \brief Volume integrand for the ADM mass calculation. + * + * We cast the ADM mass as an infinite volume integral by applying Gauss' law on + * the surface integral defined in `adm_mass_surface_integrand`: + * + * \begin{equation} + * M_\text{ADM} = \frac{1}{16\pi} + * \int_{V_\infty} \Big( + * \bar\gamma^{jk} \partial_i \bar\Gamma^i_{jk} + * + \bar\Gamma^i_{jk} \partial_i \bar\gamma^{jk} + * - \bar\gamma^{ij} \partial_i \bar\Gamma_j + * - \bar\Gamma_j \partial_i \bar\gamma^{ij} + * - 8 \bar D^2 \psi + * \Big) dV, + * \end{equation} + * + * where we can use the momentum constraint to replace $8 \bar D^2 \psi$ with + * + * \begin{equation} + * 8 \bar D^2 \psi = \psi \bar R + \frac{2}{3} \psi^5 K^2 + * - \frac{1}{4} \psi^5 \frac{1}{\alpha^2} + * \Big[ (\bar L \beta)_{ij} - \bar u_{ij} \Big] + * \Big[ (\bar L \beta)^{ij} - \bar u^{ij} \Big]. + * \end{equation} + * + * \note This is similar to Eq. 3.149 in \cite BaumgarteShapiro, except that + * here we don't assume $\bar\gamma = 1$. + * + * \param result output pointer + * \param conformal_factor the conformal factor + * \param conformal_ricci_scalar the conformal Ricci scalar $\bar R$ + * \param trace_extrinsic_curvature the extrinsic curvature trace $K$ + * \param longitudinal_shift_minus_dt_conformal_metric_over_lapse_square the + * quantity computed in + * `Xcts::Tags::LongitudinalShiftMinusDtConformalMetricOverLapseSquare` + * \param inv_conformal_metric the inverse conformal metric $\bar\gamma^{ij}$ + * \param deriv_inv_conformal_metric the gradient of the inverse conformal + * metric $\partial_i \bar\gamma^{jk}$ + * \param conformal_christoffel_second_kind the conformal christoffel symbol + * $\bar\Gamma^i_{jk}$ + * \param deriv_conformal_christoffel_second_kind the gradient of the conformal + * christoffel symbol $\partial_i \bar\Gamma^j_{kl}$ + */ +void adm_mass_volume_integrand( + gsl::not_null*> result, + const Scalar& conformal_factor, + const Scalar& conformal_ricci_scalar, + const Scalar& trace_extrinsic_curvature, + const Scalar& + longitudinal_shift_minus_dt_conformal_metric_over_lapse_square, + const tnsr::II& inv_conformal_metric, + const tnsr::iJK& deriv_inv_conformal_metric, + const tnsr::Ijj& conformal_christoffel_second_kind, + const tnsr::iJkk& deriv_conformal_christoffel_second_kind); + +/// Return-by-value overload +Scalar adm_mass_volume_integrand( + const Scalar& conformal_factor, + const Scalar& conformal_ricci_scalar, + const Scalar& trace_extrinsic_curvature, + const Scalar& + longitudinal_shift_minus_dt_conformal_metric_over_lapse_square, + const tnsr::II& inv_conformal_metric, + const tnsr::iJK& deriv_inv_conformal_metric, + const tnsr::Ijj& conformal_christoffel_second_kind, + const tnsr::iJkk& deriv_conformal_christoffel_second_kind); +/// @} + } // namespace Xcts diff --git a/src/PointwiseFunctions/Xcts/CenterOfMass.cpp b/src/PointwiseFunctions/Xcts/CenterOfMass.cpp index c157aa7a0589..1010407e4bdf 100644 --- a/src/PointwiseFunctions/Xcts/CenterOfMass.cpp +++ b/src/PointwiseFunctions/Xcts/CenterOfMass.cpp @@ -6,20 +6,42 @@ namespace Xcts { void center_of_mass_surface_integrand( - gsl::not_null*> result, + gsl::not_null*> result, const Scalar& conformal_factor, const tnsr::I& unit_normal) { - tenex::evaluate(result, 3. / (8. * M_PI) * pow<4>(conformal_factor()) * - unit_normal(ti::I)); + tenex::evaluate(result, + 3. / (8. * M_PI) * pow<4>(conformal_factor()) * + unit_normal(ti::I) * unit_normal(ti::J)); } -tnsr::I center_of_mass_surface_integrand( +tnsr::II center_of_mass_surface_integrand( const Scalar& conformal_factor, const tnsr::I& unit_normal) { - tnsr::I result; + tnsr::II result; center_of_mass_surface_integrand(make_not_null(&result), conformal_factor, unit_normal); return result; } +void center_of_mass_volume_integrand( + gsl::not_null*> result, + const Scalar& conformal_factor, + const tnsr::i& deriv_conformal_factor, + const tnsr::I& unit_normal) { + tenex::evaluate(result, 3. / (8. * M_PI) * 4. * + pow<3>(conformal_factor()) * + deriv_conformal_factor(ti::j) * + unit_normal(ti::I) * unit_normal(ti::J)); +} + +tnsr::I center_of_mass_volume_integrand( + const Scalar& conformal_factor, + const tnsr::i& deriv_conformal_factor, + const tnsr::I& unit_normal) { + tnsr::I result; + center_of_mass_volume_integrand(make_not_null(&result), conformal_factor, + deriv_conformal_factor, unit_normal); + return result; +} + } // namespace Xcts diff --git a/src/PointwiseFunctions/Xcts/CenterOfMass.hpp b/src/PointwiseFunctions/Xcts/CenterOfMass.hpp index 24c49f1a9608..a8250862eb3b 100644 --- a/src/PointwiseFunctions/Xcts/CenterOfMass.hpp +++ b/src/PointwiseFunctions/Xcts/CenterOfMass.hpp @@ -17,14 +17,14 @@ namespace Xcts { * \cite Ossokine2015yla): * * \begin{equation} - * C_{CoM}^i - * = \frac{1}{M_{ADM}} \int_{S_\infty} \frac{3}{8\pi} \psi^4 n^i \, d\bar{A}. + * C_\text{CoM}^i = \frac{3}{8 \pi M_\text{ADM}} + * \oint_{S_\infty} \psi^4 n^i n^j \, dS_j * \end{equation} * - * Note that we don't include the ADM mass $M_{ADM}$ in this integrand. After + * \note We don't include the ADM mass $M_{ADM}$ in this integrand. After * integrating the result of this function, you have to divide by $M_{ADM}$. - * See `Xcts::adm_mass_surface_integrand` for details on how to calculate - * $M_{ADM}$. + * + * \see `Xcts::adm_mass_surface_integrand` * * \warning This integral assumes that the conformal metric falls off to * flatness faster than $1/r^2$. That means that it cannot be directly used @@ -32,18 +32,51 @@ namespace Xcts { * for XCTS with Superposed Kerr-Schild (SKS) because of the exponential * fall-off terms. * - * \param result output buffer for the surface integrand + * \param result output pointer * \param conformal_factor the conformal factor $\psi$ * \param unit_normal the outward-pointing unit normal $n^i = x^i / r$ */ void center_of_mass_surface_integrand( + gsl::not_null*> result, + const Scalar& conformal_factor, + const tnsr::I& unit_normal); + +/// Return-by-value overload +tnsr::II center_of_mass_surface_integrand( + const Scalar& conformal_factor, + const tnsr::I& unit_normal); +/// @} + +/// @{ +/*! + * \brief Volume integrand for the center of mass calculation. + * + * We cast the center of mass as an infinite volume integral by applying Gauss' + * law on the surface integral defined in `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. + * \end{equation} + * + * \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$ + */ +void center_of_mass_volume_integrand( gsl::not_null*> result, const Scalar& conformal_factor, + const tnsr::i& deriv_conformal_factor, const tnsr::I& unit_normal); /// Return-by-value overload -tnsr::I center_of_mass_surface_integrand( +tnsr::I center_of_mass_volume_integrand( const Scalar& conformal_factor, + const tnsr::i& deriv_conformal_factor, const tnsr::I& unit_normal); /// @} diff --git a/tests/Unit/Elliptic/Systems/Xcts/Events/Test_ObserveAdmIntegrals.cpp b/tests/Unit/Elliptic/Systems/Xcts/Events/Test_ObserveAdmIntegrals.cpp index 0393adae01f0..6551333a7532 100644 --- a/tests/Unit/Elliptic/Systems/Xcts/Events/Test_ObserveAdmIntegrals.cpp +++ b/tests/Unit/Elliptic/Systems/Xcts/Events/Test_ObserveAdmIntegrals.cpp @@ -25,11 +25,13 @@ #include "PointwiseFunctions/GeneralRelativity/Christoffel.hpp" #include "PointwiseFunctions/GeneralRelativity/ExtrinsicCurvature.hpp" #include "PointwiseFunctions/GeneralRelativity/Lapse.hpp" +#include "PointwiseFunctions/GeneralRelativity/Ricci.hpp" #include "PointwiseFunctions/GeneralRelativity/Shift.hpp" #include "PointwiseFunctions/GeneralRelativity/SpacetimeMetric.hpp" #include "PointwiseFunctions/GeneralRelativity/SpatialMetric.hpp" #include "PointwiseFunctions/SpecialRelativity/LorentzBoostMatrix.hpp" #include "PointwiseFunctions/Xcts/ExtrinsicCurvature.hpp" +#include "PointwiseFunctions/Xcts/LongitudinalOperator.hpp" namespace { @@ -180,8 +182,7 @@ void test_local_adm_integrals(const double& distance, const auto inv_conformal_metric = tenex::evaluate( inv_spatial_metric(ti::I, ti::J) * pow<4>(conformal_factor())); - // Compute spatial derivatives (needed for extrinsic curvature). - const auto deriv_lapse = partial_derivative(lapse, mesh, inv_jacobian); + // Compute spatial derivatives. const auto deriv_shift = partial_derivative(shift, mesh, inv_jacobian); const auto deriv_spatial_metric = partial_derivative(spatial_metric, mesh, inv_jacobian); @@ -195,6 +196,15 @@ void test_local_adm_integrals(const double& distance, deriv_conformal_metric, inv_conformal_metric); const auto conformal_christoffel_contracted = tenex::evaluate( conformal_christoffel_second_kind(ti::J, ti::i, ti::j)); + const auto deriv_conformal_christoffel_second_kind = partial_derivative( + conformal_christoffel_second_kind, mesh, inv_jacobian); + + // Compute conformal Ricci scalar. + const auto conformal_ricci_tensor = + gr::ricci_tensor(conformal_christoffel_second_kind, + deriv_conformal_christoffel_second_kind); + const auto conformal_ricci_scalar = + gr::ricci_scalar(conformal_ricci_tensor, inv_conformal_metric); // Define variables that appear in the formulas of dt_spatial_metric. const auto& x = get<0>(inertial_coords); @@ -203,7 +213,7 @@ void test_local_adm_integrals(const double& distance, const auto barred_r = sqrt(square(x) + square(y) + square(lorentz_factor * z)); - // Compute spatial metric time derivative (needed for extrinsic curvature). + // Compute spatial metric time derivative. // Note that these formulas were derived in a Mathematica notebook // specifically for this problem. Here, we are evaluating them at t = 0. auto dt_spatial_metric = @@ -221,6 +231,17 @@ void test_local_adm_integrals(const double& distance, (2 * mass * cube(boost_speed) * z * cube(1 + mass / barred_r)) / (square(1 - square(boost_speed)) * cube(barred_r)); + // Compute longitudinal operators. + tnsr::II longitudinal_shift; + Xcts::longitudinal_operator(make_not_null(&longitudinal_shift), shift, + deriv_shift, inv_conformal_metric, + conformal_christoffel_second_kind); + const auto longitudinal_shift_background_minus_dt_conformal_metric = + tenex::evaluate(longitudinal_shift(ti::I, ti::J)); + const auto longitudinal_shift_excess = + make_with_value>( + inertial_coords, 0.0); + // Compute extrinsic curvature and its trace. const auto extrinsic_curvature = gr::extrinsic_curvature(lapse, shift, deriv_shift, spatial_metric, diff --git a/tests/Unit/PointwiseFunctions/Xcts/Test_AdmLinearMomentum.cpp b/tests/Unit/PointwiseFunctions/Xcts/Test_AdmLinearMomentum.cpp index f8835f237812..f2a15b11aa60 100644 --- a/tests/Unit/PointwiseFunctions/Xcts/Test_AdmLinearMomentum.cpp +++ b/tests/Unit/PointwiseFunctions/Xcts/Test_AdmLinearMomentum.cpp @@ -6,6 +6,7 @@ #include "DataStructures/DataVector.hpp" #include "DataStructures/Tensor/EagerMath/Determinant.hpp" #include "DataStructures/Tensor/EagerMath/Magnitude.hpp" +#include "DataStructures/Tensor/Slice.hpp" #include "DataStructures/Tensor/Tensor.hpp" #include "Domain/AreaElement.hpp" #include "Domain/CreateInitialElement.hpp" @@ -31,15 +32,14 @@ using Schwarzschild = Xcts::Solutions::Schwarzschild; using KerrSchild = Xcts::Solutions::WrappedGr; template -void test_linear_momentum_surface_integral(const double distance, - const double mass, - const double boost_speed, - const Solution& solution, - const double horizon_radius) { +void test_infinite_surface_integral(const double distance, const double mass, + const double horizon_radius, + const double boost_speed, + const Solution& solution) { // Set up domain const size_t h_refinement = 1; const size_t p_refinement = 6; - domain::creators::Sphere shell{ + const domain::creators::Sphere shell{ /* inner_radius */ horizon_radius, /* outer_radius */ distance, /* interior */ domain::creators::Sphere::Excision{}, @@ -167,10 +167,220 @@ void test_linear_momentum_surface_integral(const double distance, custom_approx(lorentz_factor * mass * boost_speed)); } +template +void test_infinite_volume_integral(const double distance, const double mass, + const double horizon_radius, + const double boost_speed, + const Solution& solution) { + // Set up domain + const size_t h_refinement = 1; + const size_t p_refinement = 6; + const domain::creators::Sphere shell{ + /* inner_radius */ 2 * horizon_radius, + /* outer_radius */ distance, + /* interior */ domain::creators::Sphere::Excision{}, + /* initial_refinement */ h_refinement, + /* initial_number_of_grid_points */ p_refinement + 1, + /* use_equiangular_map */ true, + /* equatorial_compression */ {}, + /* radial_partitioning */ {}, + /* radial_distribution */ domain::CoordinateMaps::Distribution::Inverse}; + const auto shell_domain = shell.create_domain(); + const auto& blocks = shell_domain.blocks(); + const auto& initial_ref_levels = shell.initial_refinement_levels(); + const auto element_ids = initial_element_ids(initial_ref_levels); + const Mesh<3> mesh{p_refinement + 1, Spectral::Basis::Legendre, + Spectral::Quadrature::GaussLobatto}; + const Mesh<2> face_mesh{p_refinement + 1, Spectral::Basis::Legendre, + Spectral::Quadrature::GaussLobatto}; + + // Initialize surface integral + tnsr::I total_integral({0., 0., 0.}); + + // Compute integrals by summing over each element + for (const auto& element_id : element_ids) { + // Get element information + const auto& current_block = blocks.at(element_id.block_id()); + const auto current_element = domain::Initialization::create_initial_element( + element_id, current_block, initial_ref_levels); + const ElementMap<3, Frame::Inertial> logical_to_inertial_map( + element_id, current_block.stationary_map().get_clone()); + + // Get 3D coordinates + const auto logical_coords = logical_coordinates(mesh); + const auto inertial_coords = logical_to_inertial_map(logical_coords); + const auto jacobian = logical_to_inertial_map.jacobian(logical_coords); + const auto det_jacobian = determinant(jacobian); + const auto inv_jacobian = + logical_to_inertial_map.inv_jacobian(logical_coords); + + // Get required fields + const auto solution_fields = solution.variables( + inertial_coords, + tmpl::list< + Xcts::Tags::ConformalFactor, + ::Tags::deriv, + tmpl::size_t<3>, Frame::Inertial>, + Xcts::Tags::ConformalMetric, + Xcts::Tags::InverseConformalMetric, + gr::Tags::InverseSpatialMetric, + gr::Tags::ExtrinsicCurvature, + gr::Tags::TraceExtrinsicCurvature, + Xcts::Tags::ConformalChristoffelSecondKind, + Xcts::Tags::ConformalChristoffelContracted>{}); + const auto& conformal_factor = + get>(solution_fields); + const auto& deriv_conformal_factor = + get<::Tags::deriv, + tmpl::size_t<3>, Frame::Inertial>>(solution_fields); + const auto& conformal_metric = + get>( + solution_fields); + const auto& inv_conformal_metric = + get>( + solution_fields); + const auto& inv_spatial_metric = + get>( + solution_fields); + const auto& extrinsic_curvature = + get>( + solution_fields); + const auto& trace_extrinsic_curvature = + get>(solution_fields); + const auto& conformal_christoffel_second_kind = + get>( + solution_fields); + const auto& conformal_christoffel_contracted = + get>( + solution_fields); + + // Compute the inverse extrinsic curvature + tnsr::II inv_extrinsic_curvature; + tenex::evaluate(make_not_null(&inv_extrinsic_curvature), + inv_spatial_metric(ti::I, ti::K) * + inv_spatial_metric(ti::J, ti::L) * + extrinsic_curvature(ti::k, ti::l)); + + const auto surface_integrand = Xcts::adm_linear_momentum_surface_integrand( + conformal_factor, inv_spatial_metric, inv_extrinsic_curvature, + trace_extrinsic_curvature); + + const auto volume_integrand = Xcts::adm_linear_momentum_volume_integrand( + surface_integrand, conformal_factor, deriv_conformal_factor, + conformal_metric, inv_conformal_metric, + conformal_christoffel_second_kind, conformal_christoffel_contracted); + for (int I = 0; I < 3; I++) { + total_integral.get(I) += + definite_integral(volume_integrand.get(I) * get(det_jacobian), mesh); + } + + // Loop over external boundaries + for (auto boundary_direction : current_element.external_boundaries()) { + // Skip interfaces not at the outer boundary + if (boundary_direction != Direction<3>::lower_zeta()) { + continue; + } + + // Get interface coordinates. + const auto face_logical_coords = + interface_logical_coordinates(face_mesh, boundary_direction); + const auto face_inv_jacobian = + logical_to_inertial_map.inv_jacobian(face_logical_coords); + + // Slice required fields to the interface + const size_t slice_index = + index_to_slice_at(mesh.extents(), boundary_direction); + const auto& face_surface_integrand = + data_on_slice(surface_integrand, mesh.extents(), + boundary_direction.dimension(), slice_index); + + // Compute Euclidean area element + const auto area_element = + euclidean_area_element(face_inv_jacobian, boundary_direction); + + // Compute Euclidean face normal + auto euclidean_face_normal = unnormalized_face_normal( + face_mesh, logical_to_inertial_map, boundary_direction); + const auto face_normal_magnitude = magnitude(euclidean_face_normal); + for (size_t d = 0; d < 3; ++d) { + euclidean_face_normal.get(d) /= get(face_normal_magnitude); + } + + // Compute surface integrand + const auto contracted_integrand = tenex::evaluate( + -face_surface_integrand(ti::I, ti::J) * euclidean_face_normal(ti::j)); + + // Compute contribution to surface integral + for (int I = 0; I < 3; I++) { + total_integral.get(I) += definite_integral( + contracted_integrand.get(I) * get(area_element), face_mesh); + } + } + } + + // Check result + const double lorentz_factor = 1. / sqrt(1. - square(boost_speed)); + auto custom_approx = Approx::custom().epsilon(10. / distance).scale(1.0); + CHECK(get<0>(total_integral) == custom_approx(0.)); + CHECK(get<1>(total_integral) == custom_approx(0.)); + CHECK(get<2>(total_integral) == + custom_approx(lorentz_factor * mass * boost_speed)); +} + } // namespace SPECTRE_TEST_CASE("Unit.PointwiseFunctions.Xcts.AdmLinearMomentum", "[Unit][PointwiseFunctions]") { + { + INFO("Schwarzschild in Kerr-Schild coordinates"); + const double mass = 1.; + const double horizon_radius = 2. * mass; + const double boost_speed = 0.; + const Xcts::Solutions::WrappedGr solution( + mass, std::array{{0., 0., 0.}}, + std::array{{0., 0., 0.}}, + std::array{{0., 0., boost_speed}}); + for (const double distance : std::array({1.e4, 1.e5, 1.e6})) { + test_infinite_surface_integral(distance, mass, horizon_radius, + boost_speed, solution); + test_infinite_volume_integral(distance, mass, horizon_radius, boost_speed, + solution); + } + } + { + INFO("Boosted Schwarzschild in Kerr-Schild coordinates"); + const double mass = 1.; + const double horizon_radius = 2. * mass; + const double boost_speed = 0.5; + const Xcts::Solutions::WrappedGr solution( + mass, std::array{{0., 0., 0.}}, + std::array{{0., 0., 0.}}, + std::array{{0., 0., boost_speed}}); + for (const double distance : std::array({1.e4, 1.e5, 1.e6})) { + test_infinite_surface_integral(distance, mass, horizon_radius, + boost_speed, solution); + // Note: the volume integral currently doesn't work with this test case. + } + } + { + INFO("Schwarzschild in isotropic coordinates"); + const double mass = 1.; + const double horizon_radius = 0.5 * mass; + const double boost_speed = 0.; + const Xcts::Solutions::Schwarzschild solution( + mass, Xcts::Solutions::SchwarzschildCoordinates::Isotropic); + for (const double distance : std::array({1.e4, 1.e5, 1.e6})) { + test_infinite_surface_integral(distance, mass, horizon_radius, + boost_speed, solution); + test_infinite_volume_integral(distance, mass, horizon_radius, boost_speed, + solution); + } + } + // Test integrands against Python implementation with random values. const pypp::SetupLocalPythonEnvironment local_python_env{ "PointwiseFunctions/Xcts"}; @@ -193,32 +403,4 @@ SPECTRE_TEST_CASE("Unit.PointwiseFunctions.Xcts.AdmLinearMomentum", &Xcts::adm_linear_momentum_volume_integrand), "AdmLinearMomentum", {"adm_linear_momentum_volume_integrand"}, {{{-1, 1.}}}, used_for_size); - - // Test that integral converges with two analytic solutions. - { - INFO("Boosted Kerr-Schild"); - const double mass = 1.; - const double horizon_radius = 2. * mass; - const double boost_speed = 0.5; - const std::array boost_velocity({0., 0., boost_speed}); - const std::array dimensionless_spin({0., 0., 0.}); - const std::array center({0., 0., 0.}); - const KerrSchild solution(mass, dimensionless_spin, center, boost_velocity); - for (const double distance : std::array({1.e3, 1.e5, 1.e10})) { - test_linear_momentum_surface_integral(distance, mass, boost_speed, - solution, horizon_radius); - } - } - { - INFO("Isotropic Schwarzschild"); - const double mass = 1.; - const double horizon_radius = 0.5 * mass; - const double boost_speed = 0.; - const Schwarzschild solution( - mass, Xcts::Solutions::SchwarzschildCoordinates::Isotropic); - for (const double distance : std::array({1.e3, 1.e5, 1.e10})) { - test_linear_momentum_surface_integral(distance, mass, boost_speed, - solution, horizon_radius); - } - } } diff --git a/tests/Unit/PointwiseFunctions/Xcts/Test_AdmMass.cpp b/tests/Unit/PointwiseFunctions/Xcts/Test_AdmMass.cpp index d10dfc6a7708..a09c7d34c718 100644 --- a/tests/Unit/PointwiseFunctions/Xcts/Test_AdmMass.cpp +++ b/tests/Unit/PointwiseFunctions/Xcts/Test_AdmMass.cpp @@ -6,6 +6,7 @@ #include "DataStructures/DataVector.hpp" #include "DataStructures/Tensor/EagerMath/Determinant.hpp" #include "DataStructures/Tensor/EagerMath/Magnitude.hpp" +#include "DataStructures/Tensor/Slice.hpp" #include "DataStructures/Tensor/Tensor.hpp" #include "Domain/AreaElement.hpp" #include "Domain/CoordinateMaps/CoordinateMap.hpp" @@ -27,19 +28,16 @@ namespace { -using Schwarzschild = Xcts::Solutions::Schwarzschild; -using KerrSchild = Xcts::Solutions::WrappedGr; - template -void test_mass_surface_integral(const double distance, const double mass, - const double boost_speed, - const Solution& solution, - const double horizon_radius) { +void test_infinite_surface_integral(const double distance, const double mass, + const double horizon_radius, + const double boost_speed, + const Solution& solution) { // Set up domain const size_t h_refinement = 1; const size_t p_refinement = 6; - domain::creators::Sphere shell{ - /* inner_radius */ horizon_radius, + const domain::creators::Sphere shell{ + /* inner_radius */ 1.1 * horizon_radius, /* outer_radius */ distance, /* interior */ domain::creators::Sphere::Excision{}, /* initial_refinement */ h_refinement, @@ -156,35 +154,242 @@ void test_mass_surface_integral(const double distance, const double mass, CHECK(get(surface_integral) == custom_approx(lorentz_factor * mass)); } +template +void test_infinite_volume_integral(const double distance, const double mass, + const double horizon_radius, + const double boost_speed, + const Solution& solution) { + // Set up domain + const size_t h_refinement = 1; + const size_t p_refinement = 6; + const domain::creators::Sphere shell{ + /* inner_radius */ 1.1 * horizon_radius, + /* outer_radius */ distance, + /* interior */ domain::creators::Sphere::Excision{}, + /* initial_refinement */ h_refinement, + /* initial_number_of_grid_points */ p_refinement + 1, + /* use_equiangular_map */ true, + /* equatorial_compression */ {}, + /* radial_partitioning */ {}, + /* radial_distribution */ domain::CoordinateMaps::Distribution::Inverse}; + const auto shell_domain = shell.create_domain(); + const auto& blocks = shell_domain.blocks(); + const auto& initial_ref_levels = shell.initial_refinement_levels(); + const auto element_ids = initial_element_ids(initial_ref_levels); + const Mesh<3> mesh{p_refinement + 1, Spectral::Basis::Legendre, + Spectral::Quadrature::GaussLobatto}; + const Mesh<2> face_mesh{p_refinement + 1, Spectral::Basis::Legendre, + Spectral::Quadrature::GaussLobatto}; + + // Initialize "reduced" integral. + Scalar total_integral(0.); + + // Compute integrals by summing over each element + for (const auto& element_id : element_ids) { + // Get element information + const auto& current_block = blocks.at(element_id.block_id()); + const auto current_element = domain::Initialization::create_initial_element( + element_id, current_block, initial_ref_levels); + const ElementMap<3, Frame::Inertial> logical_to_inertial_map( + element_id, current_block.stationary_map().get_clone()); + + // Get 3D coordinates + const auto logical_coords = logical_coordinates(mesh); + const auto inertial_coords = logical_to_inertial_map(logical_coords); + const auto jacobian = logical_to_inertial_map.jacobian(logical_coords); + const auto det_jacobian = determinant(jacobian); + const auto inv_jacobian = + logical_to_inertial_map.inv_jacobian(logical_coords); + + // Get required fields + const auto background_fields = solution.variables( + inertial_coords, mesh, inv_jacobian, + tmpl::list< + Xcts::Tags::ConformalFactor, + ::Tags::deriv, + tmpl::size_t<3>, Frame::Inertial>, + Xcts::Tags::ConformalRicciScalar, + gr::Tags::TraceExtrinsicCurvature, + Xcts::Tags::LongitudinalShiftMinusDtConformalMetricOverLapseSquare< + DataVector>, + Xcts::Tags::ConformalMetric, + ::Tags::deriv< + Xcts::Tags::ConformalMetric, + tmpl::size_t<3>, Frame::Inertial>, + Xcts::Tags::InverseConformalMetric, + Xcts::Tags::ConformalChristoffelSecondKind, + ::Tags::deriv, + tmpl::size_t<3>, Frame::Inertial>, + Xcts::Tags::ConformalChristoffelContracted>{}); + const auto& conformal_factor = + get>(background_fields); + const auto& deriv_conformal_factor = + get<::Tags::deriv, + tmpl::size_t<3>, Frame::Inertial>>(background_fields); + const auto& conformal_ricci_scalar = + get>(background_fields); + const auto& trace_extrinsic_curvature = + get>(background_fields); + const auto& longitudinal_shift_minus_dt_conformal_metric_over_lapse_square = + get>(background_fields); + const auto& conformal_metric = + get>( + background_fields); + const auto& deriv_conformal_metric = get<::Tags::deriv< + Xcts::Tags::ConformalMetric, + tmpl::size_t<3>, Frame::Inertial>>(background_fields); + const auto& inv_conformal_metric = + get>( + background_fields); + const auto& conformal_christoffel_second_kind = + get>( + background_fields); + const auto& deriv_conformal_christoffel_second_kind = + get<::Tags::deriv, + tmpl::size_t<3>, Frame::Inertial>>(background_fields); + const auto& conformal_christoffel_contracted = + get>( + background_fields); + + const auto deriv_inv_conformal_metric = + tenex::evaluate( + inv_conformal_metric(ti::J, ti::L) * + inv_conformal_metric(ti::K, ti::M) * + (deriv_conformal_metric(ti::i, ti::l, ti::m) - + conformal_christoffel_second_kind(ti::N, ti::i, ti::l) * + conformal_metric(ti::n, ti::m) - + conformal_christoffel_second_kind(ti::N, ti::i, ti::m) * + conformal_metric(ti::l, ti::n)) - + conformal_christoffel_second_kind(ti::J, ti::i, ti::l) * + inv_conformal_metric(ti::L, ti::K) - + conformal_christoffel_second_kind(ti::K, ti::i, ti::l) * + inv_conformal_metric(ti::J, ti::L)); + + // Evaluate volume integral. + const auto volume_integrand = Xcts::adm_mass_volume_integrand( + conformal_factor, conformal_ricci_scalar, trace_extrinsic_curvature, + longitudinal_shift_minus_dt_conformal_metric_over_lapse_square, + inv_conformal_metric, deriv_inv_conformal_metric, + conformal_christoffel_second_kind, + deriv_conformal_christoffel_second_kind); + total_integral.get() += + definite_integral(get(volume_integrand) * get(det_jacobian), mesh); + + // Loop over external boundaries. + for (auto boundary_direction : current_element.external_boundaries()) { + // Skip interfaces not at the inner boundary. + if (boundary_direction != Direction<3>::lower_zeta()) { + continue; + } + + // Get interface coordinates. + const auto face_logical_coords = + interface_logical_coordinates(face_mesh, boundary_direction); + const auto face_inv_jacobian = + logical_to_inertial_map.inv_jacobian(face_logical_coords); + + // Slice required fields to the interface + const size_t slice_index = + index_to_slice_at(mesh.extents(), boundary_direction); + const auto& face_deriv_conformal_factor = + data_on_slice(deriv_conformal_factor, mesh.extents(), + boundary_direction.dimension(), slice_index); + const auto& face_inv_conformal_metric = + data_on_slice(inv_conformal_metric, mesh.extents(), + boundary_direction.dimension(), slice_index); + const auto& face_conformal_christoffel_second_kind = + data_on_slice(conformal_christoffel_second_kind, mesh.extents(), + boundary_direction.dimension(), slice_index); + const auto& face_conformal_christoffel_contracted = + data_on_slice(conformal_christoffel_contracted, mesh.extents(), + boundary_direction.dimension(), slice_index); + + // Compute Euclidean area element + const auto area_element = + euclidean_area_element(face_inv_jacobian, boundary_direction); + + // Compute Euclidean face normal + auto euclidean_face_normal = unnormalized_face_normal( + face_mesh, logical_to_inertial_map, boundary_direction); + const auto face_normal_magnitude = magnitude(euclidean_face_normal); + for (size_t d = 0; d < 3; ++d) { + euclidean_face_normal.get(d) /= get(face_normal_magnitude); + } + + // Evaluate surface integral. + const auto surface_integrand = Xcts::adm_mass_surface_integrand( + face_deriv_conformal_factor, face_inv_conformal_metric, + face_conformal_christoffel_second_kind, + face_conformal_christoffel_contracted); + const auto contracted_integrand = tenex::evaluate( + -surface_integrand(ti::I) * euclidean_face_normal(ti::i)); + + // Compute contribution to surface integral + total_integral.get() += definite_integral( + get(contracted_integrand) * get(area_element), face_mesh); + } + } + + // Check result + const double lorentz_factor = 1. / sqrt(1. - square(boost_speed)); + auto custom_approx = Approx::custom().epsilon(10. / distance).scale(1.0); + CHECK(get(total_integral) == custom_approx(lorentz_factor * mass)); +} + } // namespace SPECTRE_TEST_CASE("Unit.PointwiseFunctions.Xcts.AdmMass", "[Unit][PointwiseFunctions]") { - // Test that integral converges with two analytic solutions. { - INFO("Boosted Kerr-Schild"); + INFO("Schwarzschild in Kerr-Schild coordinates"); + const double mass = 1.; + const double horizon_radius = 2. * mass; + const double boost_speed = 0.; + const Xcts::Solutions::WrappedGr solution( + mass, std::array{{0., 0., 0.}}, + std::array{{0., 0., 0.}}, + std::array{{0., 0., boost_speed}}); + for (const double distance : std::array({1.e4, 1.e5, 1.e6})) { + test_infinite_surface_integral(distance, mass, horizon_radius, + boost_speed, solution); + test_infinite_volume_integral(distance, mass, horizon_radius, boost_speed, + solution); + } + } + { + INFO("Boosted Schwarzschild in Kerr-Schild coordinates"); const double mass = 1.; const double horizon_radius = 2. * mass; const double boost_speed = 0.5; - const std::array boost_velocity({0., 0., boost_speed}); - const std::array dimensionless_spin({0., 0., 0.}); - const std::array center({0., 0., 0.}); - const KerrSchild solution(mass, dimensionless_spin, center, boost_velocity); - for (const double distance : std::array({1.e3, 1.e5, 1.e10})) { - test_mass_surface_integral(distance, mass, boost_speed, solution, - horizon_radius); + const Xcts::Solutions::WrappedGr solution( + mass, std::array{{0., 0., 0.}}, + std::array{{0., 0., 0.}}, + std::array{{0., 0., boost_speed}}); + for (const double distance : std::array({1.e4, 1.e5, 1.e6})) { + test_infinite_surface_integral(distance, mass, horizon_radius, + boost_speed, solution); + // Note: the volume integral currently doesn't work with this test case. } } { - INFO("Isotropic Schwarzschild"); + INFO("Schwarzschild in isotropic coordinates"); const double mass = 1.; const double horizon_radius = 0.5 * mass; const double boost_speed = 0.; - const Schwarzschild solution( + const Xcts::Solutions::Schwarzschild solution( mass, Xcts::Solutions::SchwarzschildCoordinates::Isotropic); - for (const double distance : std::array({1.e3, 1.e5, 1.e10})) { - test_mass_surface_integral(distance, mass, boost_speed, solution, - horizon_radius); + for (const double distance : std::array({1.e4, 1.e5, 1.e6})) { + test_infinite_surface_integral(distance, mass, horizon_radius, + boost_speed, solution); + test_infinite_volume_integral(distance, mass, horizon_radius, boost_speed, + solution); } } } diff --git a/tests/Unit/PointwiseFunctions/Xcts/Test_CenterOfMass.cpp b/tests/Unit/PointwiseFunctions/Xcts/Test_CenterOfMass.cpp index cd7edfe9bf75..3476ccc363c3 100644 --- a/tests/Unit/PointwiseFunctions/Xcts/Test_CenterOfMass.cpp +++ b/tests/Unit/PointwiseFunctions/Xcts/Test_CenterOfMass.cpp @@ -6,6 +6,7 @@ #include "DataStructures/DataVector.hpp" #include "DataStructures/Tensor/EagerMath/Determinant.hpp" #include "DataStructures/Tensor/EagerMath/Magnitude.hpp" +#include "DataStructures/Tensor/Slice.hpp" #include "DataStructures/Tensor/Tensor.hpp" #include "Domain/AreaElement.hpp" #include "Domain/CoordinateMaps/CoordinateMap.hpp" @@ -24,25 +25,31 @@ namespace { -/** - * This test shifts the isotropic Schwarzschild solution and checks that the - * center of mass corresponds to the coordinate shift. - */ -void test_center_of_mass_surface_integral(const double distance) { +/* + The tests below shift the isotropic Schwarzschild solution and check that the + center of mass corresponds to the coordinate shift. + + We have found that the center of mass integral often diverges from the + expected result as we increase the outer radius. This could be due to + round-off errors in the calculation, making the result less accurate. Here, + we're simply checking that this deviation is below some fixed tolerance. +*/ + +constexpr double TOLERANCE = 1.e-4; + +void test_infinite_surface_integral(const double distance, + const double z_shift) { // Get Schwarzschild solution in isotropic coordinates. const double mass = 1; const Xcts::Solutions::Schwarzschild solution( mass, Xcts::Solutions::SchwarzschildCoordinates::Isotropic); const double horizon_radius = 0.5 * mass; - // Define z-shift applied to the coordinates. - const double z_shift = 2. * mass; - // Set up domain const size_t h_refinement = 1; const size_t p_refinement = 6; const domain::creators::Sphere shell{ - /* inner_radius */ z_shift + horizon_radius, + /* inner_radius */ z_shift + 1.1 * horizon_radius, /* outer_radius */ distance, /* interior */ domain::creators::Sphere::Excision{}, /* initial_refinement */ h_refinement, @@ -125,29 +132,176 @@ void test_center_of_mass_surface_integral(const double distance) { area_element(inv_jacobian, boundary_direction, inv_conformal_metric, sqrt_det_conformal_metric); - // Integrate + // Compute face normal + auto conformal_face_normal = unnormalized_face_normal( + face_mesh, logical_to_inertial_map, boundary_direction); + const auto face_normal_magnitude = + magnitude(conformal_face_normal, inv_conformal_metric); + for (size_t d = 0; d < 3; ++d) { + conformal_face_normal.get(d) /= get(face_normal_magnitude); + } + + // Evaluate surface integral const auto surface_integrand = Xcts::center_of_mass_surface_integrand( conformal_factor, conformal_unit_normal); + const auto contracted_integrand = tenex::evaluate( + surface_integrand(ti::I, ti::J) * conformal_face_normal(ti::j)); for (int I = 0; I < 3; I++) { surface_integral.get(I) += definite_integral( - surface_integrand.get(I) * get(conformal_area_element), face_mesh); + contracted_integrand.get(I) * get(conformal_area_element), + face_mesh); } } } // Check result - auto custom_approx = Approx::custom().epsilon(10. / distance).scale(1.0); + auto custom_approx = Approx::custom().epsilon(TOLERANCE).scale(1.0); CHECK(get<0>(surface_integral) == custom_approx(0.)); CHECK(get<1>(surface_integral) == custom_approx(0.)); CHECK(get<2>(surface_integral) / mass == custom_approx(z_shift)); } +void test_infinite_volume_integral(const double distance, + const double z_shift) { + // Get Schwarzschild solution in isotropic coordinates. + const double mass = 1; + const Xcts::Solutions::Schwarzschild solution( + mass, Xcts::Solutions::SchwarzschildCoordinates::Isotropic); + const double horizon_radius = 0.5 * mass; + + // Set up domain + const size_t h_refinement = 1; + const size_t p_refinement = 6; + const domain::creators::Sphere shell{ + /* inner_radius */ z_shift + 1.1 * horizon_radius, + /* outer_radius */ distance, + /* interior */ domain::creators::Sphere::Excision{}, + /* initial_refinement */ h_refinement, + /* initial_number_of_grid_points */ p_refinement + 1, + /* use_equiangular_map */ true, + /* equatorial_compression */ {}, + /* radial_partitioning */ {}, + /* radial_distribution */ domain::CoordinateMaps::Distribution::Inverse}; + const auto shell_domain = shell.create_domain(); + const auto& blocks = shell_domain.blocks(); + const auto& initial_ref_levels = shell.initial_refinement_levels(); + const auto element_ids = initial_element_ids(initial_ref_levels); + const Mesh<3> mesh{p_refinement + 1, Spectral::Basis::Legendre, + Spectral::Quadrature::GaussLobatto}; + const Mesh<2> face_mesh{p_refinement + 1, Spectral::Basis::Legendre, + Spectral::Quadrature::GaussLobatto}; + + // Initialize "reduced" integral + tnsr::I total_integral({0., 0., 0.}); + + // Compute integrals by summing over each element + for (const auto& element_id : element_ids) { + // Get element information + const auto& current_block = blocks.at(element_id.block_id()); + const auto current_element = domain::Initialization::create_initial_element( + element_id, current_block, initial_ref_levels); + const ElementMap<3, Frame::Inertial> logical_to_inertial_map( + element_id, current_block.stationary_map().get_clone()); + + // Get 3D coordinates + const auto logical_coords = logical_coordinates(mesh); + const auto inertial_coords = logical_to_inertial_map(logical_coords); + const auto jacobian = logical_to_inertial_map.jacobian(logical_coords); + const auto det_jacobian = determinant(jacobian); + const auto inv_jacobian = + logical_to_inertial_map.inv_jacobian(logical_coords); + + // Shift coordinates used to get analytic solution + auto shifted_coords = inertial_coords; + shifted_coords.get(2) -= z_shift; + + // Get required fields + const auto shifted_fields = solution.variables( + shifted_coords, + tmpl::list< + Xcts::Tags::ConformalFactor, + ::Tags::deriv, + tmpl::size_t<3>, Frame::Inertial>>{}); + const auto& conformal_factor = + get>(shifted_fields); + const auto& deriv_conformal_factor = + get<::Tags::deriv, + tmpl::size_t<3>, Frame::Inertial>>(shifted_fields); + + // Compute conformal unit normal vector + const auto r_coordinate = magnitude(inertial_coords); + const auto unit_normal = + tenex::evaluate(inertial_coords(ti::I) / r_coordinate()); + + // Evaluate volume integral. + const auto volume_integrand = Xcts::center_of_mass_volume_integrand( + conformal_factor, deriv_conformal_factor, unit_normal); + for (int I = 0; I < 3; I++) { + total_integral.get(I) += + definite_integral(volume_integrand.get(I) * get(det_jacobian), mesh); + } + + // Loop over external boundaries + for (auto boundary_direction : current_element.external_boundaries()) { + // Skip interfaces not at the inner boundary + if (boundary_direction != Direction<3>::lower_zeta()) { + continue; + } + + // Get interface coordinates. + const auto face_logical_coords = + interface_logical_coordinates(face_mesh, boundary_direction); + const auto face_inv_jacobian = + logical_to_inertial_map.inv_jacobian(face_logical_coords); + + // Slice required fields to the interface + const size_t slice_index = + index_to_slice_at(mesh.extents(), boundary_direction); + const auto& face_conformal_factor = + data_on_slice(conformal_factor, mesh.extents(), + boundary_direction.dimension(), slice_index); + const auto& face_unit_normal = + data_on_slice(unit_normal, mesh.extents(), + boundary_direction.dimension(), slice_index); + + // Compute Euclidean area element + const auto area_element = + euclidean_area_element(face_inv_jacobian, boundary_direction); + + // Compute Euclidean face normal + auto euclidean_face_normal = unnormalized_face_normal( + face_mesh, logical_to_inertial_map, boundary_direction); + const auto face_normal_magnitude = magnitude(euclidean_face_normal); + for (size_t d = 0; d < 3; ++d) { + euclidean_face_normal.get(d) /= get(face_normal_magnitude); + } + + // Evaluate surface integral. + const auto surface_integrand = Xcts::center_of_mass_surface_integrand( + face_conformal_factor, face_unit_normal); + const auto contracted_integrand = tenex::evaluate( + -surface_integrand(ti::I, ti::J) * euclidean_face_normal(ti::j)); + for (int I = 0; I < 3; I++) { + total_integral.get(I) += definite_integral( + contracted_integrand.get(I) * get(area_element), face_mesh); + } + } + } + + // Check result + auto custom_approx = Approx::custom().epsilon(TOLERANCE).scale(1.0); + CHECK(get<0>(total_integral) == custom_approx(0.)); + CHECK(get<1>(total_integral) == custom_approx(0.)); + CHECK(get<2>(total_integral) / mass == custom_approx(z_shift)); +} + } // namespace SPECTRE_TEST_CASE("Unit.PointwiseFunctions.Xcts.CenterOfMass", "[Unit][PointwiseFunctions]") { - // Test that integral converges with distance. - for (const double distance : std::array({1.e3, 1.e4, 1.e5})) { - test_center_of_mass_surface_integral(distance); + for (const double distance : std::array({1.e4, 1.e5, 1.e6})) { + test_infinite_surface_integral(distance, 0.1); + // Note: the infinite volume integral blows up for non-zero shifts. + test_infinite_volume_integral(distance, 0.); } }