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

Compute ADM momenta as volume integrals #6433

Open
wants to merge 2 commits into
base: develop
Choose a base branch
from
+121 −77
Open

Compute ADM momenta as volume integrals #6433

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
c656ebb
Compute ADM momenta as volume integrals
iago-mendes Oct 10, 2024
fe89e58
fixup
iago-mendes Jan 24, 2025
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
162 changes: 98 additions & 64 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 @@ -45,19 +45,42 @@ void local_adm_integrals(
center_of_mass->get(I) = 0;
}

// Skip elements not at the outer boundary
const size_t radial_dimension = 2;
const auto radial_segment_id = element.id().segment_id(radial_dimension);
if (radial_segment_id.index() !=
two_to_the(radial_segment_id.refinement_level()) - 1) {
return;
// Compute quantities not in the data box
tnsr::II<DataVector, 3> inv_extrinsic_curvature;
tenex::evaluate<ti::I, ti::J>(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));

// Get integrands
const auto linear_momentum_surface_integrand =
Xcts::adm_linear_momentum_surface_integrand(
conformal_factor, inv_spatial_metric, inv_extrinsic_curvature,
trace_extrinsic_curvature);
const auto linear_momentum_volume_integrand =
Xcts::adm_linear_momentum_volume_integrand(
linear_momentum_surface_integrand, conformal_factor,
deriv_conformal_factor, conformal_metric, inv_conformal_metric,
conformal_christoffel_second_kind, conformal_christoffel_contracted);
const auto angular_momentum_z_volume_integrand =
Xcts::adm_angular_momentum_z_volume_integrand(
linear_momentum_volume_integrand, inertial_coords);

// Evaluate volume integrals.
const auto det_jacobian =
Scalar<DataVector>(1. / get(determinant(inv_jacobian)));
adm_angular_momentum_z->get() += definite_integral(
get(angular_momentum_z_volume_integrand) * get(det_jacobian), mesh);
for (int I = 0; I < 3; I++) {
adm_linear_momentum->get(I) += definite_integral(
linear_momentum_volume_integrand.get(I) * get(det_jacobian), mesh);
}

// Loop over external boundaries
for (const auto boundary_direction : element.external_boundaries()) {
// Skip interfaces not at the outer boundary
if (boundary_direction != Direction<3>::upper_zeta()) {
continue;
// Skip non-zeta boundaries
if (boundary_direction.dimension() != 2) {
continue;
}

// Project fields to the boundary
Expand All @@ -70,37 +93,31 @@ void local_adm_integrals(
const auto face_inv_conformal_metric = dg::project_tensor_to_boundary(
inv_conformal_metric, mesh, boundary_direction);
const auto face_conformal_christoffel_second_kind =
dg::project_tensor_to_boundary(conformal_christoffel_second_kind, mesh,
boundary_direction);
dg::project_tensor_to_boundary(conformal_christoffel_second_kind,
mesh, boundary_direction);
const auto face_conformal_christoffel_contracted =
dg::project_tensor_to_boundary(conformal_christoffel_contracted, mesh,
boundary_direction);
boundary_direction);
const auto face_spatial_metric = dg::project_tensor_to_boundary(
spatial_metric, mesh, boundary_direction);
const auto face_inv_spatial_metric = dg::project_tensor_to_boundary(
inv_spatial_metric, mesh, boundary_direction);
const auto face_extrinsic_curvature = dg::project_tensor_to_boundary(
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_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
// meshes that have no grid point at the boundary. The DetSurfaceJacobian
// could then be multiplied by the (conformal) metric determinant to form
// the area element. Note that the DetSurfaceJacobian is computed using the
// conformal metric.
const auto face_inv_jacobian =
dg::project_tensor_to_boundary(inv_jacobian, mesh, boundary_direction);

// Compute the inverse extrinsic curvature
tnsr::II<DataVector, 3> face_inv_extrinsic_curvature;
tenex::evaluate<ti::I, ti::J>(make_not_null(&face_inv_extrinsic_curvature),
face_inv_spatial_metric(ti::I, ti::K) *
face_inv_spatial_metric(ti::J, ti::L) *
face_extrinsic_curvature(ti::k, ti::l));
// the area element. Note that the DetSurfaceJacobian is computed using
// the conformal metric.
const auto face_inv_jacobian = dg::project_tensor_to_boundary(
inv_jacobian, mesh, boundary_direction);

// Get interface mesh and normal
const auto& face_mesh = mesh.slice_away(boundary_direction.dimension());
Expand All @@ -113,49 +130,66 @@ void local_adm_integrals(
// 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 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(
face_deriv_conformal_factor, face_inv_conformal_metric,
face_conformal_christoffel_second_kind,
face_conformal_christoffel_contracted);
const auto linear_momentum_integrand =
Xcts::adm_linear_momentum_surface_integrand(
face_conformal_factor, face_inv_spatial_metric,
face_inv_extrinsic_curvature, face_trace_extrinsic_curvature);
const auto angular_momentum_z_integrand =
Xcts::adm_angular_momentum_z_surface_integrand(
linear_momentum_integrand, face_inertial_coords);
const auto center_of_mass_integrand =
Xcts::center_of_mass_surface_integrand(face_conformal_factor,
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) * flat_face_normal(ti::j));
const auto contracted_angular_momentum_z_integrand = tenex::evaluate(
angular_momentum_z_integrand(ti::I) * flat_face_normal(ti::i));

// Take integrals
adm_mass->get() += definite_integral(
get(contracted_mass_integrand) * get(conformal_area_element),
face_mesh);
adm_angular_momentum_z->get() += definite_integral(
get(contracted_angular_momentum_z_integrand) * get(flat_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(flat_area_element),
// Interfaces at the inner boundary
if (boundary_direction == Direction<3>::lower_zeta()) {
// Compute surface integrands
const auto face_linear_momentum_surface_integrand =
dg::project_tensor_to_boundary(linear_momentum_surface_integrand,
mesh, boundary_direction);
const auto angular_momentum_z_surface_integrand =
Xcts::adm_angular_momentum_z_surface_integrand(
face_linear_momentum_surface_integrand, face_inertial_coords);

// Contract surface integrands with face normal
const auto contracted_linear_momentum_integrand = tenex::evaluate<ti::I>(
-face_linear_momentum_surface_integrand(ti::I, ti::J) *
flat_face_normal(ti::j));
const auto contracted_angular_momentum_z_integrand =
tenex::evaluate(-angular_momentum_z_surface_integrand(ti::I) *
flat_face_normal(ti::i));

// Take integrals
adm_angular_momentum_z->get() += definite_integral(
get(contracted_angular_momentum_z_integrand) * get(flat_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(flat_area_element),
face_mesh);
}
}

// Interfaces at the outer boundary
if (boundary_direction == Direction<3>::upper_zeta()) {
// Compute surface integrands
const auto mass_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 center_of_mass_surface_integrand =
Xcts::center_of_mass_surface_integrand(face_conformal_factor,
face_inertial_coords);

// Contract surface integrands with face normal
const auto contracted_mass_integrand = tenex::evaluate(
mass_surface_integrand(ti::I) * conformal_face_normal(ti::i));

// Take integrals
adm_mass->get() += definite_integral(
get(contracted_mass_integrand) * get(conformal_area_element),
face_mesh);
center_of_mass->get(I) += definite_integral(
center_of_mass_integrand.get(I) * get(flat_area_element), face_mesh);
for (int I = 0; I < 3; I++) {
center_of_mass->get(I) += definite_integral(
center_of_mass_surface_integrand.get(I) * get(flat_area_element),
face_mesh);
}
}
}
}
Expand Down
36 changes: 23 additions & 13 deletions tests/Unit/Elliptic/Systems/Xcts/Events/Test_ObserveAdmIntegrals.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,13 +25,11 @@
#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 {

Expand All @@ -54,7 +52,7 @@ void test_local_adm_integrals(const double& distance,
const std::vector<double>& prev_distances) {
// Define black hole parameters.
const double mass = 1;
const double boost_speed = 0.5;
const double boost_speed = 0.;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this change needed?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, because the volume integrals give incorrect results for a boosted Schwarzschild case. In the past, we've wondered if this is due to some falloff assumption from Gauss Law being broken.

const double lorentz_factor = 1. / sqrt(1. - square(boost_speed));
const std::array<double, 3> boost_velocity{{0., 0., boost_speed}};

Expand Down Expand Up @@ -231,19 +229,30 @@ void test_local_adm_integrals(const double& distance,
const auto trace_extrinsic_curvature = tenex::evaluate(
inv_spatial_metric(ti::I, ti::J) * extrinsic_curvature(ti::i, ti::j));

// Compute face normal (related to the conformal metric).
auto direction = Direction<3>::upper_zeta();
auto conformal_face_normal =
unnormalized_face_normal(face_mesh, logical_to_inertial_map, direction);
const auto& face_inv_conformal_metric =
dg::project_tensor_to_boundary(inv_conformal_metric, mesh, direction);
const auto face_normal_magnitude =
magnitude(conformal_face_normal, face_inv_conformal_metric);
// Compute face normals (related to the conformal metric)
auto lower_conformal_face_normal = unnormalized_face_normal(
face_mesh, logical_to_inertial_map, Direction<3>::lower_zeta());
const auto& lower_face_inv_conformal_metric =
dg::project_tensor_to_boundary(inv_conformal_metric, mesh,
Direction<3>::lower_zeta());
const auto lower_face_normal_magnitude =
magnitude(lower_conformal_face_normal, lower_face_inv_conformal_metric);
auto upper_conformal_face_normal = unnormalized_face_normal(
face_mesh, logical_to_inertial_map, Direction<3>::upper_zeta());
const auto& upper_face_inv_conformal_metric =
dg::project_tensor_to_boundary(inv_conformal_metric, mesh,
Direction<3>::upper_zeta());
const auto upper_face_normal_magnitude =
magnitude(upper_conformal_face_normal, upper_face_inv_conformal_metric);
for (size_t d = 0; d < 3; ++d) {
conformal_face_normal.get(d) /= get(face_normal_magnitude);
lower_conformal_face_normal.get(d) /= get(lower_face_normal_magnitude);
upper_conformal_face_normal.get(d) /= get(upper_face_normal_magnitude);
}
const DirectionMap<3, tnsr::i<DataVector, 3>> conformal_face_normals(
{std::make_pair(direction, conformal_face_normal)});
{std::make_pair(Direction<3>::lower_zeta(),
lower_conformal_face_normal),
std::make_pair(Direction<3>::upper_zeta(),
upper_conformal_face_normal)});

// Compute local integrals.
Scalar<double> local_adm_mass;
Expand Down Expand Up @@ -283,6 +292,7 @@ void test_local_adm_integrals(const double& distance,

} // namespace

// [[TimeOut, 10]]
SPECTRE_TEST_CASE("Unit.PointwiseFunctions.Xcts.ObserveAdmIntegrals",
"[Unit][PointwiseFunctions]") {
// Test convergence with distance
Expand Down
Loading