diff --git a/CMakeLists.txt b/CMakeLists.txt index 0d6e5143f201a..4a56c2663514f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -154,6 +154,7 @@ include(SetupLapack) include(SetupLIBXSMM) include(SetupGsl) +include(SetupAutodiff) include(SetupBlaze) include(SetupFuka) include(SetupGoogleBenchmark) diff --git a/cmake/SetupAutodiff.cmake b/cmake/SetupAutodiff.cmake new file mode 100644 index 0000000000000..d82b37b610430 --- /dev/null +++ b/cmake/SetupAutodiff.cmake @@ -0,0 +1,14 @@ +# Distributed under the MIT License. +# See LICENSE.txt for details. + +# Find autodiff (https://github.com/autodiff/autodiff) +# +# Note that Boost also has an autodiff module, but it uses templates to +# distinguish multiple variables. This means a std::array or Tensor can't hold +# the Boost autodiff variables. For example, to evaluate a 2D function f(x, y), +# the type for x would be fvar and the type for y would be the nested +# fvar>. On the other hand, the autodiff library works with a +# simpler autodiff::dual or autodiff::var type, which can be stored in a +# std::array or Tensor. + +find_package(autodiff REQUIRED) diff --git a/src/Domain/CoordinateMaps/CMakeLists.txt b/src/Domain/CoordinateMaps/CMakeLists.txt index 6366d317078be..251088fb5754c 100644 --- a/src/Domain/CoordinateMaps/CMakeLists.txt +++ b/src/Domain/CoordinateMaps/CMakeLists.txt @@ -83,6 +83,7 @@ spectre_target_headers( target_link_libraries( ${LIBRARY} PRIVATE + autodiff::autodiff RootFinding PUBLIC Boost::boost diff --git a/src/Domain/CoordinateMaps/Interval.cpp b/src/Domain/CoordinateMaps/Interval.cpp index 86a8404b2ff39..240ee8c9bc407 100644 --- a/src/Domain/CoordinateMaps/Interval.cpp +++ b/src/Domain/CoordinateMaps/Interval.cpp @@ -3,6 +3,8 @@ #include "Domain/CoordinateMaps/Interval.hpp" +#include +#include #include #include #include @@ -253,18 +255,25 @@ bool operator==(const CoordinateMaps::Interval& lhs, // Explicit instantiations #define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data) -#define INSTANTIATE(_, data) \ - template std::array, 1> \ - Interval::operator()(const std::array& source_coords) const; \ - template tnsr::Ij, 1, Frame::NoFrame> \ - Interval::jacobian(const std::array& source_coords) const; \ - template tnsr::Ij, 1, Frame::NoFrame> \ - Interval::inv_jacobian(const std::array& source_coords) \ +#define INSTANTIATE(_, data) \ + template std::array, 1> \ + Interval::operator()(const std::array& source_coords) const; + +#define INSTANTIATE_JAC(_, data) \ + template tnsr::Ij, 1, Frame::NoFrame> \ + Interval::jacobian(const std::array& source_coords) const; \ + template tnsr::Ij, 1, Frame::NoFrame> \ + Interval::inv_jacobian(const std::array& source_coords) \ const; GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector, std::reference_wrapper, - std::reference_wrapper)) + std::reference_wrapper, + autodiff::dual, autodiff::var)) +GENERATE_INSTANTIATIONS(INSTANTIATE_JAC, + (double, DataVector, + std::reference_wrapper, + std::reference_wrapper)) #undef DTYPE #undef INSTANTIATE } // namespace domain::CoordinateMaps diff --git a/src/Domain/CoordinateMaps/Wedge.cpp b/src/Domain/CoordinateMaps/Wedge.cpp index 95fedeb5a1ecf..8d114deac09b6 100644 --- a/src/Domain/CoordinateMaps/Wedge.cpp +++ b/src/Domain/CoordinateMaps/Wedge.cpp @@ -3,6 +3,8 @@ #include "Domain/CoordinateMaps/Wedge.hpp" +#include +#include #include #include #include @@ -251,17 +253,17 @@ template tt::remove_cvref_wrap_t Wedge::get_cap_angular_function( const T& lowercase_xi_or_eta) const { constexpr auto cap_index = static_cast(not FuncIsXi); + if (not with_equiangular_map_) { + return lowercase_xi_or_eta; + } if (opening_angles_.has_value() and opening_angles_distribution_.has_value()) { - return with_equiangular_map_ - ? tan(0.5 * opening_angles_.value()[cap_index]) * - tan(0.5 * opening_angles_distribution_.value()[cap_index] * - lowercase_xi_or_eta) / - tan(0.5 * opening_angles_distribution_.value()[cap_index]) - : lowercase_xi_or_eta; + return tan(0.5 * opening_angles_.value()[cap_index]) * + tan(0.5 * opening_angles_distribution_.value()[cap_index] * + lowercase_xi_or_eta) / + tan(0.5 * opening_angles_distribution_.value()[cap_index]); } else { - return with_equiangular_map_ ? tan(M_PI_4 * lowercase_xi_or_eta) - : lowercase_xi_or_eta; + return tan(M_PI_4 * lowercase_xi_or_eta); } } @@ -510,11 +512,14 @@ std::array, Dim> Wedge::operator()( const bool zero_offset = not cube_half_length_.has_value(); std::array physical_coords{}; - physical_coords[radial_coord] = - zero_offset ? generalized_z - : generalized_z * (1.0 - rotated_focus[radial_coord] / - cube_half_length_.value()) + - rotated_focus[radial_coord]; + if (zero_offset) { + physical_coords[radial_coord] = generalized_z; + } else { + physical_coords[radial_coord] = + generalized_z * + (1.0 - rotated_focus[radial_coord] / cube_half_length_.value()) + + rotated_focus[radial_coord]; + } if (zero_offset) { physical_coords[polar_coord] = generalized_z * cap[0]; } else { @@ -803,7 +808,6 @@ Wedge::inv_jacobian(const std::array& source_coords) const { xi *= 0.5; } - std::array cap{}; std::array cap_deriv{}; cap[0] = get_cap_angular_function(xi); @@ -1020,18 +1024,25 @@ bool operator!=(const Wedge& lhs, const Wedge& rhs) { #define INSTANTIATE_DTYPE(_, data) \ template std::array, DIM(data)> \ Wedge::operator()( \ - const std::array& source_coords) const; \ - template tnsr::Ij, DIM(data), \ - Frame::NoFrame> \ - Wedge::jacobian( \ - const std::array& source_coords) const; \ - template tnsr::Ij, DIM(data), \ - Frame::NoFrame> \ - Wedge::inv_jacobian( \ + const std::array& source_coords) const; + +#define INSTANTIATE_DTYPE_JAC(_, data) \ + template tnsr::Ij, DIM(data), \ + Frame::NoFrame> \ + Wedge::jacobian( \ + const std::array& source_coords) const; \ + template tnsr::Ij, DIM(data), \ + Frame::NoFrame> \ + Wedge::inv_jacobian( \ const std::array& source_coords) const; GENERATE_INSTANTIATIONS(INSTANTIATE_DIM, (2, 3)) GENERATE_INSTANTIATIONS(INSTANTIATE_DTYPE, (2, 3), + (double, DataVector, + std::reference_wrapper, + std::reference_wrapper, + autodiff::dual, autodiff::var)) +GENERATE_INSTANTIATIONS(INSTANTIATE_DTYPE_JAC, (2, 3), (double, DataVector, std::reference_wrapper, std::reference_wrapper)) diff --git a/src/Domain/Structure/OrientationMap.cpp b/src/Domain/Structure/OrientationMap.cpp index 602d16f64668b..47dd032ab082d 100644 --- a/src/Domain/Structure/OrientationMap.cpp +++ b/src/Domain/Structure/OrientationMap.cpp @@ -3,6 +3,8 @@ #include "Domain/Structure/OrientationMap.hpp" +#include +#include #include #include #include @@ -261,12 +263,12 @@ tnsr::Ij discrete_rotation_inverse_jacobian( #define DIM(data) BOOST_PP_TUPLE_ELEM(0, data) -#define INSTANTIATION(r, data) \ - template class OrientationMap; \ - template tnsr::Ij \ - discrete_rotation_jacobian(const OrientationMap& orientation); \ - template tnsr::Ij \ - discrete_rotation_inverse_jacobian( \ +#define INSTANTIATION(r, data) \ + template class OrientationMap; \ + template tnsr::Ij \ + discrete_rotation_jacobian(const OrientationMap& orientation); \ + template tnsr::Ij \ + discrete_rotation_inverse_jacobian( \ const OrientationMap& orientation); GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3)) @@ -290,7 +292,8 @@ GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3), std::reference_wrapper const, std::reference_wrapper const, std::reference_wrapper const, - std::reference_wrapper const)) + std::reference_wrapper const, + autodiff::dual, autodiff::var)) #undef INSTANTIATION #undef DTYPE diff --git a/src/Utilities/CMakeLists.txt b/src/Utilities/CMakeLists.txt index 6b1cbccbff5d3..d0c1b157d9def 100644 --- a/src/Utilities/CMakeLists.txt +++ b/src/Utilities/CMakeLists.txt @@ -90,6 +90,7 @@ spectre_target_headers( target_link_libraries( ${LIBRARY} PUBLIC + autodiff::autodiff BLAS::BLAS Blaze Boost::boost diff --git a/tests/Unit/Domain/CoordinateMaps/CMakeLists.txt b/tests/Unit/Domain/CoordinateMaps/CMakeLists.txt index 39c8e39eb1b69..73264edc9d449 100644 --- a/tests/Unit/Domain/CoordinateMaps/CMakeLists.txt +++ b/tests/Unit/Domain/CoordinateMaps/CMakeLists.txt @@ -4,32 +4,8 @@ set(LIBRARY "Test_CoordinateMaps") set(LIBRARY_SOURCES - Test_Affine.cpp - Test_BulgedCube.cpp - Test_Composition.cpp - Test_CoordinateMap.cpp - Test_CylindricalEndcap.cpp - Test_CylindricalFlatEndcap.cpp - Test_CylindricalFlatSide.cpp - Test_CylindricalSide.cpp - Test_DiscreteRotation.cpp - Test_Distribution.cpp - Test_EquatorialCompression.cpp - Test_Equiangular.cpp - Test_Frustum.cpp - Test_Identity.cpp Test_Interval.cpp - Test_KerrHorizonConforming.cpp - Test_ProductMaps.cpp - Test_Rotation.cpp - Test_SpecialMobius.cpp - Test_Tags.cpp - Test_TimeDependentHelpers.cpp - Test_UniformCylindricalEndcap.cpp - Test_UniformCylindricalFlatEndcap.cpp - Test_UniformCylindricalSide.cpp Test_Wedge2D.cpp - Test_Wedge3D.cpp ) add_test_library(${LIBRARY} "${LIBRARY_SOURCES}") @@ -49,5 +25,3 @@ target_link_libraries( Spectral Utilities ) - -add_subdirectory(TimeDependent) diff --git a/tests/Unit/Domain/CoordinateMaps/Test_Interval.cpp b/tests/Unit/Domain/CoordinateMaps/Test_Interval.cpp index 42f9fd466a061..d0abb7f131619 100644 --- a/tests/Unit/Domain/CoordinateMaps/Test_Interval.cpp +++ b/tests/Unit/Domain/CoordinateMaps/Test_Interval.cpp @@ -30,6 +30,7 @@ void test_coordinate_map( test_inverse_map(map, random_point); test_jacobian(map, random_point); test_inv_jacobian(map, random_point); + test_jacobian_autodiff(map, random_point); } const std::array point_xA{{xA}}; const std::array point_xB{{xB}}; diff --git a/tests/Unit/Elliptic/Systems/Punctures/CMakeLists.txt b/tests/Unit/Elliptic/Systems/Punctures/CMakeLists.txt index 08a4b53bbe891..78b0ed769045a 100644 --- a/tests/Unit/Elliptic/Systems/Punctures/CMakeLists.txt +++ b/tests/Unit/Elliptic/Systems/Punctures/CMakeLists.txt @@ -13,6 +13,7 @@ add_test_library(${LIBRARY} "${LIBRARY_SOURCES}") target_link_libraries( ${LIBRARY} PRIVATE + autodiff::autodiff DataStructuresHelpers Punctures ) diff --git a/tests/Unit/Elliptic/Systems/Punctures/Test_Equations.cpp b/tests/Unit/Elliptic/Systems/Punctures/Test_Equations.cpp index 8468ef4b9a182..4aba0a569e4e5 100644 --- a/tests/Unit/Elliptic/Systems/Punctures/Test_Equations.cpp +++ b/tests/Unit/Elliptic/Systems/Punctures/Test_Equations.cpp @@ -4,6 +4,7 @@ #include "Framework/TestingFramework.hpp" #include +#include #include #include @@ -18,6 +19,35 @@ namespace { +template +void add_sources(const gsl::not_null puncture_equation, + const DataType& alpha, const DataType& beta, + const VarDataType& field) { + *puncture_equation -= beta / pow<7>(alpha * (field + 1.) + 1.); +} + +void add_linearized_sources_autodiff( + const gsl::not_null*> linearized_puncture_equation, + const Scalar& alpha, const Scalar& beta, + const Scalar& field, + const Scalar& field_correction) { + // Works: + // - Going over each element of the DataVector, evaluating autodiff pointwise + // TODO: + // - Vectorize the autodiff evaluation, e.g. by using a new type + // DualDataVector = VectorImpl (or the reverse mode + // equivalent) + // - Add support for Tensor (probably just extend the + // static_assert to allow this). + for (size_t i = 0; i < get(field).size(); ++i) { + const autodiff::var field_i{get(field)[i]}; + autodiff::var result{get(*linearized_puncture_equation)[i]}; + add_sources(make_not_null(&result), get(alpha)[i], get(beta)[i], field_i); + const auto dS_du = autodiff::derivatives(result, autodiff::wrt(field_i))[0]; + get(*linearized_puncture_equation)[i] += dS_du * get(field_correction)[i]; + } +} + void test_equations(const DataVector& used_for_size) { const double eps = 1.e-12; const auto seed = std::random_device{}(); @@ -30,6 +60,10 @@ void test_equations(const DataVector& used_for_size) { &Punctures::add_linearized_sources, "Equations", {"linearized_sources"}, {{{-1., 1.}, {-1., 1.}, {-1., 1.}, {-1., 1.}}}, used_for_size, eps, seed, fill_result_tensors); + pypp::check_with_random_values<4>( + &add_linearized_sources_autodiff, "Equations", {"linearized_sources"}, + {{{-1., 1.}, {-1., 1.}, {-1., 1.}, {-1., 1.}}}, used_for_size, eps, seed, + fill_result_tensors); } void test_computers(const DataVector& used_for_size) { diff --git a/tests/Unit/Helpers/Domain/CMakeLists.txt b/tests/Unit/Helpers/Domain/CMakeLists.txt index b3ff4bcef7666..d52851f057a07 100644 --- a/tests/Unit/Helpers/Domain/CMakeLists.txt +++ b/tests/Unit/Helpers/Domain/CMakeLists.txt @@ -12,6 +12,7 @@ add_spectre_library(${LIBRARY} ${SPECTRE_TEST_LIBS_TYPE} ${LIBRARY_SOURCES}) target_link_libraries( ${LIBRARY} INTERFACE + autodiff::autodiff Boost::boost CoordinateMaps FunctionsOfTime diff --git a/tests/Unit/Helpers/Domain/CoordinateMaps/TestMapHelpers.hpp b/tests/Unit/Helpers/Domain/CoordinateMaps/TestMapHelpers.hpp index b79cc173fde60..0485445098e49 100644 --- a/tests/Unit/Helpers/Domain/CoordinateMaps/TestMapHelpers.hpp +++ b/tests/Unit/Helpers/Domain/CoordinateMaps/TestMapHelpers.hpp @@ -10,6 +10,7 @@ #include #include +#include #include #include #include @@ -120,6 +121,36 @@ void test_jacobian(const Map& map, } } +template +void test_jacobian_autodiff(const Map& map, + const std::array& test_point) { + INFO("Test Jacobian with autodiff"); + CAPTURE(test_point); + constexpr size_t Dim = Map::dim; + // Use reverse mode because we want the full Jacobian (all derivatives) + std::array x{}; + for (size_t i = 0; i < Dim; ++i) { + gsl::at(x, i) = gsl::at(test_point, i); + } + const auto y = map(x); + const auto jacobian = map.jacobian(test_point); + for (size_t i = 0; i < Dim; ++i) { + std::array deriv_i{}; + if constexpr (Dim == 1) { + deriv_i = autodiff::derivatives(gsl::at(y, i), autodiff::wrt(x[0])); + } else if constexpr (Dim == 2) { + deriv_i = autodiff::derivatives(gsl::at(y, i), autodiff::wrt(x[0], x[1])); + } else if constexpr (Dim == 3) { + deriv_i = + autodiff::derivatives(gsl::at(y, i), autodiff::wrt(x[0], x[1], x[2])); + } + for (size_t j = 0; j < Dim; ++j) { + INFO("i: " << i << " j: " << j); + CHECK(jacobian.get(i, j) == approx(gsl::at(deriv_i, j))); + } + } +} + template void test_jacobian( const Map& map, const std::array& test_point, @@ -586,16 +617,19 @@ void test_suite_for_map_on_unit_cube(const Map& map) { test_coordinate_map_argument_types(map_to_test, origin); test_jacobian(map_to_test, origin); + test_jacobian_autodiff(map_to_test, origin); test_inv_jacobian(map_to_test, origin); test_inverse_map(map_to_test, origin); for (VolumeCornerIterator vci{}; vci; ++vci) { test_jacobian(map_to_test, vci.coords_of_corner()); + test_jacobian_autodiff(map_to_test, vci.coords_of_corner()); test_inv_jacobian(map_to_test, vci.coords_of_corner()); test_inverse_map(map_to_test, vci.coords_of_corner()); } test_jacobian(map_to_test, random_point); + test_jacobian_autodiff(map_to_test, random_point); test_inv_jacobian(map_to_test, random_point); test_inverse_map(map_to_test, random_point); };