From 1c3a65756f1d60fea27782bc18d5de61a96cda5f Mon Sep 17 00:00:00 2001 From: Mark Scheel Date: Thu, 28 Feb 2019 09:53:58 -0800 Subject: [PATCH] Add ability to do flat-space integrals over a Strahlkorper. --- src/ApparentHorizons/StrahlkorperGr.cpp | 35 ++++++++++++ src/ApparentHorizons/StrahlkorperGr.hpp | 34 +++++++++++- src/ApparentHorizons/Tags.hpp | 29 +++++++++- src/ApparentHorizons/TagsDeclarations.hpp | 4 ++ .../ApparentHorizons/Test_StrahlkorperGr.cpp | 54 +++++++++++++++---- 5 files changed, 144 insertions(+), 12 deletions(-) diff --git a/src/ApparentHorizons/StrahlkorperGr.cpp b/src/ApparentHorizons/StrahlkorperGr.cpp index d4eda855d329..9ca7b67edb50 100644 --- a/src/ApparentHorizons/StrahlkorperGr.cpp +++ b/src/ApparentHorizons/StrahlkorperGr.cpp @@ -527,6 +527,34 @@ Scalar area_element( return area_element; } +template +Scalar euclidean_area_element( + const StrahlkorperTags::aliases::Jacobian& jacobian, + const tnsr::i& normal_one_form, + const DataVector& radius, + const tnsr::i& r_hat) noexcept { + auto cap_theta = make_with_value>(r_hat, 0.0); + auto cap_phi = make_with_value>(r_hat, 0.0); + + for (size_t i = 0; i < 3; ++i) { + cap_theta.get(i) = jacobian.get(i, 0); + cap_phi.get(i) = jacobian.get(i, 1); + for (size_t j = 0; j < 3; ++j) { + cap_theta.get(i) += r_hat.get(i) * + (r_hat.get(j) - normal_one_form.get(j)) * + jacobian.get(j, 0); + cap_phi.get(i) += r_hat.get(i) * (r_hat.get(j) - normal_one_form.get(j)) * + jacobian.get(j, 1); + } + } + + auto area_element = Scalar{square(radius)}; + get(area_element) *= sqrt(get(dot_product(cap_theta, cap_theta)) * + get(dot_product(cap_phi, cap_phi)) - + square(get(dot_product(cap_theta, cap_phi)))); + return area_element; +} + template double surface_integral_of_scalar( const Scalar& area_element, const Scalar& scalar, @@ -709,6 +737,13 @@ template Scalar StrahlkorperGr::area_element( const DataVector& radius, const tnsr::i& r_hat) noexcept; +template Scalar +StrahlkorperGr::euclidean_area_element( + const StrahlkorperTags::aliases::Jacobian& jacobian, + const tnsr::i& normal_one_form, + const DataVector& radius, + const tnsr::i& r_hat) noexcept; + template double StrahlkorperGr::surface_integral_of_scalar( const Scalar& area_element, const Scalar& scalar, const Strahlkorper& strahlkorper) noexcept; diff --git a/src/ApparentHorizons/StrahlkorperGr.hpp b/src/ApparentHorizons/StrahlkorperGr.hpp index 29f9edfd4d4e..551d014b02af 100644 --- a/src/ApparentHorizons/StrahlkorperGr.hpp +++ b/src/ApparentHorizons/StrahlkorperGr.hpp @@ -150,7 +150,8 @@ Scalar ricci_scalar( * these input arguments depend only on the Strahlkorper, not on the * metric, and can be computed from a Strahlkorper using ComputeItems * in `StrahlkorperTags`. Note that this does not include the factor - * of \f$\sin\theta\f$, i.e., this returns \f$r^2\f$ for flat space. + * of \f$\sin\theta\f$, i.e., this returns \f$r^2\f$ for a spherical + * `Strahlkorper` in flat space. * This choice makes the area element returned here compatible with * `definite_integral` defined in `YlmSpherePack.hpp`. */ @@ -162,6 +163,37 @@ Scalar area_element( const DataVector& radius, const tnsr::i& r_hat) noexcept; +/*! + * \ingroup SurfacesGroup + * \brief Euclidean area element of a 2D `Strahlkorper`. + * + * This is useful for computing a flat-space integral over an + * arbitrarily-shaped `Strahlkorper`. + * + * \details Implements Eq. (D.13), using Eqs. (D.4) and (D.5), + * of \cite Baumgarte1996hh. Specifically, computes + * \f$\sqrt{(\Theta^i\Theta_i)(\Phi^j\Phi_j)-(\Theta^i\Phi_i)^2}\f$, + * \f$\Theta^i=\left(n^i(n_j-s_j) r J^j_\theta + r J^i_\theta\right)\f$, + * \f$\Phi^i=\left(n^i(n_j-s_j)r J^j_\phi + r J^i_\phi\right)\f$, + * and \f$\Theta^i\f$ and \f$\Phi^i\f$ are lowered by the + * Euclidean spatial metric. Here \f$J^i_\alpha\f$, \f$s_j\f$, + * \f$r\f$, and \f$n^i=n_i\f$ correspond to the input arguments + * `jacobian`, `normal_one_form`, `radius`, and `r_hat`, respectively; + * these input arguments depend only on the Strahlkorper, not on the + * metric, and can be computed from a Strahlkorper using ComputeItems + * in `StrahlkorperTags`. Note that this does not include the factor + * of \f$\sin\theta\f$, i.e., this returns \f$r^2\f$ for a spherical + * `Strahlkorper`. + * This choice makes the area element returned here compatible with + * `definite_integral` defined in `YlmSpherePack.hpp`. + */ +template +Scalar euclidean_area_element( + const StrahlkorperTags::aliases::Jacobian& jacobian, + const tnsr::i& normal_one_form, + const DataVector& radius, + const tnsr::i& r_hat) noexcept; + /*! * \ingroup SurfacesGroup * \brief Surface integral of a scalar on a 2D `Strahlkorper` diff --git a/src/ApparentHorizons/Tags.hpp b/src/ApparentHorizons/Tags.hpp index 0a9ec7c4d4c3..339e359ba65e 100644 --- a/src/ApparentHorizons/Tags.hpp +++ b/src/ApparentHorizons/Tags.hpp @@ -7,11 +7,11 @@ #include "ApparentHorizons/Strahlkorper.hpp" #include "ApparentHorizons/StrahlkorperGr.hpp" -#include "ApparentHorizons/TagsDeclarations.hpp" // IWYU pragma: keep +#include "ApparentHorizons/TagsDeclarations.hpp" // IWYU pragma: keep #include "ApparentHorizons/TagsTypeAliases.hpp" #include "DataStructures/DataBox/DataBoxTag.hpp" #include "DataStructures/DataVector.hpp" -#include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp" // IWYU pragma: keep +#include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp" // IWYU pragma: keep #include "Utilities/ForceInline.hpp" #include "Utilities/TMPL.hpp" @@ -215,6 +215,30 @@ struct Tangents : db::ComputeTag { Rhat, Jacobian>; }; +/// Computes the Euclidean area element on a Strahlkorper. +/// Useful for flat space integrals. +template +struct EuclideanAreaElement : db::ComputeTag { + static std::string name() noexcept { return "EuclideanAreaElement"; } + static constexpr auto function = + ::StrahlkorperGr::euclidean_area_element; + using argument_tags = tmpl::list< + StrahlkorperTags::Jacobian, StrahlkorperTags::NormalOneForm, + StrahlkorperTags::Radius, StrahlkorperTags::Rhat>; +}; + +/// Computes the flat-space integral of a scalar over a Strahlkorper. +template +struct EuclideanSurfaceIntegral : db::ComputeTag { + static std::string name() noexcept { + return "EuclideanSurfaceIntegral" + IntegrandTag::name(); + } + static constexpr auto function = + ::StrahlkorperGr::surface_integral_of_scalar; + using argument_tags = tmpl::list, IntegrandTag, + StrahlkorperTags::Strahlkorper>; +}; + template using items_tags = tmpl::list>; @@ -254,5 +278,6 @@ struct SurfaceIntegral : db::ComputeTag { using argument_tags = tmpl::list, IntegrandTag, StrahlkorperTags::Strahlkorper>; }; + } // namespace Tags } // namespace StrahlkorperGr diff --git a/src/ApparentHorizons/TagsDeclarations.hpp b/src/ApparentHorizons/TagsDeclarations.hpp index 4ad4f16330bf..d21795cbc6ce 100644 --- a/src/ApparentHorizons/TagsDeclarations.hpp +++ b/src/ApparentHorizons/TagsDeclarations.hpp @@ -31,6 +31,10 @@ template struct NormalOneForm; template struct Tangents; +template +struct EuclideanAreaElement; +template +struct EuclideanSurfaceIntegral; } // namespace StrahlkorperTags namespace StrahlkorperGr { diff --git a/tests/Unit/ApparentHorizons/Test_StrahlkorperGr.cpp b/tests/Unit/ApparentHorizons/Test_StrahlkorperGr.cpp index fe2470398c67..ba2c6f14320e 100644 --- a/tests/Unit/ApparentHorizons/Test_StrahlkorperGr.cpp +++ b/tests/Unit/ApparentHorizons/Test_StrahlkorperGr.cpp @@ -247,6 +247,40 @@ void test_area_element(const Solution& solution, const double& surface_radius, CHECK_ITERABLE_APPROX(get(area_element), expected(get(area_element).size())); } +void test_euclidean_area_element( + const Strahlkorper& strahlkorper) noexcept { + const auto box = db::create< + db::AddSimpleTags>, + db::AddComputeTags< + StrahlkorperTags::compute_items_tags>>(strahlkorper); + + // Create a Minkowski metric. + const gr::Solutions::Minkowski<3> solution{}; + const double t = 0.0; + const auto& cart_coords = + db::get>(box); + const auto vars = solution.variables( + cart_coords, t, gr::Solutions::Minkowski<3>::tags{}); + const auto& spatial_metric = + get>(vars); + + const auto& normal_one_form = + db::get>(box); + const auto& r_hat = db::get>(box); + const auto& radius = db::get>(box); + const auto& jacobian = + db::get>(box); + + // We are using a flat metric, so area_element and euclidean_area_element + // should be the same, and this is what we test. + const auto area_element = StrahlkorperGr::area_element( + spatial_metric, jacobian, normal_one_form, radius, r_hat); + const auto euclidean_area_element = StrahlkorperGr::euclidean_area_element( + jacobian, normal_one_form, radius, r_hat); + + CHECK_ITERABLE_APPROX(get(euclidean_area_element), get(area_element)); +} + template void test_area(const Solution& solution, const Strahlkorper& strahlkorper, const double expected, @@ -596,15 +630,15 @@ SPECTRE_TEST_CASE("Unit.ApparentHorizons.StrahlkorperGr.ExtrinsicCurvature", SPECTRE_TEST_CASE("Unit.ApparentHorizons.StrahlkorperGr.RicciScalar", "[ApparentHorizons][Unit]") { const double mass = 1.0; - test_ricci_scalar(gr::Solutions::KerrSchild(mass, {{0.0, 0.0, 0.0}}, - {{0.0, 0.0, 0.0}}), - [&mass](const auto& cartesian_coords) noexcept { - return TestHelpers::Schwarzschild::spatial_ricci( - cartesian_coords, mass); - }, - [&mass](const size_t size) noexcept { - return DataVector(size, 0.5 / square(mass)); - }); + test_ricci_scalar( + gr::Solutions::KerrSchild(mass, {{0.0, 0.0, 0.0}}, {{0.0, 0.0, 0.0}}), + [&mass](const auto& cartesian_coords) noexcept { + return TestHelpers::Schwarzschild::spatial_ricci(cartesian_coords, + mass); + }, + [&mass](const size_t size) noexcept { + return DataVector(size, 0.5 / square(mass)); + }); test_ricci_scalar( gr::Solutions::Minkowski<3>{}, [](const auto& cartesian_coords) noexcept { @@ -649,6 +683,8 @@ SPECTRE_TEST_CASE("Unit.ApparentHorizons.StrahlkorperGr.AreaElement", test_area(gr::Solutions::KerrSchild{mass, spin, center}, kerr_horizon, expected_area, expected_irreducible_mass); + + test_euclidean_area_element(kerr_horizon); } SPECTRE_TEST_CASE("Unit.ApparentHorizons.StrahlkorperGr.SurfaceInteralOfScalar",