diff --git a/src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.cpp b/src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.cpp index 5078a37bbf12..cfbaaea52421 100644 --- a/src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.cpp +++ b/src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.cpp @@ -31,12 +31,11 @@ void local_adm_integrals( const tnsr::II& inv_spatial_metric, const tnsr::ii& extrinsic_curvature, const Scalar& trace_extrinsic_curvature, + const tnsr::I& inertial_coords, const InverseJacobian& inv_jacobian, const Mesh<3>& mesh, const Element<3>& element, - const DirectionMap<3, tnsr::i>& conformal_face_normals, - const DirectionMap<3, tnsr::I>& - conformal_face_normal_vectors) { + const DirectionMap<3, tnsr::i>& conformal_face_normals) { // Initialize integrals to 0 adm_mass->get() = 0; for (int I = 0; I < 3; I++) { @@ -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 @@ -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(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( @@ -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( linear_momentum_integrand(ti::I, ti::J) * conformal_face_normal(ti::j)); - const auto contracted_center_of_mass_integrand = tenex::evaluate( - center_of_mass_integrand(ti::I, ti::J) * conformal_face_normal(ti::j)); // Take integrals adm_mass->get() += definite_integral( @@ -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); } } } diff --git a/src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.hpp b/src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.hpp index 36d6eb07843f..d0d91eeabb8a 100644 --- a/src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.hpp +++ b/src/Elliptic/Systems/Xcts/Events/ObserveAdmIntegrals.hpp @@ -56,12 +56,11 @@ void local_adm_integrals( const tnsr::II& inv_spatial_metric, const tnsr::ii& extrinsic_curvature, const Scalar& trace_extrinsic_curvature, + const tnsr::I& inertial_coords, const InverseJacobian& inv_jacobian, const Mesh<3>& mesh, const Element<3>& element, - const DirectionMap<3, tnsr::i>& conformal_face_normals, - const DirectionMap<3, tnsr::I>& - conformal_face_normal_vectors); + const DirectionMap<3, tnsr::i>& conformal_face_normals); /// @} /// @{ @@ -135,10 +134,10 @@ class ObserveAdmIntegrals : public Event { gr::Tags::InverseSpatialMetric, gr::Tags::ExtrinsicCurvature, gr::Tags::TraceExtrinsicCurvature, + 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 & inv_spatial_metric, const tnsr::ii& extrinsic_curvature, const Scalar& trace_extrinsic_curvature, + const tnsr::I& inertial_coords, const InverseJacobian& inv_jacobian, const Mesh<3>& mesh, const Element<3>& element, const DirectionMap<3, tnsr::i>& conformal_face_normals, - const DirectionMap<3, tnsr::I>& - conformal_face_normal_vectors, const ObservationBox& box, Parallel::GlobalCache& cache, const ArrayIndex& array_index, const ParallelComponent* const /*meta*/, @@ -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), diff --git a/tests/Unit/Elliptic/Systems/Xcts/Events/Test_ObserveAdmIntegrals.cpp b/tests/Unit/Elliptic/Systems/Xcts/Events/Test_ObserveAdmIntegrals.cpp index 6551333a7532..be164f695a27 100644 --- a/tests/Unit/Elliptic/Systems/Xcts/Events/Test_ObserveAdmIntegrals.cpp +++ b/tests/Unit/Elliptic/Systems/Xcts/Events/Test_ObserveAdmIntegrals.cpp @@ -263,12 +263,6 @@ void test_local_adm_integrals(const double& distance, const DirectionMap<3, tnsr::i> conformal_face_normals( {std::make_pair(direction, conformal_face_normal)}); - // Compute face normal vector. - const auto conformal_face_normal_vector = tenex::evaluate( - face_inv_conformal_metric(ti::I, ti::J) * conformal_face_normal(ti::j)); - const DirectionMap<3, tnsr::I> conformal_face_normal_vectors( - {std::make_pair(direction, conformal_face_normal_vector)}); - // Compute local integrals. Scalar local_adm_mass; tnsr::I local_adm_linear_momentum; @@ -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);