Skip to content

Commit

Permalink
Merge pull request #6276 from iago-mendes/volume-integrands
Browse files Browse the repository at this point in the history 
Add volume integrands for XCTS asymptotic quantities.
  • Loading branch information
@nilsvu
nilsvu authored Sep 23, 2024
2 parents 32dd542 + f7a04b0 commit 20b3671
Show file tree
Hide file tree
Showing 12 changed files with 928 additions and 219 deletions.
30 changes: 16 additions & 14 deletions src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,12 +31,11 @@ void local_adm_integrals(
const tnsr::II<DataVector, 3>& inv_spatial_metric,
const tnsr::ii<DataVector, 3>& extrinsic_curvature,
const Scalar<DataVector>& trace_extrinsic_curvature,
const tnsr::I<DataVector, 3, Frame::Inertial>& inertial_coords,
const InverseJacobian<DataVector, 3, Frame::ElementLogical,
Frame::Inertial>& inv_jacobian,
const Mesh<3>& mesh, const Element<3>& element,
const DirectionMap<3, tnsr::i<DataVector, 3>>& conformal_face_normals,
const DirectionMap<3, tnsr::I<DataVector, 3>>&
conformal_face_normal_vectors) {
const DirectionMap<3, tnsr::i<DataVector, 3>>& conformal_face_normals) {
// Initialize integrals to 0
adm_mass->get() = 0;
for (int I = 0; I < 3; I++) {
Expand Down Expand Up @@ -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
Expand All @@ -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<ti::i>(
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<DataVector>(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(
Expand All @@ -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<ti::I>(
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);
}
}
}
Expand Down
14 changes: 6 additions & 8 deletions src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,12 +56,11 @@ void local_adm_integrals(
const tnsr::II<DataVector, 3>& inv_spatial_metric,
const tnsr::ii<DataVector, 3>& extrinsic_curvature,
const Scalar<DataVector>& trace_extrinsic_curvature,
const tnsr::I<DataVector, 3, Frame::Inertial>& inertial_coords,
const InverseJacobian<DataVector, 3, Frame::ElementLogical,
Frame::Inertial>& inv_jacobian,
const Mesh<3>& mesh, const Element<3>& element,
const DirectionMap<3, tnsr::i<DataVector, 3>>& conformal_face_normals,
const DirectionMap<3, tnsr::I<DataVector, 3>>&
conformal_face_normal_vectors);
const DirectionMap<3, tnsr::i<DataVector, 3>>& conformal_face_normals);
/// @}

/// @{
Expand Down Expand Up @@ -135,10 +134,10 @@ class ObserveAdmIntegrals : public Event {
gr::Tags::InverseSpatialMetric<DataVector, 3, Frame::Inertial>,
gr::Tags::ExtrinsicCurvature<DataVector, 3, Frame::Inertial>,
gr::Tags::TraceExtrinsicCurvature<DataVector>,
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 <typename DataBoxType, typename ComputeTagsList,
Expand All @@ -155,12 +154,11 @@ class ObserveAdmIntegrals : public Event {
const tnsr::II<DataVector, 3>& inv_spatial_metric,
const tnsr::ii<DataVector, 3>& extrinsic_curvature,
const Scalar<DataVector>& trace_extrinsic_curvature,
const tnsr::I<DataVector, 3, Frame::Inertial>& inertial_coords,
const InverseJacobian<DataVector, 3, Frame::ElementLogical,
Frame::Inertial>& inv_jacobian,
const Mesh<3>& mesh, const Element<3>& element,
const DirectionMap<3, tnsr::i<DataVector, 3>>& conformal_face_normals,
const DirectionMap<3, tnsr::I<DataVector, 3>>&
conformal_face_normal_vectors,
const ObservationBox<DataBoxType, ComputeTagsList>& box,
Parallel::GlobalCache<Metavariables>& cache,
const ArrayIndex& array_index, const ParallelComponent* const /*meta*/,
Expand All @@ -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),
Expand Down
23 changes: 6 additions & 17 deletions src/PointwiseFunctions/Xcts/AdmLinearMomentum.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,19 +3,6 @@

#include "PointwiseFunctions/Xcts/AdmLinearMomentum.hpp"

#include <cmath>

#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(
Expand Down Expand Up @@ -47,34 +34,36 @@ void adm_linear_momentum_volume_integrand(
gsl::not_null<tnsr::I<DataVector, 3>*> result,
const tnsr::II<DataVector, 3>& surface_integrand,
const Scalar<DataVector>& conformal_factor,
const tnsr::i<DataVector, 3>& conformal_factor_deriv,
const tnsr::i<DataVector, 3>& deriv_conformal_factor,
const tnsr::ii<DataVector, 3>& conformal_metric,
const tnsr::II<DataVector, 3>& inv_conformal_metric,
const tnsr::Ijj<DataVector, 3>& conformal_christoffel_second_kind,
const tnsr::i<DataVector, 3>& conformal_christoffel_contracted) {
// Note: we can ignore the $1/(8\pi)$ term below because it is already
// included in `surface_integrand`.
tenex::evaluate<ti::I>(
result,
-(conformal_christoffel_second_kind(ti::I, ti::j, ti::k) *
surface_integrand(ti::J, ti::K) +
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<DataVector, 3> adm_linear_momentum_volume_integrand(
const tnsr::II<DataVector, 3>& surface_integrand,
const Scalar<DataVector>& conformal_factor,
const tnsr::i<DataVector, 3>& conformal_factor_deriv,
const tnsr::i<DataVector, 3>& deriv_conformal_factor,
const tnsr::ii<DataVector, 3>& conformal_metric,
const tnsr::II<DataVector, 3>& inv_conformal_metric,
const tnsr::Ijj<DataVector, 3>& conformal_christoffel_second_kind,
const tnsr::i<DataVector, 3>& conformal_christoffel_contracted) {
tnsr::I<DataVector, 3> 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;
}
Expand Down
55 changes: 29 additions & 26 deletions src/PointwiseFunctions/Xcts/AdmLinearMomentum.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 <cmath>

#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}$
Expand All @@ -54,23 +54,26 @@ tnsr::II<DataVector, 3> 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}$
Expand All @@ -83,7 +86,7 @@ void adm_linear_momentum_volume_integrand(
gsl::not_null<tnsr::I<DataVector, 3>*> result,
const tnsr::II<DataVector, 3>& surface_integrand,
const Scalar<DataVector>& conformal_factor,
const tnsr::i<DataVector, 3>& conformal_factor_deriv,
const tnsr::i<DataVector, 3>& deriv_conformal_factor,
const tnsr::ii<DataVector, 3>& conformal_metric,
const tnsr::II<DataVector, 3>& inv_conformal_metric,
const tnsr::Ijj<DataVector, 3>& conformal_christoffel_second_kind,
Expand All @@ -93,7 +96,7 @@ void adm_linear_momentum_volume_integrand(
tnsr::I<DataVector, 3> adm_linear_momentum_volume_integrand(
const tnsr::II<DataVector, 3>& surface_integrand,
const Scalar<DataVector>& conformal_factor,
const tnsr::i<DataVector, 3>& conformal_factor_deriv,
const tnsr::i<DataVector, 3>& deriv_conformal_factor,
const tnsr::ii<DataVector, 3>& conformal_metric,
const tnsr::II<DataVector, 3>& inv_conformal_metric,
const tnsr::Ijj<DataVector, 3>& conformal_christoffel_second_kind,
Expand Down
63 changes: 63 additions & 0 deletions src/PointwiseFunctions/Xcts/AdmMass.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,4 +33,67 @@ tnsr::I<DataVector, 3> adm_mass_surface_integrand(
return result;
}

void adm_mass_volume_integrand(
gsl::not_null<Scalar<DataVector>*> result,
const Scalar<DataVector>& conformal_factor,
const Scalar<DataVector>& conformal_ricci_scalar,
const Scalar<DataVector>& trace_extrinsic_curvature,
const Scalar<DataVector>&
longitudinal_shift_minus_dt_conformal_metric_over_lapse_square,
const Scalar<DataVector>& energy_density,
const tnsr::II<DataVector, 3>& inv_conformal_metric,
const tnsr::iJK<DataVector, 3>& deriv_inv_conformal_metric,
const tnsr::Ijj<DataVector, 3>& conformal_christoffel_second_kind,
const tnsr::i<DataVector, 3>& conformal_christoffel_contracted,
const tnsr::iJkk<DataVector, 3>& 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<DataVector> adm_mass_volume_integrand(
const Scalar<DataVector>& conformal_factor,
const Scalar<DataVector>& conformal_ricci_scalar,
const Scalar<DataVector>& trace_extrinsic_curvature,
const Scalar<DataVector>&
longitudinal_shift_minus_dt_conformal_metric_over_lapse_square,
const Scalar<DataVector>& energy_density,
const tnsr::II<DataVector, 3>& inv_conformal_metric,
const tnsr::iJK<DataVector, 3>& deriv_inv_conformal_metric,
const tnsr::Ijj<DataVector, 3>& conformal_christoffel_second_kind,
const tnsr::i<DataVector, 3>& conformal_christoffel_contracted,
const tnsr::iJkk<DataVector, 3>& deriv_conformal_christoffel_second_kind) {
Scalar<DataVector> 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
Loading

0 comments on commit 20b3671

Please sign in to comment.