Skip to content

Commit

Permalink
Merge pull request #6223 from nikwit/speed-up-coord-map
Browse files Browse the repository at this point in the history 
Numerically invert Jacobian in combined call
  • Loading branch information
@nilsdeppe
nilsdeppe authored Aug 21, 2024
2 parents c5e1483 + e442112 commit c9ba57b
Show file tree
Hide file tree
Showing 3 changed files with 116 additions and 114 deletions.
4 changes: 3 additions & 1 deletion src/Domain/CoordinateMaps/CoordinateMap.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -407,7 +407,9 @@ class CoordinateMap

/// @{
/// Compute the mapped coordinates, frame velocity, Jacobian, and inverse
/// Jacobian
/// Jacobian. The inverse Jacobian is computed by numerically inverting the
/// Jacobian as this was measured to be quicker than computing it directly for
/// more complex map concatenations.
std::tuple<tnsr::I<double, dim, TargetFrame>,
InverseJacobian<double, dim, SourceFrame, TargetFrame>,
Jacobian<double, dim, SourceFrame, TargetFrame>,
Expand Down
39 changes: 12 additions & 27 deletions src/Domain/CoordinateMaps/CoordinateMap.tpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@
#include <utility>
#include <vector>

#include "DataStructures/Tensor/EagerMath/DeterminantAndInverse.hpp"
#include "DataStructures/Tensor/Identity.hpp"
#include "DataStructures/Tensor/Tensor.hpp"
#include "Domain/CoordinateMaps/CoordinateMapHelpers.hpp"
Expand Down Expand Up @@ -441,42 +442,31 @@ auto CoordinateMap<SourceFrame, TargetFrame, Maps...>::
tnsr::I<T, dim, TargetFrame>> {
check_functions_of_time(functions_of_time);
std::array<T, dim> mapped_point = make_array<T, dim>(std::move(source_point));
InverseJacobian<T, dim, SourceFrame, TargetFrame> inv_jac{};
Jacobian<T, dim, SourceFrame, TargetFrame> jac{};
tnsr::I<T, dim, TargetFrame> frame_velocity{};

tuple_transform(
maps_,
[&frame_velocity, &inv_jac, &jac, &mapped_point, time,
&functions_of_time](const auto& map, auto index,
const std::tuple<Maps...>& maps) {
[&frame_velocity, &jac, &mapped_point, time, &functions_of_time](
const auto& map, auto index, const std::tuple<Maps...>& maps) {
constexpr size_t count = decltype(index)::value;
using Map = std::decay_t<decltype(map)>;

tnsr::Ij<T, dim, Frame::NoFrame> noframe_jac{};
tnsr::Ij<T, dim, Frame::NoFrame> noframe_inv_jac{};

if (UNLIKELY(count == 0)) {
// Set Jacobian and inverse Jacobian
detail::get_inv_jacobian(
make_not_null(&noframe_inv_jac), map, mapped_point, time,
functions_of_time,
domain::is_jacobian_time_dependent_t<Map, T>{});
if constexpr (count == 0) {
// Set Jacobian
detail::get_jacobian(make_not_null(&noframe_jac), map, mapped_point,
time, functions_of_time,
domain::is_jacobian_time_dependent_t<Map, T>{});
for (size_t target = 0; target < dim; ++target) {
for (size_t source = 0; source < dim; ++source) {
jac.get(target, source) =
std::move(noframe_jac.get(target, source));
inv_jac.get(source, target) =
std::move(noframe_inv_jac.get(source, target));
}
}

// Set frame velocity
if (domain::is_map_time_dependent_v<
std::tuple_element_t<0, std::decay_t<decltype(maps)>>>) {
if constexpr (domain::is_map_time_dependent_v<std::tuple_element_t<
0, std::decay_t<decltype(maps)>>>) {
std::array<T, dim> noframe_frame_velocity =
detail::get_frame_velocity(map, mapped_point, time,
functions_of_time);
Expand All @@ -496,22 +486,16 @@ auto CoordinateMap<SourceFrame, TargetFrame, Maps...>::
// velocity is also zero. That is, we do not optimize for the map
// being instantaneously zero.

detail::get_inv_jacobian(
make_not_null(&noframe_inv_jac), map, mapped_point, time,
functions_of_time,
domain::is_jacobian_time_dependent_t<Map, T>{});
detail::get_jacobian(make_not_null(&noframe_jac), map, mapped_point,
time, functions_of_time,
domain::is_jacobian_time_dependent_t<Map, T>{});

// Perform matrix multiplication for Jacobian and inverse Jacobian
detail::multiply_inv_jacobian(make_not_null(&inv_jac),
noframe_inv_jac);
// Perform matrix multiplication for Jacobian
detail::multiply_jacobian(make_not_null(&jac), noframe_jac);

// Set frame velocity, only if map is time-dependent
std::array<T, dim> noframe_frame_velocity{};
if (domain::is_map_time_dependent_v<
if constexpr (domain::is_map_time_dependent_v<
std::tuple_element_t<count, std::decay_t<decltype(maps)>>>) {
noframe_frame_velocity = detail::get_frame_velocity(
map, mapped_point, time, functions_of_time);
Expand Down Expand Up @@ -557,8 +541,9 @@ auto CoordinateMap<SourceFrame, TargetFrame, Maps...>::
InverseJacobian<T, dim, SourceFrame, TargetFrame>,
Jacobian<T, dim, SourceFrame, TargetFrame>,
tnsr::I<T, dim, TargetFrame>>{
tnsr::I<T, dim, TargetFrame>(std::move(mapped_point)), std::move(inv_jac),
std::move(jac), std::move(frame_velocity)};
tnsr::I<T, dim, TargetFrame>(std::move(mapped_point)),
determinant_and_inverse(jac).second, std::move(jac),
std::move(frame_velocity)};
}

template <typename SourceFrame, typename TargetFrame, typename... Maps>
Expand Down
187 changes: 101 additions & 86 deletions tests/Unit/Domain/CoordinateMaps/Test_CoordinateMap.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -141,13 +141,14 @@ void test_single_coordinate_map() {

const auto coords_jacs_velocity =
map_base->coords_frame_velocity_jacobians(local_source_points);
CHECK(std::get<0>(coords_jacs_velocity) == map(local_source_points));
CHECK(std::get<1>(coords_jacs_velocity) ==
map.inv_jacobian(local_source_points));
CHECK(std::get<2>(coords_jacs_velocity) ==
map.jacobian(local_source_points));
CHECK(std::get<3>(coords_jacs_velocity) ==
tnsr::I<double, dim, Frame::Grid>{0.0});
CHECK_ITERABLE_APPROX(std::get<0>(coords_jacs_velocity),
map(local_source_points));
CHECK_ITERABLE_APPROX(std::get<1>(coords_jacs_velocity),
map.inv_jacobian(local_source_points));
CHECK_ITERABLE_APPROX(std::get<2>(coords_jacs_velocity),
map.jacobian(local_source_points));
CHECK_ITERABLE_APPROX(std::get<3>(coords_jacs_velocity),
(tnsr::I<double, dim, Frame::Grid>{0.0}));
};

for (const auto& coord : coords1d) {
Expand Down Expand Up @@ -349,11 +350,14 @@ void test_coordinate_map_with_affine_map() {

const auto coords_jacs_velocity =
map.coords_frame_velocity_jacobians(source_points);
CHECK(std::get<0>(coords_jacs_velocity) == map(source_points));
CHECK(std::get<1>(coords_jacs_velocity) == map.inv_jacobian(source_points));
CHECK(std::get<2>(coords_jacs_velocity) == map.jacobian(source_points));
CHECK(std::get<3>(coords_jacs_velocity) ==
tnsr::I<double, 1, Frame::Grid>{0.0});
CHECK_ITERABLE_APPROX(std::get<0>(coords_jacs_velocity),
map(source_points));
CHECK_ITERABLE_APPROX(std::get<1>(coords_jacs_velocity),
map.inv_jacobian(source_points));
CHECK_ITERABLE_APPROX(std::get<2>(coords_jacs_velocity),
map.jacobian(source_points));
CHECK_ITERABLE_APPROX(std::get<3>(coords_jacs_velocity),
(tnsr::I<double, 1, Frame::Grid>{0.0}));
}

// Test 2D
Expand Down Expand Up @@ -396,13 +400,14 @@ void test_coordinate_map_with_affine_map() {

const auto coords_jacs_velocity =
prod_map2d.coords_frame_velocity_jacobians(source_points);
CHECK(std::get<0>(coords_jacs_velocity) == prod_map2d(source_points));
CHECK(std::get<1>(coords_jacs_velocity) ==
prod_map2d.inv_jacobian(source_points));
CHECK(std::get<2>(coords_jacs_velocity) ==
prod_map2d.jacobian(source_points));
CHECK(std::get<3>(coords_jacs_velocity) ==
tnsr::I<double, 2, Frame::Grid>{0.0});
CHECK_ITERABLE_APPROX(std::get<0>(coords_jacs_velocity),
prod_map2d(source_points));
CHECK_ITERABLE_APPROX(std::get<1>(coords_jacs_velocity),
prod_map2d.inv_jacobian(source_points));
CHECK_ITERABLE_APPROX(std::get<2>(coords_jacs_velocity),
prod_map2d.jacobian(source_points));
CHECK_ITERABLE_APPROX(std::get<3>(coords_jacs_velocity),
(tnsr::I<double, 2, Frame::Grid>{0.0}));
}

// Test 3D
Expand Down Expand Up @@ -462,13 +467,14 @@ void test_coordinate_map_with_affine_map() {

const auto coords_jacs_velocity =
prod_map3d.coords_frame_velocity_jacobians(source_points);
CHECK(std::get<0>(coords_jacs_velocity) == prod_map3d(source_points));
CHECK(std::get<1>(coords_jacs_velocity) ==
prod_map3d.inv_jacobian(source_points));
CHECK(std::get<2>(coords_jacs_velocity) ==
prod_map3d.jacobian(source_points));
CHECK(std::get<3>(coords_jacs_velocity) ==
tnsr::I<double, 3, Frame::Grid>{0.0});
CHECK_ITERABLE_APPROX(std::get<0>(coords_jacs_velocity),
prod_map3d(source_points));
CHECK_ITERABLE_APPROX(std::get<1>(coords_jacs_velocity),
prod_map3d.inv_jacobian(source_points));
CHECK_ITERABLE_APPROX(std::get<2>(coords_jacs_velocity),
prod_map3d.jacobian(source_points));
CHECK_ITERABLE_APPROX(std::get<3>(coords_jacs_velocity),
(tnsr::I<double, 3, Frame::Grid>{0.0}));
}
}

Expand Down Expand Up @@ -516,13 +522,14 @@ void test_coordinate_map_with_rotation_map() {

const auto coords_jacs_velocity =
double_rotated2d.coords_frame_velocity_jacobians(source_points);
CHECK(std::get<0>(coords_jacs_velocity) == double_rotated2d(source_points));
CHECK(std::get<1>(coords_jacs_velocity) ==
double_rotated2d.inv_jacobian(source_points));
CHECK(std::get<2>(coords_jacs_velocity) ==
double_rotated2d.jacobian(source_points));
CHECK(std::get<3>(coords_jacs_velocity) ==
tnsr::I<double, 2, Frame::Grid>{0.0});
CHECK_ITERABLE_APPROX(std::get<0>(coords_jacs_velocity),
double_rotated2d(source_points));
CHECK_ITERABLE_APPROX(std::get<1>(coords_jacs_velocity),
double_rotated2d.inv_jacobian(source_points));
CHECK_ITERABLE_APPROX(std::get<2>(coords_jacs_velocity),
double_rotated2d.jacobian(source_points));
CHECK_ITERABLE_APPROX(std::get<3>(coords_jacs_velocity),
(tnsr::I<double, 2, Frame::Grid>{0.0}));
}

// Test 3D
Expand Down Expand Up @@ -564,13 +571,14 @@ void test_coordinate_map_with_rotation_map() {

const auto coords_jacs_velocity =
double_rotated3d.coords_frame_velocity_jacobians(source_points);
CHECK(std::get<0>(coords_jacs_velocity) == double_rotated3d(source_points));
CHECK(std::get<1>(coords_jacs_velocity) ==
double_rotated3d.inv_jacobian(source_points));
CHECK(std::get<2>(coords_jacs_velocity) ==
double_rotated3d.jacobian(source_points));
CHECK(std::get<3>(coords_jacs_velocity) ==
tnsr::I<double, 3, Frame::Grid>{0.0});
CHECK_ITERABLE_APPROX(std::get<0>(coords_jacs_velocity),
double_rotated3d(source_points));
CHECK_ITERABLE_APPROX(std::get<1>(coords_jacs_velocity),
double_rotated3d.inv_jacobian(source_points));
CHECK_ITERABLE_APPROX(std::get<2>(coords_jacs_velocity),
double_rotated3d.jacobian(source_points));
CHECK_ITERABLE_APPROX(std::get<3>(coords_jacs_velocity),
(tnsr::I<double, 3, Frame::Grid>{0.0}));
}

// Check inequivalence operator
Expand Down Expand Up @@ -614,14 +622,15 @@ void test_coordinate_map_with_rotation_map_datavector() {

const auto coords_jacs_velocity =
double_rotated2d.coords_frame_velocity_jacobians(coords2d);
CHECK(std::get<0>(coords_jacs_velocity) == double_rotated2d(coords2d));
CHECK(std::get<1>(coords_jacs_velocity) ==
double_rotated2d.inv_jacobian(coords2d));
CHECK(std::get<2>(coords_jacs_velocity) ==
double_rotated2d.jacobian(coords2d));
CHECK(std::get<3>(coords_jacs_velocity) ==
tnsr::I<DataVector, 2, Frame::Grid>{
DataVector{coords2d.get(0).size(), 0.0}});
CHECK_ITERABLE_APPROX(std::get<0>(coords_jacs_velocity),
double_rotated2d(coords2d));
CHECK_ITERABLE_APPROX(std::get<1>(coords_jacs_velocity),
double_rotated2d.inv_jacobian(coords2d));
CHECK_ITERABLE_APPROX(std::get<2>(coords_jacs_velocity),
double_rotated2d.jacobian(coords2d));
CHECK_ITERABLE_APPROX(std::get<3>(coords_jacs_velocity),
(tnsr::I<DataVector, 2, Frame::Grid>{
DataVector{coords2d.get(0).size(), 0.0}}));
}

// Test 3D
Expand Down Expand Up @@ -671,14 +680,15 @@ void test_coordinate_map_with_rotation_map_datavector() {

const auto coords_jacs_velocity =
double_rotated3d.coords_frame_velocity_jacobians(coords3d);
CHECK(std::get<0>(coords_jacs_velocity) == double_rotated3d(coords3d));
CHECK(std::get<1>(coords_jacs_velocity) ==
double_rotated3d.inv_jacobian(coords3d));
CHECK(std::get<2>(coords_jacs_velocity) ==
double_rotated3d.jacobian(coords3d));
CHECK(std::get<3>(coords_jacs_velocity) ==
tnsr::I<DataVector, 3, Frame::Grid>{
DataVector{coords3d.get(0).size(), 0.0}});
CHECK_ITERABLE_APPROX(std::get<0>(coords_jacs_velocity),
double_rotated3d(coords3d));
CHECK_ITERABLE_APPROX(std::get<1>(coords_jacs_velocity),
double_rotated3d.inv_jacobian(coords3d));
CHECK_ITERABLE_APPROX(std::get<2>(coords_jacs_velocity),
double_rotated3d.jacobian(coords3d));
CHECK_ITERABLE_APPROX(std::get<3>(coords_jacs_velocity),
(tnsr::I<DataVector, 3, Frame::Grid>{
DataVector{coords3d.get(0).size(), 0.0}}));
}
}

Expand Down Expand Up @@ -718,13 +728,14 @@ void test_coordinate_map_with_rotation_wedge() {

const auto coords_jacs_velocity =
composed_map.coords_frame_velocity_jacobians(test_point_vector);
CHECK(std::get<0>(coords_jacs_velocity) == composed_map(test_point_vector));
CHECK(std::get<1>(coords_jacs_velocity) ==
composed_map.inv_jacobian(test_point_vector));
CHECK(std::get<2>(coords_jacs_velocity) ==
composed_map.jacobian(test_point_vector));
CHECK(std::get<3>(coords_jacs_velocity) ==
tnsr::I<double, 2, Frame::Grid>{0.0});
CHECK_ITERABLE_APPROX(std::get<0>(coords_jacs_velocity),
composed_map(test_point_vector));
CHECK_ITERABLE_APPROX(std::get<1>(coords_jacs_velocity),
composed_map.inv_jacobian(test_point_vector));
CHECK_ITERABLE_APPROX(std::get<2>(coords_jacs_velocity),
composed_map.jacobian(test_point_vector));
CHECK_ITERABLE_APPROX(std::get<3>(coords_jacs_velocity),
(tnsr::I<double, 2, Frame::Grid>{0.0}));
}

void test_make_vector_coordinate_map_base() {
Expand Down Expand Up @@ -1077,14 +1088,16 @@ void test_coordinate_maps_are_identity() {
const auto coords_jacs_velocity =
wedge_composed_with_giant_identity.coords_frame_velocity_jacobians(
source_point);
CHECK(std::get<0>(coords_jacs_velocity) ==
wedge_composed_with_giant_identity(source_point));
CHECK(std::get<1>(coords_jacs_velocity) ==
wedge_composed_with_giant_identity.inv_jacobian(source_point));
CHECK(std::get<2>(coords_jacs_velocity) ==
wedge_composed_with_giant_identity.jacobian(source_point));
CHECK(std::get<3>(coords_jacs_velocity) ==
tnsr::I<double, 3, Frame::Inertial>{0.0});
CHECK_ITERABLE_APPROX(std::get<0>(coords_jacs_velocity),
wedge_composed_with_giant_identity(source_point));
CHECK_ITERABLE_APPROX(
std::get<1>(coords_jacs_velocity),
wedge_composed_with_giant_identity.inv_jacobian(source_point));
CHECK_ITERABLE_APPROX(
std::get<2>(coords_jacs_velocity),
wedge_composed_with_giant_identity.jacobian(source_point));
CHECK_ITERABLE_APPROX(std::get<3>(coords_jacs_velocity),
(tnsr::I<double, 3, Frame::Inertial>{0.0}));
}
}

Expand Down Expand Up @@ -1230,30 +1243,32 @@ void test_time_dependent_map() {
const auto velocity =
functions_of_time.at("Translation")->func_and_deriv(final_time)[1];
// The 1.5 factor comes from the Jacobian
CHECK(std::get<3>(coords_jacs_velocity) ==
tnsr::I<double, 1, Frame::Inertial>{1.5 *
velocity[velocity.size() - 1]});
CHECK_ITERABLE_APPROX(std::get<3>(coords_jacs_velocity),
(tnsr::I<double, 1, Frame::Inertial>{
1.5 * velocity[velocity.size() - 1]}));
}
{
const auto coords_jacs_velocity =
serialized_map.coords_frame_velocity_jacobians(
tnsr_datavector_logical, final_time, functions_of_time);
CHECK(
std::get<0>(coords_jacs_velocity) ==
CHECK_ITERABLE_APPROX(
std::get<0>(coords_jacs_velocity),
serialized_map(tnsr_datavector_logical, final_time, functions_of_time));
CHECK(std::get<1>(coords_jacs_velocity) ==
serialized_map.inv_jacobian(tnsr_datavector_logical, final_time,
functions_of_time));
CHECK(std::get<2>(coords_jacs_velocity) ==
serialized_map.jacobian(tnsr_datavector_logical, final_time,
functions_of_time));
CHECK_ITERABLE_APPROX(
std::get<1>(coords_jacs_velocity),
serialized_map.inv_jacobian(tnsr_datavector_logical, final_time,
functions_of_time));
CHECK_ITERABLE_APPROX(
std::get<2>(coords_jacs_velocity),
serialized_map.jacobian(tnsr_datavector_logical, final_time,
functions_of_time));
const auto velocity =
functions_of_time.at("Translation")->func_and_deriv(final_time)[1];
// The 1.5 factor comes from the Jacobian
CHECK(std::get<3>(coords_jacs_velocity) ==
tnsr::I<DataVector, 1, Frame::Inertial>{
DataVector{tnsr_datavector_logical.get(0).size(),
1.5 * velocity[velocity.size() - 1]}});
CHECK_ITERABLE_APPROX(std::get<3>(coords_jacs_velocity),
(tnsr::I<DataVector, 1, Frame::Inertial>{DataVector{
tnsr_datavector_logical.get(0).size(),
1.5 * velocity[velocity.size() - 1]}}));
}
}

Expand Down

0 comments on commit c9ba57b

Please sign in to comment.