Skip to content

Commit

Permalink
Merge pull request #6366 from nilsvu/py_binary_traj
Browse files Browse the repository at this point in the history 
Bind PN BinaryTrajectories in Python
  • Loading branch information
@knelli2
knelli2 authored Nov 11, 2024
2 parents e7fa9dd + 0eea3c1 commit c3c9143
Show file tree
Hide file tree
Showing 17 changed files with 319 additions and 291 deletions.
26 changes: 15 additions & 11 deletions tests/Unit/ControlSystem/Systems/Test_Expansion.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@
#include "Utilities/TMPL.hpp"

namespace Frame {
struct Distorted;
struct Grid;
struct Inertial;
} // namespace Frame

Expand All @@ -33,7 +33,7 @@ namespace {
using ExpansionMap = domain::CoordinateMaps::TimeDependent::CubicScale<3>;

using CoordMap =
domain::CoordinateMap<Frame::Distorted, Frame::Inertial, ExpansionMap>;
domain::CoordinateMap<Frame::Grid, Frame::Inertial, ExpansionMap>;

template <size_t DerivOrder>
void test_expansion_control_system() {
Expand Down Expand Up @@ -141,8 +141,14 @@ void test_expansion_control_system() {

const auto position_function = [&binary_trajectories](const double time) {
const double separation = binary_trajectories.separation(time);
return std::pair<std::array<double, 3>, std::array<double, 3>>{
{0.5 * separation, 0.0, 0.0}, {-0.5 * separation, 0.0, 0.0}};
std::array<tnsr::I<double, 3>, 2> result{};
get<0>(result[0]) = 0.5 * separation;
get<1>(result[0]) = 0.0;
get<2>(result[0]) = 0.0;
get<0>(result[1]) = -0.5 * separation;
get<1>(result[1]) = 0.0;
get<2>(result[1]) = 0.0;
return result;
};

const auto horizon_function = [&position_function, &runner,
Expand All @@ -156,16 +162,14 @@ void test_expansion_control_system() {
horizon_function);

// Grab results
std::array<double, 3> grid_position_of_a{};
std::array<double, 3> grid_position_of_b{};
std::tie(grid_position_of_a, grid_position_of_b) =
const auto grid_positions =
TestHelpers::grid_frame_horizon_centers_for_basic_control_systems<
element_component>(final_time, runner, position_function, coord_map);

// Our expected positions are just the initial positions
const std::array<double, 3> expected_grid_position_of_a{
const tnsr::I<double, 3, Frame::Grid> expected_grid_position_of_a{
{0.5 * initial_separation, 0.0, 0.0}};
const std::array<double, 3> expected_grid_position_of_b{
const tnsr::I<double, 3, Frame::Grid> expected_grid_position_of_b{
{-0.5 * initial_separation, 0.0, 0.0}};

const auto& expansion_f_of_t =
Expand All @@ -181,9 +185,9 @@ void test_expansion_control_system() {
binary_trajectories.separation(final_time) / initial_separation;
CHECK(custom_approx(expected_exp_factor) == exp_factor);

CHECK_ITERABLE_CUSTOM_APPROX(expected_grid_position_of_a, grid_position_of_a,
CHECK_ITERABLE_CUSTOM_APPROX(expected_grid_position_of_a, grid_positions[0],
custom_approx);
CHECK_ITERABLE_CUSTOM_APPROX(expected_grid_position_of_b, grid_position_of_b,
CHECK_ITERABLE_CUSTOM_APPROX(expected_grid_position_of_b, grid_positions[1],
custom_approx);
}

Expand Down
14 changes: 6 additions & 8 deletions tests/Unit/ControlSystem/Systems/Test_RotScaleTrans.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ using RotationMap = domain::CoordinateMaps::TimeDependent::Rotation<3>;
using ExpansionMap = domain::CoordinateMaps::TimeDependent::CubicScale<3>;

using CoordMap =
domain::CoordinateMap<Frame::Distorted, Frame::Inertial, ExpansionMap,
domain::CoordinateMap<Frame::Grid, Frame::Inertial, ExpansionMap,
RotationMap, TranslationMap>;

std::string create_input_string(const std::string& name) {
Expand Down Expand Up @@ -217,16 +217,14 @@ void test_rotscaletrans_control_system(const double rotation_eps = 5.0e-5) {
horizon_function);

// Grab results
std::array<double, 3> grid_position_of_a{};
std::array<double, 3> grid_position_of_b{};
std::tie(grid_position_of_a, grid_position_of_b) =
const auto grid_positions =
TestHelpers::grid_frame_horizon_centers_for_basic_control_systems<
element_component>(final_time, runner, position_function, coord_map);

// Our expected positions are just the initial positions
const std::array<double, 3> expected_grid_position_of_a{
const tnsr::I<double, 3, Frame::Grid> expected_grid_position_of_a{
{0.5 * initial_separation, 0.0, 0.0}};
const std::array<double, 3> expected_grid_position_of_b{
const tnsr::I<double, 3, Frame::Grid> expected_grid_position_of_b{
{-0.5 * initial_separation, 0.0, 0.0}};

const auto& rotation_f_of_t = dynamic_cast<
Expand Down Expand Up @@ -264,9 +262,9 @@ void test_rotscaletrans_control_system(const double rotation_eps = 5.0e-5) {
binary_trajectories.separation(final_time) / initial_separation,
expansion_factor, custom_approx1);

CHECK_ITERABLE_CUSTOM_APPROX(expected_grid_position_of_a, grid_position_of_a,
CHECK_ITERABLE_CUSTOM_APPROX(expected_grid_position_of_a, grid_positions[0],
custom_approx2);
CHECK_ITERABLE_CUSTOM_APPROX(expected_grid_position_of_b, grid_position_of_b,
CHECK_ITERABLE_CUSTOM_APPROX(expected_grid_position_of_b, grid_positions[1],
custom_approx2);
}
} // namespace
Expand Down
16 changes: 7 additions & 9 deletions tests/Unit/ControlSystem/Systems/Test_Rotation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@
#include "Utilities/TMPL.hpp"

namespace Frame {
struct Distorted;
struct Grid;
struct Inertial;
} // namespace Frame

Expand All @@ -37,7 +37,7 @@ namespace {
using RotationMap = domain::CoordinateMaps::TimeDependent::Rotation<3>;

using CoordMap =
domain::CoordinateMap<Frame::Distorted, Frame::Inertial, RotationMap>;
domain::CoordinateMap<Frame::Grid, Frame::Inertial, RotationMap>;

template <size_t DerivOrder>
void test_rotation_control_system(const bool newtonian) {
Expand Down Expand Up @@ -156,16 +156,14 @@ void test_rotation_control_system(const bool newtonian) {
horizon_function);

// Grab results
std::array<double, 3> grid_position_of_a{};
std::array<double, 3> grid_position_of_b{};
std::tie(grid_position_of_a, grid_position_of_b) =
const auto grid_positions =
TestHelpers::grid_frame_horizon_centers_for_basic_control_systems<
element_component>(final_time, runner, position_function, coord_map);

// Our expected positions are just the initial positions
const std::array<double, 3> expected_grid_position_of_a{
const tnsr::I<double, 3, Frame::Grid> expected_grid_position_of_a{
{0.5 * initial_separation, 0.0, 0.0}};
const std::array<double, 3> expected_grid_position_of_b{
const tnsr::I<double, 3, Frame::Grid> expected_grid_position_of_b{
{-0.5 * initial_separation, 0.0, 0.0}};

const auto& rotation_f_of_t = dynamic_cast<
Expand All @@ -189,9 +187,9 @@ void test_rotation_control_system(const bool newtonian) {
{0.0, 0.0, binary_trajectories.angular_velocity(final_time)}};
CHECK_ITERABLE_CUSTOM_APPROX(expected_omega, omega, custom_approx1);

CHECK_ITERABLE_CUSTOM_APPROX(expected_grid_position_of_a, grid_position_of_a,
CHECK_ITERABLE_CUSTOM_APPROX(expected_grid_position_of_a, grid_positions[0],
custom_approx2);
CHECK_ITERABLE_CUSTOM_APPROX(expected_grid_position_of_b, grid_position_of_b,
CHECK_ITERABLE_CUSTOM_APPROX(expected_grid_position_of_b, grid_positions[1],
custom_approx2);
}

Expand Down
27 changes: 15 additions & 12 deletions tests/Unit/ControlSystem/Systems/Test_Translation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@
#include "Utilities/TMPL.hpp"

namespace Frame {
struct Distorted;
struct Grid;
struct Inertial;
} // namespace Frame

Expand All @@ -39,7 +39,7 @@ namespace {
using TranslationMap = domain::CoordinateMaps::TimeDependent::Translation<3>;

using CoordMap =
domain::CoordinateMap<Frame::Distorted, Frame::Inertial, TranslationMap>;
domain::CoordinateMap<Frame::Grid, Frame::Inertial, TranslationMap>;

template <size_t DerivOrder>
void test_translation_control_system() {
Expand Down Expand Up @@ -144,9 +144,14 @@ void test_translation_control_system() {

const auto position_function = [&initial_separation,
&velocity](const double time) {
const std::array<double, 3> init_pos{{0.5 * initial_separation, 0.0, 0.0}};
return std::pair<std::array<double, 3>, std::array<double, 3>>{
init_pos + velocity * time, -init_pos + velocity * time};
std::array<tnsr::I<double, 3>, 2> result{};
get<0>(result[0]) = 0.5 * initial_separation + velocity[0] * time;
get<1>(result[0]) = velocity[1] * time;
get<2>(result[0]) = velocity[2] * time;
get<0>(result[1]) = -0.5 * initial_separation + velocity[0] * time;
get<1>(result[1]) = velocity[1] * time;
get<2>(result[1]) = velocity[2] * time;
return result;
};

const auto horizon_function = [&position_function, &runner,
Expand All @@ -160,16 +165,14 @@ void test_translation_control_system() {
horizon_function);

// Grab results
std::array<double, 3> grid_position_of_a{};
std::array<double, 3> grid_position_of_b{};
std::tie(grid_position_of_a, grid_position_of_b) =
const auto grid_positions =
TestHelpers::grid_frame_horizon_centers_for_basic_control_systems<
element_component>(final_time, runner, position_function, coord_map);

// Our expected positions are just the initial positions
const std::array<double, 3> expected_grid_position_of_a{
const tnsr::I<double, 3, Frame::Grid> expected_grid_position_of_a{
{0.5 * initial_separation, 0.0, 0.0}};
const std::array<double, 3> expected_grid_position_of_b{
const tnsr::I<double, 3, Frame::Grid> expected_grid_position_of_b{
{-0.5 * initial_separation, 0.0, 0.0}};

const auto& translation_f_of_t =
Expand All @@ -190,9 +193,9 @@ void test_translation_control_system() {
CHECK_ITERABLE_CUSTOM_APPROX(
trans_and_2_derivs[0], DataVector(velocity * final_time), custom_approx);

CHECK_ITERABLE_CUSTOM_APPROX(expected_grid_position_of_a, grid_position_of_a,
CHECK_ITERABLE_CUSTOM_APPROX(expected_grid_position_of_a, grid_positions[0],
custom_approx);
CHECK_ITERABLE_CUSTOM_APPROX(expected_grid_position_of_b, grid_position_of_b,
CHECK_ITERABLE_CUSTOM_APPROX(expected_grid_position_of_b, grid_positions[1],
custom_approx);
}

Expand Down
46 changes: 14 additions & 32 deletions tests/Unit/Helpers/ControlSystem/SystemHelpers.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -417,42 +417,16 @@ double initialize_translation_functions_of_time(

template <typename ElementComponent, typename Metavars, typename F,
typename CoordMap>
std::pair<std::array<double, 3>, std::array<double, 3>>
std::array<tnsr::I<double, 3, Frame::Grid>, 2>
grid_frame_horizon_centers_for_basic_control_systems(
const double time, ActionTesting::MockRuntimeSystem<Metavars>& runner,
const F position_function, const CoordMap& coord_map) {
auto& cache = ActionTesting::cache<ElementComponent>(runner, 0);
const auto& functions_of_time =
Parallel::get<domain::Tags::FunctionsOfTime>(cache);

// This whole switching between tensors and arrays is annoying and
// clunky, but it's the best that could be done at the moment without
// changing BinaryTrajectories to return tensors, which doesn't seem like
// a good idea.

std::pair<std::array<double, 3>, std::array<double, 3>> positions =
position_function(time);

// Covert arrays to tensor so we can pass them into the coordinate map
const tnsr::I<double, 3, Frame::Inertial> inertial_position_of_a(
positions.first);
const tnsr::I<double, 3, Frame::Inertial> inertial_position_of_b(
positions.second);

// Convert to "grid coordinates"
const auto grid_position_of_a_tnsr =
*coord_map.inverse(inertial_position_of_a, time, functions_of_time);
const auto grid_position_of_b_tnsr =
*coord_map.inverse(inertial_position_of_b, time, functions_of_time);

// Convert tensors back to arrays so we can pass them to the control
// systems. Just reuse `positions`
for (size_t i = 0; i < 3; i++) {
gsl::at(positions.first, i) = grid_position_of_a_tnsr.get(i);
gsl::at(positions.second, i) = grid_position_of_b_tnsr.get(i);
}

return positions;
const auto positions = position_function(time);
return {{*coord_map.inverse(positions[0], time, functions_of_time),
*coord_map.inverse(positions[1], time, functions_of_time)}};
}

template <typename ElementComponent, typename Metavars, typename F,
Expand All @@ -468,8 +442,16 @@ build_horizons_for_basic_control_systems(

// Construct strahlkorpers to pass to control systems. Only the centers
// matter.
ylm::Strahlkorper<Frame::Distorted> horizon_a{2, 2, 1.0, positions.first};
ylm::Strahlkorper<Frame::Distorted> horizon_b{2, 2, 1.0, positions.second};
ylm::Strahlkorper<Frame::Distorted> horizon_a{
2,
2,
1.0,
{{get<0>(positions[0]), get<1>(positions[0]), get<2>(positions[0])}}};
ylm::Strahlkorper<Frame::Distorted> horizon_b{
2,
2,
1.0,
{{get<0>(positions[1]), get<1>(positions[1]), get<2>(positions[1])}}};

return std::make_pair<ylm::Strahlkorper<Frame::Distorted>,
ylm::Strahlkorper<Frame::Distorted>>(
Expand Down
Loading

0 comments on commit c3c9143

Please sign in to comment.