Skip to content

Commit

Permalink
Instantiate inertial->distorted for strahlkorper_in_different_frame
Browse files Browse the repository at this point in the history
  • Loading branch information
@geoffrey4444
geoffrey4444 committed Jan 24, 2024
1 parent 112d0f6 commit bfeb959
Show file tree
Hide file tree
Showing 3 changed files with 115 additions and 50 deletions.
114 changes: 85 additions & 29 deletions src/Domain/StrahlkorperTransformations.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,10 +45,11 @@ void coords_to_different_frame(
std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
functions_of_time,
const double time) {
// If we ever want a source in a different frame than Grid, we
// need to add 'if constexpr's for that case.
static_assert(std::is_same_v<SrcFrame, ::Frame::Grid>,
"Source frame must currently be Grid frame");
// If ever additional cases besides Grid->Inertial, Grid->Distorted,
// and Inertial->Distorted are needed, add if constexprs below
static_assert(std::is_same_v<SrcFrame, ::Frame::Grid> or
std::is_same_v<SrcFrame, ::Frame::Inertial>,
"Source frame must currently be Grid frame or Inertial frame");
const auto block_logical_coords = block_logical_coordinates(
domain, src_cartesian_coords, time, functions_of_time);

Expand All @@ -67,18 +68,36 @@ void coords_to_different_frame(
const auto& block_id_and_coords = block_logical_coords[s].value();
const auto& block = domain.blocks()[block_id_and_coords.id.get_index()];

if constexpr (std::is_same_v<DestFrame, ::Frame::Distorted>) {
if constexpr (std::is_same_v<DestFrame, ::Frame::Distorted> and
std::is_same_v<SrcFrame, ::Frame::Grid>) {
if (not block.has_distorted_frame()) {
throw std::runtime_error(
"Strahlkorper lies outside of distorted-frame region");
}
const auto& grid_to_distorted_map =
block.moving_mesh_grid_to_distorted_map();
x_dest = grid_to_distorted_map(x_src, time, functions_of_time);
} else if constexpr (std::is_same_v<DestFrame, ::Frame::Distorted> and
std::is_same_v<SrcFrame, ::Frame::Inertial>) {
if (not block.has_distorted_frame()) {
throw std::runtime_error(
"Strahlkorper lies outside of distorted-frame region");
}
const auto& distorted_to_inertial_map =
block.moving_mesh_distorted_to_inertial_map();
const auto& inv_point =
distorted_to_inertial_map.inverse(x_src, time, functions_of_time);
if (inv_point.has_value()) {
x_dest = *inv_point;
} else {
throw std::runtime_error(
"Map from Frame::Distorted to Frame::Inertial is not invertible");
}
} else {
static_assert(
std::is_same_v<DestFrame, ::Frame::Inertial>,
"Destination frame must currently be Distorted or Inertial");
static_assert(std::is_same_v<DestFrame, ::Frame::Inertial> and
std::is_same_v<SrcFrame, ::Frame::Grid>,
"Source frame -> destination frame must be Grid -> "
"Distorted, Grid -> Inertial, or Inertial -> Distorted");
const auto& grid_to_inertial_map =
block.moving_mesh_grid_to_inertial_map();
x_dest = grid_to_inertial_map(x_src, time, functions_of_time);
Expand Down Expand Up @@ -183,27 +202,50 @@ void strahlkorper_in_different_frame(
const auto& block_id_and_coords = block_logical_coords[0].value();
const auto& block = domain.blocks()[block_id_and_coords.id.get_index()];

static_assert(std::is_same_v<SrcFrame, ::Frame::Grid>,
"Source frame must currently be Grid frame");
static_assert(std::is_same_v<DestFrame, ::Frame::Inertial>,
"Destination frame must currently be Inertial frame");
const auto& x_src =
block.moving_mesh_logical_to_grid_map()(block_id_and_coords.data);

const double r_src =
std::hypot(get<0>(x_src) - center_src[0], get<1>(x_src) - center_src[1],
get<2>(x_src) - center_src[2]);
const double theta_src = atan2(std::hypot(get<0>(x_src) - center_src[0],
get<1>(x_src) - center_src[1]),
get<2>(x_src) - center_src[2]);
const double phi_src =
atan2(get<1>(x_src) - center_src[1], get<0>(x_src) - center_src[0]);

// Evaluate the radius of the surface at (theta_src,phi_src).
const double r_surf_src = src_strahlkorper.radius(theta_src, phi_src);

// If r_surf_src = r_src, then r_dest is on the surface.
return r_surf_src - r_src;
static_assert(
std::is_same_v<SrcFrame, ::Frame::Grid> or
std::is_same_v<SrcFrame, ::Frame::Inertial>,
"Source frame must currently be Grid frame or Inertial frame");
static_assert(std::is_same_v<DestFrame, ::Frame::Inertial> or
std::is_same_v<DestFrame, ::Frame::Distorted>,
"Destination frame must currently be Inertial frame or "
"Distorted frame");
if constexpr (std::is_same_v<SrcFrame, ::Frame::Grid>) {
const auto& x_src =
block.moving_mesh_logical_to_grid_map()(block_id_and_coords.data);
const double r_src = std::hypot(get<0>(x_src) - center_src[0],
get<1>(x_src) - center_src[1],
get<2>(x_src) - center_src[2]);
const double theta_src = atan2(std::hypot(get<0>(x_src) - center_src[0],
get<1>(x_src) - center_src[1]),
get<2>(x_src) - center_src[2]);
const double phi_src =
atan2(get<1>(x_src) - center_src[1], get<0>(x_src) - center_src[0]);

// Evaluate the radius of the surface at (theta_src,phi_src).
const double r_surf_src = src_strahlkorper.radius(theta_src, phi_src);

// If r_surf_src = r_src, then r_dest is on the surface.
return r_surf_src - r_src;
} else {
const auto& x_src = block.moving_mesh_grid_to_inertial_map()(
block.moving_mesh_logical_to_grid_map()(block_id_and_coords.data),
time, functions_of_time);
const double r_src = std::hypot(get<0>(x_src) - center_src[0],
get<1>(x_src) - center_src[1],
get<2>(x_src) - center_src[2]);
const double theta_src = atan2(std::hypot(get<0>(x_src) - center_src[0],
get<1>(x_src) - center_src[1]),
get<2>(x_src) - center_src[2]);
const double phi_src =
atan2(get<1>(x_src) - center_src[1], get<0>(x_src) - center_src[0]);

// Evaluate the radius of the surface at (theta_src,phi_src).
const double r_surf_src = src_strahlkorper.radius(theta_src, phi_src);

// If r_surf_src = r_src, then r_dest is on the surface.
return r_surf_src - r_src;
}
};

// This version of radius_function returns a double, not a
Expand Down Expand Up @@ -406,6 +448,20 @@ void strahlkorper_coords_in_different_frame(
GENERATE_INSTANTIATIONS(INSTANTIATEGENERAL, (::Frame::Grid),
(::Frame::Inertial))

// Generate a specific instantiation: inertial -> distorted
// (needed, e.g., for initializing a binary-black-hole ringdown).
// Don't just add to GENERATE_INSTANTIATIONS above, to avoid also generating
// unwanted instantiations, such as inertial -> inertial
template void strahlkorper_in_different_frame(
const gsl::not_null<ylm::Strahlkorper<::Frame::Distorted>*>
dest_strahlkorper,
const ylm::Strahlkorper<::Frame::Inertial>& src_strahlkorper,
const Domain<3>& domain,
const std::unordered_map<
std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
functions_of_time,
const double time);

#define INSTANTIATEALIGNED(_, data) \
template void strahlkorper_in_different_frame_aligned( \
const gsl::not_null<ylm::Strahlkorper<DESTFRAME(data)>*> \
Expand Down
9 changes: 5 additions & 4 deletions src/Domain/StrahlkorperTransformations.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,10 +37,11 @@ class not_null;
/// requires adding member functions to Block.
///
/// \note strahlkorper_in_different_frame is currently instantiated
/// only for SrcFrame=Grid and DestFrame=Inertial. In particular, we
/// intentionally do not instantiate strahlkorper_in_different_frame
/// for SrcFrame=Grid and DestFrame=Distorted, because our
/// Grid->Distorted maps preserve centers and angles, so for
/// only for i) SrcFrame=Grid and DestFrame=Inertial and ii) SrcFrame=Inertial
/// and DestFrame=Distorted (necessary for initializing a binary-black-hole
/// ringdown). In particular, we intentionally do not instantiate
/// strahlkorper_in_different_frame for SrcFrame=Grid and DestFrame=Distorted,
/// because our Grid->Distorted maps preserve centers and angles, so for
/// Grid->Distorted one should use
/// strahlkorper_in_different_frame_aligned because it is more
/// efficient.
Expand Down
42 changes: 25 additions & 17 deletions tests/Unit/Domain/Test_StrahlkorperTransformations.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,18 +25,18 @@

namespace {

template <bool Aligned, typename DestFrame>
template <bool Aligned, typename SrcFrame, typename DestFrame>
void test_strahlkorper_in_different_frame() {
const size_t grid_points_each_dimension = 5;

// Set up a Strahlkorper corresponding to a Schwarzschild hole of
// mass 1, in the grid frame.
// mass 1, in the source frame.
// Center the Strahlkorper at (0.03,0.02,0.01) so that we test a
// nonzero center.
const std::array<double, 3> strahlkorper_grid_center = {0.03, 0.02, 0.01};
const std::array<double, 3> strahlkorper_src_center = {0.03, 0.02, 0.01};
const size_t l_max = 8;
const ylm::Strahlkorper<Frame::Grid> strahlkorper_grid(
l_max, 2.0, strahlkorper_grid_center);
const ylm::Strahlkorper<SrcFrame> strahlkorper_src(l_max, 2.0,
strahlkorper_src_center);

// Create a Domain.
// We choose a spherical shell domain extending from radius 1.9M to
Expand All @@ -56,15 +56,16 @@ void test_strahlkorper_in_different_frame() {
std::make_unique<domain::creators::time_dependence::Shape<
domain::ObjectLabel::None>>(0.0, l_max, 1.0,
std::array<double, 3>{{0.1, 0.2, 0.3}},
strahlkorper_grid_center, 2.0, 12.0)));
strahlkorper_src_center, 2.0, 12.0)));
} else {
domain_creator.reset(new domain::creators::Sphere(
1.9, 2.9, domain::creators::Sphere::Excision{}, 1_st,
grid_points_each_dimension, false, std::nullopt, radial_partitioning,
radial_distribution, ShellWedges::All,
std::make_unique<
domain::creators::time_dependence::UniformTranslation<3>>(
0.0, std::array<double, 3>({{0.01, 0.02, 0.03}}))));
0.0, std::array<double, 3>({{0.0, 0.0, 0.0}}),
std::array<double, 3>({{0.01, 0.02, 0.03}}))));
}
Domain<3> domain = domain_creator->create_domain();
const auto functions_of_time = domain_creator->functions_of_time();
Expand All @@ -74,13 +75,13 @@ void test_strahlkorper_in_different_frame() {
ylm::Strahlkorper<DestFrame> strahlkorper_dest{};
if constexpr (Aligned) {
strahlkorper_in_different_frame_aligned(make_not_null(&strahlkorper_dest),
strahlkorper_grid, domain,
strahlkorper_src, domain,
functions_of_time, time);

} else {
strahlkorper_in_different_frame(make_not_null(&strahlkorper_dest),
strahlkorper_grid, domain,
functions_of_time, time);
strahlkorper_src, domain, functions_of_time,
time);
}

