Skip to content

Commit

Permalink
Merge pull request #1418 from markscheel/StrahlkorperFlatSpaceIntegral
Browse files Browse the repository at this point in the history 
Add ability to do flat-space integrals over a Strahlkorper.
  • Loading branch information
@kidder
kidder committed Mar 15, 2019
2 parents 79b75f2 + 1c3a657 commit 7642ff5
Show file tree
Hide file tree
Showing 5 changed files with 144 additions and 12 deletions.
35 changes: 35 additions & 0 deletions src/ApparentHorizons/StrahlkorperGr.cpp
View file
Expand Up @@ -528,6 +528,34 @@ Scalar<DataVector> area_element(
return area_element;
}

template <typename Frame>
Scalar<DataVector> euclidean_area_element(
const StrahlkorperTags::aliases::Jacobian<Frame>& jacobian,
const tnsr::i<DataVector, 3, Frame>& normal_one_form,
const DataVector& radius,
const tnsr::i<DataVector, 3, Frame>& r_hat) noexcept {
auto cap_theta = make_with_value<tnsr::I<DataVector, 3, Frame>>(r_hat, 0.0);
auto cap_phi = make_with_value<tnsr::I<DataVector, 3, Frame>>(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<DataVector>{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 <typename Frame>
double surface_integral_of_scalar(
const Scalar<DataVector>& area_element, const Scalar<DataVector>& scalar,
Expand Down Expand Up @@ -710,6 +738,13 @@ template Scalar<DataVector> StrahlkorperGr::area_element<Frame::Inertial>(
const DataVector& radius,
const tnsr::i<DataVector, 3, Frame::Inertial>& r_hat) noexcept;

template Scalar<DataVector>
StrahlkorperGr::euclidean_area_element<Frame::Inertial>(
const StrahlkorperTags::aliases::Jacobian<Frame::Inertial>& jacobian,
const tnsr::i<DataVector, 3, Frame::Inertial>& normal_one_form,
const DataVector& radius,
const tnsr::i<DataVector, 3, Frame::Inertial>& r_hat) noexcept;

template double StrahlkorperGr::surface_integral_of_scalar(
const Scalar<DataVector>& area_element, const Scalar<DataVector>& scalar,
const Strahlkorper<Frame::Inertial>& strahlkorper) noexcept;
Expand Down
34 changes: 33 additions & 1 deletion src/ApparentHorizons/StrahlkorperGr.hpp
View file
Expand Up @@ -150,7 +150,8 @@ Scalar<DataVector> 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`.
*/
Expand All @@ -162,6 +163,37 @@ Scalar<DataVector> area_element(
const DataVector& radius,
const tnsr::i<DataVector, 3, Frame>& 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 <typename Frame>
Scalar<DataVector> euclidean_area_element(
const StrahlkorperTags::aliases::Jacobian<Frame>& jacobian,
const tnsr::i<DataVector, 3, Frame>& normal_one_form,
const DataVector& radius,
const tnsr::i<DataVector, 3, Frame>& r_hat) noexcept;

/*!
* \ingroup SurfacesGroup
* \brief Surface integral of a scalar on a 2D `Strahlkorper`
Expand Down
29 changes: 27 additions & 2 deletions src/ApparentHorizons/Tags.hpp
View file
Expand Up @@ -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"

Expand Down Expand Up @@ -215,6 +215,30 @@ struct Tangents : db::ComputeTag {
Rhat<Frame>, Jacobian<Frame>>;
};

/// Computes the Euclidean area element on a Strahlkorper.
/// Useful for flat space integrals.
template <typename Frame>
struct EuclideanAreaElement : db::ComputeTag {
static std::string name() noexcept { return "EuclideanAreaElement"; }
static constexpr auto function =
::StrahlkorperGr::euclidean_area_element<Frame>;
using argument_tags = tmpl::list<
StrahlkorperTags::Jacobian<Frame>, StrahlkorperTags::NormalOneForm<Frame>,
StrahlkorperTags::Radius<Frame>, StrahlkorperTags::Rhat<Frame>>;
};

/// Computes the flat-space integral of a scalar over a Strahlkorper.
template <typename IntegrandTag, typename Frame>
struct EuclideanSurfaceIntegral : db::ComputeTag {
static std::string name() noexcept {
return "EuclideanSurfaceIntegral" + IntegrandTag::name();
}
static constexpr auto function =
::StrahlkorperGr::surface_integral_of_scalar<Frame>;
using argument_tags = tmpl::list<EuclideanAreaElement<Frame>, IntegrandTag,
StrahlkorperTags::Strahlkorper<Frame>>;
};

template <typename Frame>
using items_tags = tmpl::list<Strahlkorper<Frame>>;

Expand Down Expand Up @@ -254,5 +278,6 @@ struct SurfaceIntegral : db::ComputeTag {
using argument_tags = tmpl::list<AreaElement<Frame>, IntegrandTag,
StrahlkorperTags::Strahlkorper<Frame>>;
};

} // namespace Tags
} // namespace StrahlkorperGr
4 changes: 4 additions & 0 deletions src/ApparentHorizons/TagsDeclarations.hpp
View file
Expand Up @@ -31,6 +31,10 @@ template <typename Frame>
struct NormalOneForm;
template <typename Frame>
struct Tangents;
template <typename Frame>
struct EuclideanAreaElement;
template <typename IntegrandTag, typename Frame>
struct EuclideanSurfaceIntegral;
} // namespace StrahlkorperTags

namespace StrahlkorperGr {
Expand Down
54 changes: 45 additions & 9 deletions tests/Unit/ApparentHorizons/Test_StrahlkorperGr.cpp
View file
Expand Up @@ -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<Frame::Inertial>& strahlkorper) noexcept {
const auto box = db::create<
db::AddSimpleTags<StrahlkorperTags::items_tags<Frame::Inertial>>,
db::AddComputeTags<
StrahlkorperTags::compute_items_tags<Frame::Inertial>>>(strahlkorper);

// Create a Minkowski metric.
const gr::Solutions::Minkowski<3> solution{};
const double t = 0.0;
const auto& cart_coords =
db::get<StrahlkorperTags::CartesianCoords<Frame::Inertial>>(box);
const auto vars = solution.variables(
cart_coords, t, gr::Solutions::Minkowski<3>::tags<DataVector>{});
const auto& spatial_metric =
get<gr::Tags::SpatialMetric<3, Frame::Inertial, DataVector>>(vars);

const auto& normal_one_form =
db::get<StrahlkorperTags::NormalOneForm<Frame::Inertial>>(box);
const auto& r_hat = db::get<StrahlkorperTags::Rhat<Frame::Inertial>>(box);
const auto& radius = db::get<StrahlkorperTags::Radius<Frame::Inertial>>(box);
const auto& jacobian =
db::get<StrahlkorperTags::Jacobian<Frame::Inertial>>(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 <typename Solution, typename Fr>
void test_area(const Solution& solution, const Strahlkorper<Fr>& strahlkorper,
const double expected,
Expand Down Expand Up @@ -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 {
Expand Down Expand Up @@ -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",
Expand Down

0 comments on commit 7642ff5

Please sign in to comment.