Skip to content

Commit

Permalink
fixup: update ADM observer.
Browse files Browse the repository at this point in the history
  • Loading branch information
@iago-mendes
iago-mendes committed Sep 16, 2024
1 parent 6c65dc3 commit 0b09e1f
Show file tree
Hide file tree
Showing 3 changed files with 19 additions and 30 deletions.
23 changes: 10 additions & 13 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,15 @@ 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);

// 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,15 +125,13 @@ 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));
const auto contracted_center_of_mass_integrand = tenex::evaluate<ti::I>(
center_of_mass_integrand(ti::I, ti::J) * conformal_face_normal(ti::j));

// Take integrals
adm_mass->get() += definite_integral(
Expand All @@ -143,10 +142,8 @@ void local_adm_integrals(
definite_integral(contracted_linear_momentum_integrand.get(I) *
get(conformal_area_element),
face_mesh);
center_of_mass->get(I) +=
definite_integral(contracted_center_of_mass_integrand.get(I) *
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);
}
}
}
Expand Down
16 changes: 7 additions & 9 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 @@ -132,10 +131,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>>>;
domain::Tags::Faces<3, domain::Tags::FaceNormal<3>>>;

template <typename Metavariables, typename ArrayIndex,
typename ParallelComponent>
Expand All @@ -150,12 +149,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,
Parallel::GlobalCache<Metavariables>& cache,
const ArrayIndex& array_index, const ParallelComponent* const /*meta*/,
const ObservationValue& observation_value) const {
Expand All @@ -168,8 +166,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
10 changes: 2 additions & 8 deletions tests/Unit/Elliptic/Systems/Xcts/Events/Test_ObserveAdmIntegrals.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -263,12 +263,6 @@ void test_local_adm_integrals(const double& distance,
const DirectionMap<3, tnsr::i<DataVector, 3>> conformal_face_normals(
{std::make_pair(direction, conformal_face_normal)});

// Compute face normal vector.
const auto conformal_face_normal_vector = tenex::evaluate<ti::I>(
face_inv_conformal_metric(ti::I, ti::J) * conformal_face_normal(ti::j));
const DirectionMap<3, tnsr::I<DataVector, 3>> conformal_face_normal_vectors(
{std::make_pair(direction, conformal_face_normal_vector)});

// Compute local integrals.
Scalar<double> local_adm_mass;
tnsr::I<double, 3> local_adm_linear_momentum;
Expand All @@ -280,8 +274,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);
Expand Down

0 comments on commit 0b09e1f

Please sign in to comment.