// Now compare.
Expand All @@ -92,13 +93,18 @@ void test_strahlkorper_in_different_frame() {
2.0, ylm.theta_phi_points(), 1.0,
std::array<double, 3>{{0.1, 0.2, 0.3}}));
strahlkorper_expected.reset(new ylm::Strahlkorper<DestFrame>(
l_max, l_max, new_radius, strahlkorper_grid_center));
l_max, l_max, new_radius, strahlkorper_src_center));
} else if constexpr (std::is_same_v<SrcFrame, ::Frame::Inertial> and
std::is_same_v<DestFrame, ::Frame::Distorted>) {
strahlkorper_expected.reset(new ylm::Strahlkorper<DestFrame>(

Check failure on line 99 in tests/Unit/Domain/Test_StrahlkorperTransformations.cpp

View workflow job for this annotation

GitHub Actions / Clang-tidy (Release) 

initializing non-owner argument of type 'std::unique_ptr<ylm::Strahlkorper<Frame::Distorted>>::pointer' (aka 'ylm::Strahlkorper<Frame::Distorted> *') with a newly created 'gsl::owner<>'
l_max, 2.0,
{{strahlkorper_src_center[0] - 0.005, strahlkorper_src_center[1] - 0.01,
strahlkorper_src_center[2] - 0.015}}));
} else {
strahlkorper_expected.reset(new ylm::Strahlkorper<DestFrame>(
l_max, 2.0,
{{strahlkorper_grid_center[0] + 0.005,
strahlkorper_grid_center[1] + 0.01,
strahlkorper_grid_center[2] + 0.015}}));
{{strahlkorper_src_center[0] + 0.005, strahlkorper_src_center[1] + 0.01,
strahlkorper_src_center[2] + 0.015}}));
}
CHECK_ITERABLE_APPROX(strahlkorper_expected->physical_center(),
strahlkorper_dest.physical_center());
Expand Down Expand Up @@ -190,9 +196,11 @@ SPECTRE_TEST_CASE("Unit.Domain.StrahlkorperTransformations", "[Unit]") {
domain::creators::register_derived_with_charm();
domain::creators::time_dependence::register_derived_with_charm();
domain::FunctionsOfTime::register_derived_with_charm();
test_strahlkorper_in_different_frame<false, Frame::Inertial>();
test_strahlkorper_in_different_frame<true, Frame::Inertial>();
test_strahlkorper_in_different_frame<true, Frame::Distorted>();
test_strahlkorper_in_different_frame<false, Frame::Grid, Frame::Inertial>();
test_strahlkorper_in_different_frame<false, Frame::Inertial,
Frame::Distorted>();
test_strahlkorper_in_different_frame<true, Frame::Grid, Frame::Inertial>();
test_strahlkorper_in_different_frame<true, Frame::Grid, Frame::Distorted>();
test_strahlkorper_coords_in_different_frame<true, Frame::Grid>();
test_strahlkorper_coords_in_different_frame<true, Frame::Distorted>();
test_strahlkorper_coords_in_different_frame<false, Frame::Distorted>();
Expand Down

0 comments on commit bfeb959

Please sign in to comment.