Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add volume integrands for XCTS asymptotic quantities. #6276

Merged
merged 1 commit into from
Sep 23, 2024
Merged

Add volume integrands for XCTS asymptotic quantities. #6276

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
30 changes: 16 additions & 14 deletions src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.cpp
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading