Skip to content

Commit

Permalink
Compute ADM momenta as volume integrals
Browse files Browse the repository at this point in the history
  • Loading branch information
@iago-mendes
iago-mendes committed Jan 9, 2025
1 parent 49286a3 commit 11ba5c0
Show file tree
Hide file tree
Showing 2 changed files with 228 additions and 115 deletions.
290 changes: 186 additions & 104 deletions src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,117 +45,199 @@ 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;
// Interfaces at the inner boundary
if (boundary_direction == Direction<3>::lower_zeta()) {
// Project fields to the boundary
const auto face_conformal_factor = dg::project_tensor_to_boundary(
conformal_factor, mesh, boundary_direction);
const auto face_deriv_conformal_factor = dg::project_tensor_to_boundary(
deriv_conformal_factor, mesh, boundary_direction);
const auto face_conformal_metric = dg::project_tensor_to_boundary(
conformal_metric, mesh, boundary_direction);
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);
const auto face_conformal_christoffel_contracted =
dg::project_tensor_to_boundary(conformal_christoffel_contracted, mesh,
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_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);

// Get interface mesh and normal
const auto& face_mesh = mesh.slice_away(boundary_direction.dimension());
const auto& conformal_face_normal =
conformal_face_normals.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 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 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);
}
}

// Project fields to the boundary
const auto face_conformal_factor = dg::project_tensor_to_boundary(
conformal_factor, mesh, boundary_direction);
const auto face_deriv_conformal_factor = dg::project_tensor_to_boundary(
deriv_conformal_factor, mesh, boundary_direction);
const auto face_conformal_metric = dg::project_tensor_to_boundary(
conformal_metric, mesh, boundary_direction);
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);
const auto face_conformal_christoffel_contracted =
dg::project_tensor_to_boundary(conformal_christoffel_contracted, mesh,
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_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));

// Get interface mesh and normal
const auto& face_mesh = mesh.slice_away(boundary_direction.dimension());
const auto& conformal_face_normal =
conformal_face_normals.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 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(
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 outer boundary
if (boundary_direction == Direction<3>::upper_zeta()) {
// Project fields to the boundary
const auto face_conformal_factor = dg::project_tensor_to_boundary(
conformal_factor, mesh, boundary_direction);
const auto face_deriv_conformal_factor = dg::project_tensor_to_boundary(
deriv_conformal_factor, mesh, boundary_direction);
const auto face_conformal_metric = dg::project_tensor_to_boundary(
conformal_metric, mesh, boundary_direction);
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);
const auto face_conformal_christoffel_contracted =
dg::project_tensor_to_boundary(conformal_christoffel_contracted, mesh,
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_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);

// Get interface mesh and normal
const auto& face_mesh = mesh.slice_away(boundary_direction.dimension());
const auto& conformal_face_normal =
conformal_face_normals.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 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_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
Loading

0 comments on commit 11ba5c0

Please sign in to comment.