Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Numerically invert Jacobian in combined call #6223

Merged
merged 1 commit into from
Aug 21, 2024
Merged

Numerically invert Jacobian in combined call #6223

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 3 additions & 1 deletion src/Domain/CoordinateMaps/CoordinateMap.hpp
View file
Open in desktop
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
Open in desktop
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),
nilsvu marked this conversation as resolved.
Show resolved Hide resolved
std::move(frame_velocity)};
Comment on lines +545 to +546
Copy link
Contributor

@knelli2 knelli2 Aug 16, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Add in the CoordinateMap.hpp docs for this function that we don't call the individual inverse_jacobian functions and multiply them together, but rather just compute the jacobian and numerically invert it because most of our production time-dependent maps do this already.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a reason to keep the inverse_jacobian implementations on the individual maps at all? Doesn't need to be done in this PR, but removing them would be nicer than documenting that they are unused.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well they are still used in the individual inverse_jacobian function of CoordinateMap{Base}, correct? Whether or not that function is actually used in practice, I don't know. If not, then yes we can probably get rid of it

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Or we could make it do the numerical inversion like this function now does.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah we could! Let's discuss on the Monday call

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Another thing for a separate PR: this function could also return the Jacobian determinant. We need this a bit in the hydro cases.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

From the Monday call: We can just numerically invert the jacobian in the individual function as well which means we don't need the individual inverse_jacobian functions. See #6233

}

template <typename SourceFrame, typename TargetFrame, typename... Maps>
Expand Down
187 changes: 101 additions & 86 deletions tests/Unit/Domain/CoordinateMaps/Test_CoordinateMap.cpp
View file
Open in desktop
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 @@ -1074,14 +1085,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 @@ -1227,30 +1240,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
Loading