From 41a2970531458a3b1406efd5df94899154d4aba6 Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Wed, 13 Nov 2024 11:58:04 -0800 Subject: [PATCH 1/5] Find autodiff in CMake --- CMakeLists.txt | 1 + cmake/SetupAutodiff.cmake | 48 +++++++++++++++++++++++++++++++ docs/Installation/BuildSystem.md | 2 ++ docs/Installation/Installation.md | 2 ++ 4 files changed, 53 insertions(+) create mode 100644 cmake/SetupAutodiff.cmake diff --git a/CMakeLists.txt b/CMakeLists.txt index 0d6e5143f201..4a56c2663514 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 000000000000..633fb8109da5 --- /dev/null +++ b/cmake/SetupAutodiff.cmake @@ -0,0 +1,48 @@ +# 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. + +option(SPECTRE_AUTODIFF "Enable automatic differentiation" OFF) + +if (NOT SPECTRE_AUTODIFF) + return() +endif() + +find_package(autodiff QUIET) + +if (NOT autodiff_FOUND) + if (NOT SPECTRE_FETCH_MISSING_DEPS) + message(FATAL_ERROR "Could not find autodiff. If you want to fetch " + "missing dependencies automatically, set SPECTRE_FETCH_MISSING_DEPS=ON.") + return() + endif() + message(STATUS "Fetching autodiff") + include(FetchContent) + FetchContent_Declare(autodiff + GIT_REPOSITORY https://github.com/autodiff/autodiff/ + GIT_TAG v1.1.2 + GIT_SHALLOW TRUE + ${SPECTRE_FETCHCONTENT_BASE_ARGS} + ) + set(AUTODIFF_BUILD_TESTS OFF CACHE BOOL "Build autodiff tests") + set(AUTODIFF_BUILD_EXAMPLES OFF CACHE BOOL "Build autodiff examples") + set(AUTODIFF_BUILD_PYTHON OFF CACHE BOOL "Build autodiff Python bindings") + set(AUTODIFF_BUILD_DOCS OFF CACHE BOOL "Build autodiff documentation") + FetchContent_MakeAvailable(autodiff) +endif() + +add_compile_definitions(SPECTRE_AUTODIFF) + +set_property( + GLOBAL APPEND PROPERTY SPECTRE_THIRD_PARTY_LIBS + autodiff::autodiff + ) diff --git a/docs/Installation/BuildSystem.md b/docs/Installation/BuildSystem.md index 07a005faa564..6771ac674049 100644 --- a/docs/Installation/BuildSystem.md +++ b/docs/Installation/BuildSystem.md @@ -178,6 +178,8 @@ alphabetical order): - Set to a path to a SpEC installation (the SpEC repository root) to link in SpEC libraries. In particular, the SpEC::Exporter library is linked in and enables loading SpEC data into SpECTRE. See \ref installation for details. +- SPECTRE_AUTODIFF + - Enable automatic differentation (default is `OFF`). This is experimental. - SPECTRE_DEBUG - Defines `SPECTRE_DEBUG` macro to enable `ASSERT`s and other debug checks so they can be used in Release builds. That is, you get sanity checks diff --git a/docs/Installation/Installation.md b/docs/Installation/Installation.md index e1e2044d4e75..8419d905635b 100644 --- a/docs/Installation/Installation.md +++ b/docs/Installation/Installation.md @@ -255,6 +255,8 @@ The following dependencies will be fetched automatically if you set matplotlib * [xsimd](https://github.com/xtensor-stack/xsimd) 11.0.1 or newer - for manual vectorization +* [autodiff](https://github.com/autodiff/autodiff/) 1.1.2 or newer - for + automatic differentiation * [libbacktrace](https://github.com/ianlancetaylor/libbacktrace) - to show source files and line numbers in backtraces of errors and asserts. Available by default on many systems, so you may not have to install it at all. The From 2c78d1a30eace1d4d3dd19fedd8fb309de5ad985 Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Wed, 13 Nov 2024 11:58:53 -0800 Subject: [PATCH 2/5] Test autodiff interoperability with DataStructures --- src/DataStructures/CMakeLists.txt | 4 + src/DataStructures/DynamicVector.hpp | 26 ++++++ tests/Unit/DataStructures/CMakeLists.txt | 8 ++ tests/Unit/DataStructures/Test_Autodiff.cpp | 99 +++++++++++++++++++++ 4 files changed, 137 insertions(+) create mode 100644 tests/Unit/DataStructures/Test_Autodiff.cpp diff --git a/src/DataStructures/CMakeLists.txt b/src/DataStructures/CMakeLists.txt index 8fab1a93a836..e7c6e5989401 100644 --- a/src/DataStructures/CMakeLists.txt +++ b/src/DataStructures/CMakeLists.txt @@ -83,6 +83,10 @@ target_link_libraries( Utilities ) +if (TARGET autodiff::autodiff) + target_link_libraries(${LIBRARY} PUBLIC autodiff::autodiff) +endif() + add_subdirectory(Blaze) add_subdirectory(DataBox) add_subdirectory(Python) diff --git a/src/DataStructures/DynamicVector.hpp b/src/DataStructures/DynamicVector.hpp index 24d5093a7fe1..c20d5a4a9f51 100644 --- a/src/DataStructures/DynamicVector.hpp +++ b/src/DataStructures/DynamicVector.hpp @@ -8,6 +8,10 @@ #pragma once +#ifdef SPECTRE_AUTODIFF +#include +#include +#endif // SPECTRE_AUTODIFF #include #include #include @@ -94,3 +98,25 @@ struct Options::create_from_yaml> { return result; } }; + +#ifdef SPECTRE_AUTODIFF +namespace autodiff { +namespace detail { + +template +struct VectorTraits> { + using ValueType = T; + + template + using ReplaceValueType = blaze::DynamicVector; +}; + +template +struct NumberTraits> { + using NumericType = typename NumberTraits::NumericType; + static constexpr auto Order = NumberTraits::Order; +}; + +} // namespace detail +} // namespace autodiff +#endif // SPECTRE_AUTODIFF diff --git a/tests/Unit/DataStructures/CMakeLists.txt b/tests/Unit/DataStructures/CMakeLists.txt index e5950b7adc8e..68293d1bf612 100644 --- a/tests/Unit/DataStructures/CMakeLists.txt +++ b/tests/Unit/DataStructures/CMakeLists.txt @@ -47,6 +47,10 @@ set(LIBRARY_SOURCES Test_VectorImplTestHelper.cpp ) +if (TARGET autodiff::autodiff) + list(APPEND LIBRARY_SOURCES Test_Autodiff.cpp) +endif() + add_subdirectory(DataBox) add_subdirectory(Tensor) @@ -62,6 +66,10 @@ target_link_libraries( Utilities ) +if (TARGET autodiff::autodiff) + target_link_libraries(${LIBRARY} PRIVATE autodiff::autodiff) +endif() + # [example_add_pybindings_test] spectre_add_python_bindings_test( "Unit.DataStructures.Python.DataVector" diff --git a/tests/Unit/DataStructures/Test_Autodiff.cpp b/tests/Unit/DataStructures/Test_Autodiff.cpp new file mode 100644 index 000000000000..526528f2043b --- /dev/null +++ b/tests/Unit/DataStructures/Test_Autodiff.cpp @@ -0,0 +1,99 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Framework/TestingFramework.hpp" + +#include +#include + +#include "DataStructures/DynamicVector.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "Utilities/ConstantExpressions.hpp" + +namespace { + +template +DataType f1(const DataType& x, const double param) { + return square(x) * param; +} + +template +DataType f2(const std::array& x) { + return x[0] * square(x[1]); +} +template +std::array f3(const std::array& x) { + return {{square(x[0]) * x[1], cube(x[1])}}; +} + +template +DataType f4(const tnsr::I& x, const Scalar& y) { + const auto f = tenex::evaluate(x(ti::I) * y()); + return get<0>(f); +} + +} // namespace + +SPECTRE_TEST_CASE("Unit.DataStructures.Autodiff", "[Unit][DataStructures]") { + // Test that autodiff works with the autodiff::real type in forward mode. + // The autodiff::dual type supports higher cross-derivatives in forward mode. + // Reverse mode with the autodiff::var type needs some adaptors to work with + // Blaze vectors and tensors, similar to how it is implemented for Eigen in + // the autodiff library (it works fine with single numbers though, if needed). + { + INFO("Single numbers"); + autodiff::real x = 2.0; + const auto df_dx = autodiff::derivative( + &f1, autodiff::wrt(x), autodiff::at(x, 3.0)); + CHECK(df_dx == approx(12.0)); + } + { + INFO("Vectorization"); + blaze::DynamicVector x{2.0, 3.0, 4.0}; + const auto df_dx = + autodiff::derivative(&f1>, + autodiff::wrt(x), autodiff::at(x, 3.0)); + const blaze::DynamicVector df_dx_expected{12.0, 18.0, 24.0}; + CHECK_ITERABLE_APPROX(df_dx, df_dx_expected); + } + { + INFO("Gradient"); + std::array x = {2.0, 3.0}; + const auto df_dx = autodiff::derivative( + &f2, autodiff::wrt(x[0]), autodiff::at(x)); + const auto df_dy = autodiff::derivative( + &f2, autodiff::wrt(x[1]), autodiff::at(x)); + CHECK(df_dx == approx(9.0)); + CHECK(df_dy == approx(12.0)); + // Same as above, but using the gradient convenience function + autodiff::real F{}; + std::vector grad{}; + autodiff::gradient(&f2, autodiff::wrt(x[0], x[1]), + autodiff::at(x), F, grad); + CHECK(grad[0] == approx(9.0)); + CHECK(grad[1] == approx(12.0)); + } + { + INFO("Jacobian"); + std::array x = {2.0, 3.0}; + std::array F{}; + blaze::DynamicMatrix J{}; + autodiff::jacobian(&f3, autodiff::wrt(x[0], x[1]), + autodiff::at(x), F, J); + const blaze::DynamicMatrix J_expected{{12.0, 4.0}, {0.0, 27.0}}; + CHECK_ITERABLE_APPROX(J, J_expected); + } + { + INFO("Tensors"); + tnsr::I, 2> x{}; + get<0>(x) = blaze::DynamicVector{2.0, 3.0}; + get<1>(x) = blaze::DynamicVector{4.0, 5.0}; + Scalar> y{}; + get(y) = blaze::DynamicVector{6.0, 7.0}; + const auto df_dx1 = + autodiff::derivative(&f4>, + autodiff::wrt(get<0>(x)), autodiff::at(x, y)); + const blaze::DynamicVector df_dx1_expected{6.0, 7.0}; + CHECK_ITERABLE_APPROX(df_dx1, df_dx1_expected); + } +} From 090df8bb91bd4fd92e7e4f182ed5310a40003953 Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Wed, 13 Nov 2024 13:02:09 -0800 Subject: [PATCH 3/5] Use autodiff to test map Jacobians --- src/Domain/CoordinateMaps/CMakeLists.txt | 4 ++ src/Domain/CoordinateMaps/Interval.cpp | 30 ++++++--- src/Domain/CoordinateMaps/Wedge.cpp | 61 ++++++++++++------- src/Domain/Structure/CMakeLists.txt | 4 ++ src/Domain/Structure/OrientationMap.cpp | 8 +++ .../Unit/Domain/CoordinateMaps/CMakeLists.txt | 1 + tests/Unit/Helpers/Domain/CMakeLists.txt | 4 ++ .../Domain/CoordinateMaps/TestMapHelpers.hpp | 29 +++++++++ 8 files changed, 112 insertions(+), 29 deletions(-) diff --git a/src/Domain/CoordinateMaps/CMakeLists.txt b/src/Domain/CoordinateMaps/CMakeLists.txt index 6366d317078b..b00bd98ea959 100644 --- a/src/Domain/CoordinateMaps/CMakeLists.txt +++ b/src/Domain/CoordinateMaps/CMakeLists.txt @@ -97,5 +97,9 @@ target_link_libraries( SphericalHarmonics ) +if (TARGET autodiff::autodiff) + target_link_libraries(${LIBRARY} PRIVATE autodiff::autodiff) +endif() + add_subdirectory(Python) add_subdirectory(TimeDependent) diff --git a/src/Domain/CoordinateMaps/Interval.cpp b/src/Domain/CoordinateMaps/Interval.cpp index 86a8404b2ff3..b85110fb19a2 100644 --- a/src/Domain/CoordinateMaps/Interval.cpp +++ b/src/Domain/CoordinateMaps/Interval.cpp @@ -253,18 +253,34 @@ 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)) +GENERATE_INSTANTIATIONS(INSTANTIATE_JAC, + (double, DataVector, + std::reference_wrapper, + std::reference_wrapper)) + +#ifdef SPECTRE_AUTODIFF +#include +#include +#include +GENERATE_INSTANTIATIONS(INSTANTIATE, + (autodiff::dual, autodiff::real, autodiff::var)) +#endif // SPECTRE_AUTODIFF + #undef DTYPE #undef INSTANTIATE +#undef INSTANTIATE_JAC } // namespace domain::CoordinateMaps diff --git a/src/Domain/CoordinateMaps/Wedge.cpp b/src/Domain/CoordinateMaps/Wedge.cpp index 95fedeb5a1ec..fad3a2c8b9ef 100644 --- a/src/Domain/CoordinateMaps/Wedge.cpp +++ b/src/Domain/CoordinateMaps/Wedge.cpp @@ -251,17 +251,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 +510,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 +806,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,14 +1022,16 @@ 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)) @@ -1035,9 +1039,22 @@ GENERATE_INSTANTIATIONS(INSTANTIATE_DTYPE, (2, 3), (double, DataVector, std::reference_wrapper, std::reference_wrapper)) +GENERATE_INSTANTIATIONS(INSTANTIATE_DTYPE_JAC, (2, 3), + (double, DataVector, + std::reference_wrapper, + std::reference_wrapper)) + +#ifdef SPECTRE_AUTODIFF +#include +#include +#include +GENERATE_INSTANTIATIONS(INSTANTIATE_DTYPE, (2, 3), + (autodiff::dual, autodiff::real, autodiff::var)) +#endif // SPECTRE_AUTODIFF #undef DIM #undef DTYPE #undef INSTANTIATE_DIM #undef INSTANTIATE_DTYPE +#undef INSTANTIATE_DTYPE_JAC } // namespace domain::CoordinateMaps diff --git a/src/Domain/Structure/CMakeLists.txt b/src/Domain/Structure/CMakeLists.txt index db48be9e2fef..42b1e38a19d7 100644 --- a/src/Domain/Structure/CMakeLists.txt +++ b/src/Domain/Structure/CMakeLists.txt @@ -66,3 +66,7 @@ target_link_libraries( Spectral Utilities ) + +if (TARGET autodiff::autodiff) + target_link_libraries(${LIBRARY} PRIVATE autodiff::autodiff) +endif() diff --git a/src/Domain/Structure/OrientationMap.cpp b/src/Domain/Structure/OrientationMap.cpp index 602d16f64668..435f257b6b4b 100644 --- a/src/Domain/Structure/OrientationMap.cpp +++ b/src/Domain/Structure/OrientationMap.cpp @@ -292,6 +292,14 @@ GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3), std::reference_wrapper const, std::reference_wrapper const)) +#ifdef SPECTRE_AUTODIFF +#include +#include +#include +GENERATE_INSTANTIATIONS(INSTANTIATION, (1, 2, 3), + (autodiff::dual, autodiff::real, autodiff::var)) +#endif // SPECTRE_AUTODIFF + #undef INSTANTIATION #undef DTYPE #undef DIM diff --git a/tests/Unit/Domain/CoordinateMaps/CMakeLists.txt b/tests/Unit/Domain/CoordinateMaps/CMakeLists.txt index 39c8e39eb1b6..620130b21fdc 100644 --- a/tests/Unit/Domain/CoordinateMaps/CMakeLists.txt +++ b/tests/Unit/Domain/CoordinateMaps/CMakeLists.txt @@ -41,6 +41,7 @@ target_link_libraries( DataStructures DataStructuresHelpers Domain + DomainHelpers DomainStructure ErrorHandling FunctionsOfTime diff --git a/tests/Unit/Helpers/Domain/CMakeLists.txt b/tests/Unit/Helpers/Domain/CMakeLists.txt index b3ff4bcef766..4aa9e232e3d1 100644 --- a/tests/Unit/Helpers/Domain/CMakeLists.txt +++ b/tests/Unit/Helpers/Domain/CMakeLists.txt @@ -30,6 +30,10 @@ target_link_libraries( Utilities ) +if (TARGET autodiff::autodiff) + target_link_libraries(${LIBRARY} INTERFACE autodiff::autodiff) +endif() + add_subdirectory(Amr) add_subdirectory(BoundaryConditions) add_subdirectory(Structure) diff --git a/tests/Unit/Helpers/Domain/CoordinateMaps/TestMapHelpers.hpp b/tests/Unit/Helpers/Domain/CoordinateMaps/TestMapHelpers.hpp index b79cc173fde6..71d3e3251b52 100644 --- a/tests/Unit/Helpers/Domain/CoordinateMaps/TestMapHelpers.hpp +++ b/tests/Unit/Helpers/Domain/CoordinateMaps/TestMapHelpers.hpp @@ -10,6 +10,9 @@ #include #include +#ifdef SPECTRE_AUTODIFF +#include +#endif // SPECTRE_AUTODIFF #include #include #include @@ -107,6 +110,31 @@ void test_jacobian(const Map& map, const std::array& test_point) { INFO("Test Jacobian"); CAPTURE(test_point); + constexpr size_t Dim = Map::dim; +#ifdef SPECTRE_AUTODIFF + // 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))); + } + } +#else // SPECTRE_AUTODIFF // Our default approx value is too stringent for this test Approx local_approx = Approx::custom().epsilon(1e-10).scale(1.0); const double dx = 1e-4; @@ -118,6 +146,7 @@ void test_jacobian(const Map& map, CHECK(jacobian.get(j, i) == local_approx(gsl::at(numerical_deriv_i, j))); } } +#endif // SPECTRE_AUTODIFF } template From 2565c3db3a9da8225a12d4189b2fb2ada3b81f80 Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Wed, 13 Nov 2024 13:04:29 -0800 Subject: [PATCH 4/5] WIP use autodiff to test punctures linearization --- .../Elliptic/Systems/Punctures/CMakeLists.txt | 1 + .../Systems/Punctures/Test_Equations.cpp | 25 +++++++++++++++++++ 2 files changed, 26 insertions(+) diff --git a/tests/Unit/Elliptic/Systems/Punctures/CMakeLists.txt b/tests/Unit/Elliptic/Systems/Punctures/CMakeLists.txt index 08a4b53bbe89..78b0ed769045 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 8468ef4b9a18..040475d7098d 100644 --- a/tests/Unit/Elliptic/Systems/Punctures/Test_Equations.cpp +++ b/tests/Unit/Elliptic/Systems/Punctures/Test_Equations.cpp @@ -4,10 +4,12 @@ #include "Framework/TestingFramework.hpp" #include +#include #include #include #include "DataStructures/DataVector.hpp" +#include "DataStructures/DynamicVector.hpp" #include "Elliptic/Protocols/FirstOrderSystem.hpp" #include "Elliptic/Systems/Punctures/FirstOrderSystem.hpp" #include "Elliptic/Systems/Punctures/Sources.hpp" @@ -18,6 +20,25 @@ namespace { +template +VarDataType puncture_sources(const Scalar& alpha, + const Scalar& beta, + const Scalar& field) { + return get(beta) / pow<7>(get(alpha) * (get(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) { + Scalar> field_dual(get(field)); + const auto dS_du = autodiff::derivative( + &puncture_sources>, + autodiff::wrt(get(field_dual)), autodiff::at(alpha, beta, field_dual)); + get(*linearized_puncture_equation) -= dS_du * get(field_correction); +} + void test_equations(const DataVector& used_for_size) { const double eps = 1.e-12; const auto seed = std::random_device{}(); @@ -30,6 +51,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) { From 9638b6a23efafb0824ba57bc95347091a888409f Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Wed, 13 Nov 2024 13:04:46 -0800 Subject: [PATCH 5/5] WIP make Tensor compatible --- .../Tensor/Expressions/DataTypeSupport.hpp | 5 +---- src/DataStructures/Tensor/Tensor.hpp | 22 +++++++++---------- 2 files changed, 12 insertions(+), 15 deletions(-) diff --git a/src/DataStructures/Tensor/Expressions/DataTypeSupport.hpp b/src/DataStructures/Tensor/Expressions/DataTypeSupport.hpp index e581dd59bd65..92c55422dff3 100644 --- a/src/DataStructures/Tensor/Expressions/DataTypeSupport.hpp +++ b/src/DataStructures/Tensor/Expressions/DataTypeSupport.hpp @@ -59,10 +59,7 @@ constexpr bool is_supported_number_datatype_v = /// /// \tparam X the `Tensor` data type template -struct is_supported_tensor_datatype - : std::disjunction< - std::is_same, std::is_same>, - std::is_same, std::is_same> {}; +struct is_supported_tensor_datatype : std::true_type {}; template constexpr bool is_supported_tensor_datatype_v = diff --git a/src/DataStructures/Tensor/Tensor.hpp b/src/DataStructures/Tensor/Tensor.hpp index 9cb8f8fc6968..63595706d628 100644 --- a/src/DataStructures/Tensor/Tensor.hpp +++ b/src/DataStructures/Tensor/Tensor.hpp @@ -103,17 +103,17 @@ class Tensor> { "If you are sure you need rank 5 or higher Tensor's please " "file an issue on GitHub or discuss with a core developer of " "SpECTRE."); - static_assert( - std::is_same_v> or std::is_same_v or - std::is_same_v or - std::is_same_v or - std::is_same_v or std::is_same_v or - is_spin_weighted_of_v or - is_spin_weighted_of_v or - simd::is_batch::value, - "Unsupported type. While other types are technically possible it is not " - "clear that Tensor is the correct container for them. Please seek advice " - "on the topic by discussing with the SpECTRE developers."); + // static_assert( + // std::is_same_v> or std::is_same_v or + // std::is_same_v or + // std::is_same_v or + // std::is_same_v or std::is_same_v or + // is_spin_weighted_of_v or + // is_spin_weighted_of_v or + // simd::is_batch::value, + // "Unsupported type. While other types are technically possible it is not + // " "clear that Tensor is the correct container for them. Please seek + // advice " "on the topic by discussing with the SpECTRE developers."); public: /// The type of the sequence that holds the data