WIP use autodiff to test punctures linearization

sxs-collaboration · Nov 13, 2024 · 2565c3d · 2565c3d
1 parent 090df8b
commit 2565c3d
Show file tree

Hide file tree

Showing 2 changed files with 26 additions and 0 deletions.
diff --git 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
@@ -4,10 +4,12 @@
 #include "Framework/TestingFramework.hpp"
 
 #include <array>
+#include <autodiff/forward/real.hpp>
 #include <cstddef>
 #include <string>
 
 #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 <typename DataType, typename VarDataType>
+VarDataType puncture_sources(const Scalar<DataType>& alpha,
+                        const Scalar<DataType>& beta,
+                        const Scalar<VarDataType>& field) {
+  return get(beta) / pow<7>(get(alpha) * (get(field) + 1.) + 1.);
+}
+
+void add_linearized_sources_autodiff(
+    const gsl::not_null<Scalar<DataVector>*> linearized_puncture_equation,
+    const Scalar<DataVector>& alpha, const Scalar<DataVector>& beta,
+    const Scalar<DataVector>& field,
+    const Scalar<DataVector>& field_correction) {
+  Scalar<blaze::DynamicVector<autodiff::real>> field_dual(get(field));
+  const auto dS_du = autodiff::derivative(
+      &puncture_sources<DataVector, blaze::DynamicVector<autodiff::real>>,
+      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) {