From f7a04b0e7d32c130c20fbc6efeba3ea543293b9f 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 | 30 +- .../Xcts/Events/ObserveAdmIntegrals.hpp | 14 +- .../Xcts/AdmLinearMomentum.cpp | 23 +- .../Xcts/AdmLinearMomentum.hpp | 55 ++-- src/PointwiseFunctions/Xcts/AdmMass.cpp | 63 ++++ src/PointwiseFunctions/Xcts/AdmMass.hpp | 110 ++++++- src/PointwiseFunctions/Xcts/CenterOfMass.cpp | 35 ++- src/PointwiseFunctions/Xcts/CenterOfMass.hpp | 69 ++++- .../Xcts/Events/Test_ObserveAdmIntegrals.cpp | 17 +- .../Xcts/Test_AdmLinearMomentum.cpp | 286 ++++++++++++++---- .../PointwiseFunctions/Xcts/Test_AdmMass.cpp | 265 ++++++++++++++-- .../Xcts/Test_CenterOfMass.cpp | 180 ++++++++--- 12 files changed, 928 insertions(+), 219 deletions(-) diff --git a/src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.cpp b/src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.cpp index c456903e0ff1..1a457745f76d 100644 --- a/src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.cpp +++ b/src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.cpp @@ -31,12 +31,11 @@ void local_adm_integrals( const tnsr::II& inv_spatial_metric, const tnsr::ii& extrinsic_curvature, const Scalar& trace_extrinsic_curvature, + const tnsr::I& inertial_coords, const InverseJacobian& inv_jacobian, const Mesh<3>& mesh, const Element<3>& element, - const DirectionMap<3, tnsr::i>& conformal_face_normals, - const DirectionMap<3, tnsr::I>& - conformal_face_normal_vectors) { + const DirectionMap<3, tnsr::i>& conformal_face_normals) { // Initialize integrals to 0 adm_mass->get() = 0; for (int I = 0; I < 3; I++) { @@ -82,6 +81,8 @@ void local_adm_integrals( extrinsic_curvature, mesh, boundary_direction); const auto face_trace_extrinsic_curvature = dg::project_tensor_to_boundary( trace_extrinsic_curvature, mesh, boundary_direction); + const auto face_inertial_coords = dg::project_tensor_to_boundary( + inertial_coords, mesh, boundary_direction); // This projection could be avoided by using // domain::Tags::DetSurfaceJacobian from the DataBox, which is computed // directly on the face (not projected). That would be better on Gauss @@ -103,15 +104,18 @@ void local_adm_integrals( const auto& face_mesh = mesh.slice_away(boundary_direction.dimension()); const auto& conformal_face_normal = conformal_face_normals.at(boundary_direction); - const auto& conformal_face_normal_vector = - conformal_face_normal_vectors.at(boundary_direction); + const auto face_normal_magnitude = magnitude(conformal_face_normal); + const auto flat_face_normal = tenex::evaluate( + conformal_face_normal(ti::i) / face_normal_magnitude()); - // Compute curved area element + // Compute curved and flat area elements const auto face_sqrt_det_conformal_metric = Scalar(sqrt(get(determinant(face_conformal_metric)))); const auto conformal_area_element = area_element(face_inv_jacobian, boundary_direction, face_inv_conformal_metric, face_sqrt_det_conformal_metric); + const auto flat_area_element = + euclidean_area_element(face_inv_jacobian, boundary_direction); // Compute surface integrands const auto mass_integrand = Xcts::adm_mass_surface_integrand( @@ -124,26 +128,24 @@ void local_adm_integrals( face_inv_extrinsic_curvature, face_trace_extrinsic_curvature); const auto center_of_mass_integrand = Xcts::center_of_mass_surface_integrand(face_conformal_factor, - conformal_face_normal_vector); + face_inertial_coords); // Contract surface integrands with face normal const auto contracted_mass_integrand = 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)); + linear_momentum_integrand(ti::I, ti::J) * flat_face_normal(ti::j)); // Take integrals adm_mass->get() += definite_integral( get(contracted_mass_integrand) * get(conformal_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(conformal_area_element), - face_mesh); - center_of_mass->get(I) += definite_integral( - center_of_mass_integrand.get(I) * get(conformal_area_element), + adm_linear_momentum->get(I) += definite_integral( + contracted_linear_momentum_integrand.get(I) * get(flat_area_element), face_mesh); + center_of_mass->get(I) += definite_integral( + center_of_mass_integrand.get(I) * get(flat_area_element), face_mesh); } } } diff --git a/src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.hpp b/src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.hpp index 36d6eb07843f..d0d91eeabb8a 100644 --- a/src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.hpp +++ b/src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.hpp @@ -56,12 +56,11 @@ void local_adm_integrals( const tnsr::II& inv_spatial_metric, const tnsr::ii& extrinsic_curvature, const Scalar& trace_extrinsic_curvature, + const tnsr::I& inertial_coords, const InverseJacobian& inv_jacobian, const Mesh<3>& mesh, const Element<3>& element, - const DirectionMap<3, tnsr::i>& conformal_face_normals, - const DirectionMap<3, tnsr::I>& - conformal_face_normal_vectors); + const DirectionMap<3, tnsr::i>& conformal_face_normals); /// @} /// @{ @@ -135,10 +134,10 @@ class ObserveAdmIntegrals : public Event { gr::Tags::InverseSpatialMetric, gr::Tags::ExtrinsicCurvature, gr::Tags::TraceExtrinsicCurvature, + domain::Tags::Coordinates<3, Frame::Inertial>, domain::Tags::InverseJacobian<3, Frame::ElementLogical, Frame::Inertial>, domain::Tags::Mesh<3>, domain::Tags::Element<3>, domain::Tags::Faces<3, domain::Tags::FaceNormal<3>>, - domain::Tags::Faces<3, domain::Tags::FaceNormalVector<3>>, ::Tags::ObservationBox>; template & inv_spatial_metric, const tnsr::ii& extrinsic_curvature, const Scalar& trace_extrinsic_curvature, + const tnsr::I& inertial_coords, const InverseJacobian& inv_jacobian, const Mesh<3>& mesh, const Element<3>& element, const DirectionMap<3, tnsr::i>& conformal_face_normals, - const DirectionMap<3, tnsr::I>& - conformal_face_normal_vectors, const ObservationBox& box, Parallel::GlobalCache& cache, const ArrayIndex& array_index, const ParallelComponent* const /*meta*/, @@ -182,8 +180,8 @@ class ObserveAdmIntegrals : public Event { 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, inv_jacobian, mesh, element, - conformal_face_normals, conformal_face_normal_vectors); + trace_extrinsic_curvature, inertial_coords, inv_jacobian, mesh, element, + conformal_face_normals); // Save components of linear momentum as reduction data ReductionData reduction_data{get(adm_mass), 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..df45b5521a75 100644 --- a/src/PointwiseFunctions/Xcts/AdmLinearMomentum.hpp +++ b/src/PointwiseFunctions/Xcts/AdmLinearMomentum.hpp @@ -4,30 +4,30 @@ #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 + * \note For consistency with `adm_linear_momentum_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 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 +54,26 @@ 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. + * \note For consistency with `adm_linear_momentum_surface_integrand`, this + * integrand needs to be integrated with the Euclidean volume element. * - * \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 +86,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 +96,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..e13977b2c07c 100644 --- a/src/PointwiseFunctions/Xcts/AdmMass.cpp +++ b/src/PointwiseFunctions/Xcts/AdmMass.cpp @@ -33,4 +33,67 @@ 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 Scalar& energy_density, + const tnsr::II& inv_conformal_metric, + const tnsr::iJK& deriv_inv_conformal_metric, + const tnsr::Ijj& conformal_christoffel_second_kind, + const tnsr::i& conformal_christoffel_contracted, + const tnsr::iJkk& deriv_conformal_christoffel_second_kind) { + tenex::evaluate( + result, + 1. / (16. * M_PI) * + (deriv_inv_conformal_metric(ti::i, ti::J, ti::K) * + conformal_christoffel_second_kind(ti::I, ti::j, ti::k) + + inv_conformal_metric(ti::J, ti::K) * + deriv_conformal_christoffel_second_kind(ti::i, ti::I, ti::j, + ti::k) + + conformal_christoffel_contracted(ti::l) * + inv_conformal_metric(ti::J, ti::K) * + conformal_christoffel_second_kind(ti::L, ti::j, ti::k) - + deriv_inv_conformal_metric(ti::i, ti::I, ti::J) * + conformal_christoffel_contracted(ti::j) - + inv_conformal_metric(ti::I, ti::J) * + deriv_conformal_christoffel_second_kind(ti::i, ti::K, ti::k, + ti::j) - + conformal_christoffel_contracted(ti::l) * + inv_conformal_metric(ti::L, ti::J) * + conformal_christoffel_contracted(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() + + 16. * M_PI * pow<5>(conformal_factor()) * energy_density())); +} + +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 Scalar& energy_density, + const tnsr::II& inv_conformal_metric, + const tnsr::iJK& deriv_inv_conformal_metric, + const tnsr::Ijj& conformal_christoffel_second_kind, + const tnsr::i& conformal_christoffel_contracted, + 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, + energy_density, inv_conformal_metric, deriv_inv_conformal_metric, + conformal_christoffel_second_kind, conformal_christoffel_contracted, + 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..b7469d102b0d 100644 --- a/src/PointwiseFunctions/Xcts/AdmMass.hpp +++ b/src/PointwiseFunctions/Xcts/AdmMass.hpp @@ -16,20 +16,24 @@ 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) d\bar{S}_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 + * \note For consistency with `adm_mass_volume_integrand`, this integrand needs + * to be contracted with the conformal face normal and integrated with the + * conformal area element. + * + * \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 +56,88 @@ 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( + * \partial_i \bar\gamma^{jk} \bar\Gamma^i_{jk} + * + \bar\gamma^{jk} \partial_i \bar\Gamma^i_{jk} + * + \bar\Gamma_l \bar\gamma^{jk} \bar\Gamma^l_{jk} + * - \partial_i \bar\gamma^{ij} \bar\Gamma_j + * - \bar\gamma^{ij} \partial_i \bar\Gamma_j + * - \bar\Gamma_l \bar\gamma^{lj} \bar\Gamma_j + * - 8 \bar D^2 \psi + * \Big) d\bar{V}, + * \end{equation} + * + * where we can use the Hamiltonian constraint (Eq. 3.37 in + * \cite BaumgarteShapiro) 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] + * - 16\pi \psi^5 \rho. + * \end{equation} + * + * \note This is similar to Eq. 3.149 in \cite BaumgarteShapiro, except that + * here we don't assume $\bar\gamma = 1$. + * + * \note For consistency with `adm_mass_surface_integrand`, this integrand needs + * to be integrated with the conformal volume element. + * + * \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 energy_density the energy density $\rho$ + * \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 conformal_christoffel_contracted the conformal christoffel symbol + * contracted in its first two indices $\bar\Gamma_{i} = \bar\Gamma^j_{ij}$ + * \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 Scalar& energy_density, + const tnsr::II& inv_conformal_metric, + const tnsr::iJK& deriv_inv_conformal_metric, + const tnsr::Ijj& conformal_christoffel_second_kind, + const tnsr::i& conformal_christoffel_contracted, + 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 Scalar& energy_density, + const tnsr::II& inv_conformal_metric, + const tnsr::iJK& deriv_inv_conformal_metric, + const tnsr::Ijj& conformal_christoffel_second_kind, + const tnsr::i& conformal_christoffel_contracted, + 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..15784cdd2b3a 100644 --- a/src/PointwiseFunctions/Xcts/CenterOfMass.cpp +++ b/src/PointwiseFunctions/Xcts/CenterOfMass.cpp @@ -3,22 +3,49 @@ #include "PointwiseFunctions/Xcts/CenterOfMass.hpp" +#include "DataStructures/Tensor/EagerMath/Magnitude.hpp" + namespace Xcts { void center_of_mass_surface_integrand( gsl::not_null*> result, const Scalar& conformal_factor, - const tnsr::I& unit_normal) { + const tnsr::I& coords) { + const auto euclidean_radius = magnitude(coords); tenex::evaluate(result, 3. / (8. * M_PI) * pow<4>(conformal_factor()) * - unit_normal(ti::I)); + coords(ti::I) / euclidean_radius()); } tnsr::I center_of_mass_surface_integrand( const Scalar& conformal_factor, - const tnsr::I& unit_normal) { + const tnsr::I& coords) { tnsr::I result; center_of_mass_surface_integrand(make_not_null(&result), conformal_factor, - unit_normal); + coords); + 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& coords) { + const auto euclidean_radius = magnitude(coords); + tenex::evaluate( + result, + 3. / (4. * M_PI * pow<2>(euclidean_radius())) * + (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 center_of_mass_volume_integrand( + const Scalar& conformal_factor, + const tnsr::i& deriv_conformal_factor, + const tnsr::I& coords) { + tnsr::I result; + center_of_mass_volume_integrand(make_not_null(&result), conformal_factor, + deriv_conformal_factor, coords); return result; } diff --git a/src/PointwiseFunctions/Xcts/CenterOfMass.hpp b/src/PointwiseFunctions/Xcts/CenterOfMass.hpp index 24c49f1a9608..ed913d725720 100644 --- a/src/PointwiseFunctions/Xcts/CenterOfMass.hpp +++ b/src/PointwiseFunctions/Xcts/CenterOfMass.hpp @@ -17,14 +17,19 @@ 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 \, dA, * \end{equation} * - * Note that we don't include the ADM mass $M_{ADM}$ in this integrand. After + * where $n^i = x^i / r$ and $r = \sqrt{x^2 + y^2 + z^2}$. + * + * \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}$. + * + * \note For consistency with `center_of_mass_volume_integrand`, this + * integrand needs to be integrated with the Euclidean area element. + * + * \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,19 +37,65 @@ 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$ + * \param coords the inertial coordinates $x^i$ */ void center_of_mass_surface_integrand( gsl::not_null*> result, const Scalar& conformal_factor, - const tnsr::I& unit_normal); + const tnsr::I& coords); /// Return-by-value overload tnsr::I center_of_mass_surface_integrand( const Scalar& conformal_factor, - const tnsr::I& unit_normal); + const tnsr::I& coords); +/// @} + +/// @{ +/*! + * \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} \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}$. + * + * \note For consistency with `center_of_mass_surface_integrand`, this + * integrand needs to be integrated with the Euclidean volume element. + * + * \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 coords the inertial coordinates $x^i$ + */ +void center_of_mass_volume_integrand( + gsl::not_null*> result, + const Scalar& conformal_factor, + const tnsr::i& deriv_conformal_factor, + const tnsr::I& coords); + +/// Return-by-value overload +tnsr::I center_of_mass_volume_integrand( + const Scalar& conformal_factor, + const tnsr::i& deriv_conformal_factor, + const tnsr::I& coords); /// @} } // namespace Xcts diff --git a/tests/Unit/Elliptic/Systems/Xcts/Events/Test_ObserveAdmIntegrals.cpp b/tests/Unit/Elliptic/Systems/Xcts/Events/Test_ObserveAdmIntegrals.cpp index 0393adae01f0..14a0080df979 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); @@ -203,7 +204,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 = @@ -242,12 +243,6 @@ void test_local_adm_integrals(const double& distance, const DirectionMap<3, tnsr::i> conformal_face_normals( {std::make_pair(direction, conformal_face_normal)}); - // Compute face normal vector. - const auto conformal_face_normal_vector = tenex::evaluate( - face_inv_conformal_metric(ti::I, ti::J) * conformal_face_normal(ti::j)); - const DirectionMap<3, tnsr::I> conformal_face_normal_vectors( - {std::make_pair(direction, conformal_face_normal_vector)}); - // Compute local integrals. Scalar local_adm_mass; tnsr::I local_adm_linear_momentum; @@ -259,8 +254,8 @@ void test_local_adm_integrals(const double& distance, 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, inv_jacobian, mesh, current_element, - conformal_face_normals, conformal_face_normal_vectors); + trace_extrinsic_curvature, inertial_coords, inv_jacobian, mesh, + current_element, conformal_face_normals); total_adm_mass.get() += get(local_adm_mass); for (int I = 0; I < 3; I++) { total_adm_linear_momentum.get(I) += local_adm_linear_momentum.get(I); diff --git a/tests/Unit/PointwiseFunctions/Xcts/Test_AdmLinearMomentum.cpp b/tests/Unit/PointwiseFunctions/Xcts/Test_AdmLinearMomentum.cpp index f8835f237812..0f455943c6cf 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{}, @@ -95,20 +95,11 @@ void test_linear_momentum_surface_integral(const double distance, inertial_coords, tmpl::list< Xcts::Tags::ConformalFactor, - Xcts::Tags::ConformalMetric, - Xcts::Tags::InverseConformalMetric, gr::Tags::InverseSpatialMetric, gr::Tags::ExtrinsicCurvature, gr::Tags::TraceExtrinsicCurvature>{}); const auto& conformal_factor = get>(background_fields); - const auto& conformal_metric = - get>( - background_fields); - const auto& inv_conformal_metric = get< - Xcts::Tags::InverseConformalMetric>( - background_fields); const auto& inv_spatial_metric = get>( background_fields); @@ -125,35 +116,28 @@ void test_linear_momentum_surface_integral(const double distance, inv_spatial_metric(ti::J, ti::L) * extrinsic_curvature(ti::k, ti::l)); - // Compute conformal area element - const auto sqrt_det_conformal_metric = - Scalar(sqrt(get(determinant(conformal_metric)))); - const auto conformal_area_element = - area_element(inv_jacobian, boundary_direction, inv_conformal_metric, - sqrt_det_conformal_metric); + // Compute Euclidean area element + const auto flat_area_element = + euclidean_area_element(inv_jacobian, boundary_direction); - // Compute conformal face normal - auto conformal_face_normal = unnormalized_face_normal( + // Compute Euclidean face normal + auto flat_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); + const auto face_normal_magnitude = magnitude(flat_face_normal); for (size_t d = 0; d < 3; ++d) { - conformal_face_normal.get(d) /= get(face_normal_magnitude); + flat_face_normal.get(d) /= get(face_normal_magnitude); } - // Compute surface integrand + // Evaluate surface integral const auto surface_integrand = Xcts::adm_linear_momentum_surface_integrand( conformal_factor, inv_spatial_metric, inv_extrinsic_curvature, trace_extrinsic_curvature); const auto contracted_integrand = tenex::evaluate( - surface_integrand(ti::I, ti::J) * conformal_face_normal(ti::j)); - - // Compute contribution to surface integral + surface_integrand(ti::I, ti::J) * flat_face_normal(ti::j)); for (int I = 0; I < 3; I++) { surface_integral.get(I) += definite_integral( - contracted_integrand.get(I) * get(conformal_area_element), - face_mesh); + contracted_integrand.get(I) * get(flat_area_element), face_mesh); } } } @@ -167,10 +151,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 flat_area_element = + euclidean_area_element(face_inv_jacobian, boundary_direction); + + // Compute Euclidean face normal + auto flat_face_normal = unnormalized_face_normal( + face_mesh, logical_to_inertial_map, boundary_direction); + const auto face_normal_magnitude = magnitude(flat_face_normal); + for (size_t d = 0; d < 3; ++d) { + flat_face_normal.get(d) /= get(face_normal_magnitude); + } + + // Compute surface integrand + const auto contracted_integrand = tenex::evaluate( + -face_surface_integrand(ti::I, ti::J) * flat_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(flat_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 +387,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..a7655108991f 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,256 @@ 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)); + + const auto energy_density = + make_with_value>(inertial_coords, 0.0); + + const auto sqrt_det_conformal_metric = + Scalar(sqrt(get(determinant(conformal_metric)))); + + // 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, + energy_density, inv_conformal_metric, deriv_inv_conformal_metric, + conformal_christoffel_second_kind, conformal_christoffel_contracted, + deriv_conformal_christoffel_second_kind); + total_integral.get() += definite_integral( + get(volume_integrand) * get(sqrt_det_conformal_metric) * + 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); + const auto& face_sqrt_det_conformal_metric = + data_on_slice(sqrt_det_conformal_metric, mesh.extents(), + boundary_direction.dimension(), slice_index); + + // Compute conformal area element + const auto conformal_area_element = area_element( + face_inv_jacobian, boundary_direction, face_inv_conformal_metric, + face_sqrt_det_conformal_metric); + + // Compute conformal 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, face_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::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) * conformal_face_normal(ti::i)); + + // Compute contribution to surface integral + total_integral.get() += definite_integral( + get(contracted_integrand) * get(conformal_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 +// [[TimeOut, 60]] 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..c00f3ae7b461 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-3; + +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, @@ -99,55 +106,156 @@ void test_center_of_mass_surface_integral(const double distance) { // Get required fields on the interface const auto shifted_fields = solution.variables( shifted_coords, - tmpl::list< - Xcts::Tags::ConformalFactor, - Xcts::Tags::ConformalMetric, - Xcts::Tags::InverseConformalMetric>{}); + tmpl::list>{}); const auto& conformal_factor = get>(shifted_fields); - const auto& conformal_metric = - get>( - shifted_fields); - const auto& inv_conformal_metric = get< - Xcts::Tags::InverseConformalMetric>( - shifted_fields); - - // Compute outward-pointing unit normal. - const auto conformal_r = magnitude(inertial_coords, conformal_metric); - const auto conformal_unit_normal = - tenex::evaluate(inertial_coords(ti::I) / conformal_r()); // Compute area element - const auto sqrt_det_conformal_metric = - Scalar(sqrt(get(determinant(conformal_metric)))); - const auto conformal_area_element = - area_element(inv_jacobian, boundary_direction, inv_conformal_metric, - sqrt_det_conformal_metric); + const auto flat_area_element = + euclidean_area_element(inv_jacobian, boundary_direction); - // Integrate + // Evaluate surface integral const auto surface_integrand = Xcts::center_of_mass_surface_integrand( - conformal_factor, conformal_unit_normal); + conformal_factor, inertial_coords); for (int I = 0; I < 3; I++) { surface_integral.get(I) += definite_integral( - surface_integrand.get(I) * get(conformal_area_element), face_mesh); + surface_integrand.get(I) * get(flat_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 = 11; + 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); + + // 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); + + // Evaluate volume integral. + const auto volume_integrand = Xcts::center_of_mass_volume_integrand( + 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); + } + + // 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_inertial_coords = + data_on_slice(inertial_coords, mesh.extents(), + boundary_direction.dimension(), slice_index); + + // Compute Euclidean area element + const auto flat_area_element = + euclidean_area_element(face_inv_jacobian, boundary_direction); + + // Evaluate surface integral. + const auto surface_integrand = Xcts::center_of_mass_surface_integrand( + face_conformal_factor, face_inertial_coords); + for (int I = 0; I < 3; I++) { + total_integral.get(I) += definite_integral( + surface_integrand.get(I) * get(flat_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); + test_infinite_surface_integral(distance, 0.); + test_infinite_surface_integral(distance, 0.1); + test_infinite_volume_integral(distance, 0.); + test_infinite_volume_integral(distance, 0.1); } }