Compute acceleration terms

src/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/IteratePunctureField.hpp

@@ -59,15 +59,31 @@ struct IteratePunctureField {
             const auto iterated_puncture_field) {
           tnsr::I<double, Dim> iterated_acceleration{
               {self_force_data[0], self_force_data[1], self_force_data[2]}};
+          const size_t face_size =
+              get<0>(centered_face_coordinates.value()).size();
+
           if (not iterated_puncture_field->has_value()) {
-            iterated_puncture_field->emplace(
-                get<0>(centered_face_coordinates.value()).size());
+            iterated_puncture_field->emplace(face_size);
           }
 
           puncture_field(make_not_null(&iterated_puncture_field->value()),
                          centered_face_coordinates.value(),
                          position_velocity[0], position_velocity[1],
                          iterated_acceleration, 1., order);
+          Variables<tmpl::list<CurvedScalarWave::Tags::Psi,
+                               ::Tags::dt<CurvedScalarWave::Tags::Psi>,
+                               ::Tags::deriv<CurvedScalarWave::Tags::Psi,
+                                             tmpl::size_t<3>, Frame::Inertial>>>
+              acceleration_terms(face_size);
+          acceleration_terms_1(
+              make_not_null(&acceleration_terms),
+              centered_face_coordinates.value(), position_velocity[0],
+              position_velocity[1], iterated_acceleration, self_force_data[3],
+              self_force_data[4], self_force_data[5], self_force_data[6],
+              self_force_data[7], self_force_data[8], self_force_data[9],
+              self_force_data[10], self_force_data[11], self_force_data[12],
+              self_force_data[13], self_force_data[14], 1.);
+          iterated_puncture_field->value() += acceleration_terms;
           iterated_puncture_field->value() *= charge;
         },
         make_not_null(&box));

src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/IterateAccelerationTerms.cpp

@@ -8,6 +8,7 @@
 #include <tuple>
 
 #include "Evolution/Systems/CurvedScalarWave/Worldtube/Inboxes.hpp"
+#include "Evolution/Systems/CurvedScalarWave/Worldtube/KerrSchildDerivatives.hpp"
 #include "Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.hpp"
 #include "Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp"
 #include "NumericalAlgorithms/SphericalHarmonics/Tags.hpp"
@@ -25,20 +26,183 @@ void IterateAccelerationTerms::apply(
         Tags::TimeDilationFactor>& background,
     const tnsr::I<double, Dim, Frame::Inertial>& geodesic_acc,
     const Scalar<double>& psi_monopole, const Scalar<double>& dt_psi_monopole,
-    const tnsr::i<double, Dim, Frame::Inertial>& psi_dipole, double charge,
-    std::optional<double> mass, double time, std::optional<double> turn_on_time,
-    std::optional<double> turn_on_interval) {
-  const size_t data_size = 3;
+    const tnsr::i<double, Dim, Frame::Inertial>& psi_dipole,
+    const tnsr::i<double, Dim, Frame::Inertial>& dt_psi_dipole,
+    const double charge, const std::optional<double> mass, const double time,
+    const std::optional<double> turn_on_time,
+    const std::optional<double> turn_on_interval, const size_t iteration) {
+  const size_t data_size = 15;
   acceleration_terms->get() = DataVector(data_size, 0.);
-  const auto& inverse_metric =
-      get<gr::Tags::InverseSpacetimeMetric<double, Dim>>(background);
-  const auto& dilation_factor = get<Tags::TimeDilationFactor>(background);
-  const auto& vel = pos_vel.at(1);
-  const auto self_force_acc =
-      self_force_acceleration(dt_psi_monopole, psi_dipole, vel, charge,
-                              mass.value(), inverse_metric, dilation_factor);
+  auto acc = geodesic_acc;
+  double roll_on = 0.;
+  if (time > turn_on_time.value()) {
+    const auto& particle_position = pos_vel.at(0);
+    const auto& particle_velocity = pos_vel.at(1);
+    const auto& imetric =
+        get<gr::Tags::InverseSpacetimeMetric<double, Dim>>(background);
+    const auto& metric =
+        get<gr::Tags::SpacetimeMetric<double, Dim>>(background);
+    const auto& christoffel =
+        get<gr::Tags::SpacetimeChristoffelSecondKind<double, Dim>>(background);
+    const auto& trace_christoffel =
+        get<gr::Tags::TraceSpacetimeChristoffelSecondKind<double, Dim>>(
+            background);
+    const auto di_imetric =
+        spatial_derivative_inverse_ks_metric(particle_position);
+    const auto dij_imetric =
+        second_spatial_derivative_inverse_ks_metric(particle_position);
+    const auto di_trace_christoffel =
+        spatial_derivative_ks_contracted_christoffel(particle_position);
+    const auto di_metric = spatial_derivative_ks_metric(metric, di_imetric);
+    const auto dij_metric = second_spatial_derivative_metric(
+        metric, di_metric, di_imetric, dij_imetric);
+    const auto di_christoffel = spatial_derivative_christoffel(
+        di_metric, dij_metric, imetric, di_imetric);
+    const auto& u0 = get(get<Tags::TimeDilationFactor>(background));
+
+    // The evolution of the mass is a second order effect so we only include it
+    // starting at the 2nd iteration.
+    const double evolved_mass = iteration == 1
+                                    ? mass.value()
+                                    : mass.value() - charge * get(psi_monopole);
+
+    const auto self_force_acc = self_force_acceleration(
+        dt_psi_monopole, psi_dipole, particle_velocity, charge, evolved_mass,
+        imetric, get<Tags::TimeDilationFactor>(background));
+
+    const double t_minus_turnup = time - turn_on_time.value();
+    roll_on = turn_on_function(t_minus_turnup, turn_on_interval.value());
+    const double dt_roll_on =
+        dt_turn_on_function(t_minus_turnup, turn_on_interval.value());
+    const double dt2_roll_on =
+        dt2_turn_on_function(t_minus_turnup, turn_on_interval.value());
+    for (size_t i = 0; i < Dim; ++i) {
+      acc.get(i) += roll_on * self_force_acc.get(i);
+    }
+    const tnsr::A<double, Dim> u{{u0, get<0>(particle_velocity) * u0,
+                                  get<1>(particle_velocity) * u0,
+                                  get<2>(particle_velocity) * u0}};
+    const tnsr::a<double, Dim> d_psiR{{get(dt_psi_monopole), get<0>(psi_dipole),
+                                       get<1>(psi_dipole), get<2>(psi_dipole)}};
+
+    // Eq.(55)
+    double dt2_psiR = get<0>(trace_christoffel) * get<0>(d_psiR);
+    for (size_t i = 0; i < Dim; ++i) {
+      dt2_psiR += -2. * imetric.get(0, i + 1) * dt_psi_dipole.get(i) +
+                  trace_christoffel.get(i + 1) * psi_dipole.get(i);
+    }
+    dt2_psiR /= get<0, 0>(imetric);
+
+    tnsr::a<double, Dim> dt_d_psiR{{dt2_psiR, get<0>(dt_psi_dipole),
+                                    get<1>(dt_psi_dipole),
+                                    get<2>(dt_psi_dipole)}};
+    for (size_t i = 0; i < Dim; ++i) {
+      get<0>(dt_d_psiR) += particle_velocity.get(i) * dt_psi_dipole.get(i);
+    }
+
+    double dt_evolved_mass = get(dt_psi_monopole);
+    double dt2_evolved_mass = dt2_psiR;
+    for (size_t i = 0; i < Dim; ++i) {
+      dt_evolved_mass += psi_dipole.get(i) * particle_velocity.get(i);
+      dt2_evolved_mass += 2. * dt_psi_dipole.get(i) * particle_velocity.get(i) +
+                          psi_dipole.get(i) * acc.get(i);
+    }
+    dt_evolved_mass *= iteration == 1 ? 0. : -charge;
+    dt2_evolved_mass *= iteration == 1 ? 0. : -charge;
+
+    const double dt_mass_factor = -dt_evolved_mass / evolved_mass;
+    const double dt2_mass_factor = (2. * dt_evolved_mass * dt_evolved_mass -
+                                    evolved_mass * dt2_evolved_mass) /
+                                   square(evolved_mass);
+    const auto dt_christoffel = tenex::evaluate<ti::A, ti::b, ti::c>(
+        particle_velocity(ti::I) * di_christoffel(ti::i, ti::A, ti::b, ti::c));
+
+    const auto dt_imetric = tenex::evaluate<ti::A, ti::B>(
+        particle_velocity(ti::I) * di_imetric(ti::i, ti::A, ti::B));
+    const auto dt2_imetric = tenex::evaluate<ti::A, ti::B>(
+        particle_velocity(ti::I) * particle_velocity(ti::J) *
+            dij_imetric(ti::i, ti::j, ti::A, ti::B) +
+        acc(ti::I) * di_imetric(ti::i, ti::A, ti::B));
+    // Eq.(51a)
+    const auto dt_u = tenex::evaluate<ti::A>(
+        charge / evolved_mass / u0 * imetric(ti::A, ti::B) * d_psiR(ti::b) -
+        christoffel(ti::A, ti::b, ti::c) * u(ti::B) * u(ti::C) / u0);
+    // Eq.(51b)
+    const auto dt2_u = tenex::evaluate<ti::A>(
+        charge / evolved_mass / u0 *
+            (dt_imetric(ti::A, ti::B) * d_psiR(ti::b) +
+             imetric(ti::A, ti::B) * dt_d_psiR(ti::b) +
+             dt_mass_factor * imetric(ti::A, ti::B) * d_psiR(ti::b)) -
+        (dt_christoffel(ti::A, ti::b, ti::c) * u(ti::B) * u(ti::C) +
+         2. * christoffel(ti::A, ti::b, ti::c) * dt_u(ti::B) * u(ti::C) +
+         get<0>(dt_u) * dt_u(ti::A)) /
+            u0);
+    const auto dt_u_rollon = tenex::evaluate<ti::A>(roll_on * dt_u(ti::A));
+    const auto dt2_u_rollon = tenex::evaluate<ti::A>(dt_roll_on * dt_u(ti::A) +
+                                                     roll_on * dt2_u(ti::A));
+
+    // Eq.(56)
+    tnsr::i<double, Dim> d_dt2_psiR{0.};
+    for (size_t i = 0; i < Dim; ++i) {
+      d_dt2_psiR.get(i) +=
+          -di_imetric.get(i, 0, 0) * dt2_psiR +
+          di_trace_christoffel.get(i, 0) * get(dt_psi_monopole) +
+          trace_christoffel.get(0) * dt_psi_dipole.get(i);
+      for (size_t j = 0; j < Dim; ++j) {
+        d_dt2_psiR.get(i) +=
+            -2. * di_imetric.get(i, 0, j + 1) * dt_psi_dipole.get(j) +
+            di_trace_christoffel.get(i, j + 1) * psi_dipole.get(j);
+      }
+      d_dt2_psiR.get(i) /= get<0, 0>(imetric);
+    }
+
+    // Eq.(57)
+    double dt3_psiR = trace_christoffel.get(0) * dt2_psiR;
+    for (size_t i = 0; i < Dim; ++i) {
+      dt3_psiR += -2. * imetric.get(0, i + 1) * d_dt2_psiR.get(i) +
+                  trace_christoffel.get(i + 1) * dt_psi_dipole.get(i);
+    }
+    dt3_psiR /= get<0, 0>(imetric);
+    tnsr::a<double, Dim> dt2_d_psiR;
+    dt2_d_psiR.get(0) = dt3_psiR;
+    for (size_t i = 0; i < Dim; ++i) {
+      dt2_d_psiR.get(0) += 2. * d_dt2_psiR.get(i) * particle_velocity.get(i) +
+                           dt_psi_dipole.get(i) * acc.get(i);
+      dt2_d_psiR.get(i + 1) = d_dt2_psiR.get(i);
+    }
+    auto f = self_force_per_mass(d_psiR, u, charge, evolved_mass, imetric);
+    auto dt_f =
+        dt_self_force_per_mass(d_psiR, dt_d_psiR, u, dt_u_rollon, charge,
+                               evolved_mass, imetric, dt_imetric);
+    auto dt2_f = dt2_self_force_per_mass(
+        d_psiR, dt_d_psiR, dt2_d_psiR, u, dt_u_rollon, dt2_u_rollon, charge,
+        evolved_mass, imetric, dt_imetric, dt2_imetric);
+    for (size_t i = 0; i < 4; ++i) {
+      dt_f.get(i) += dt_mass_factor * f.get(i);
+      dt2_f.get(i) +=
+          2. * dt_mass_factor * dt_f.get(i) + dt2_mass_factor * f.get(i);
+    }
+    const auto f_roll_on = tenex::evaluate<ti::A>(roll_on * f(ti::A));
+    const auto dt_f_roll_on =
+        tenex::evaluate<ti::A>(dt_roll_on * f(ti::A) + roll_on * dt_f(ti::A));
+    const auto dt2_f_roll_on = tenex::evaluate<ti::A>(
+        dt2_roll_on * f(ti::A) + 2. * dt_roll_on * dt_f(ti::A) +
+        roll_on * dt2_f(ti::A));
+
+    const auto cov_f =
+        Du_self_force_per_mass(f_roll_on, dt_f_roll_on, u, christoffel);
+    const auto dt_cov_f =
+        dt_Du_self_force_per_mass(f_roll_on, dt_f_roll_on, dt2_f_roll_on, u,
+                                  dt_u_rollon, christoffel, dt_christoffel);
+    for (size_t i = 0; i < Dim; ++i) {
+      get(*acceleration_terms)[i + 3] = f_roll_on.get(i);
+      get(*acceleration_terms)[i + 6] = dt_f_roll_on.get(i);
+      get(*acceleration_terms)[i + 9] = cov_f.get(i);
+      get(*acceleration_terms)[i + 12] = dt_cov_f.get(i);
+    }
+  }
   for (size_t i = 0; i < Dim; ++i) {
-    get(*acceleration_terms)[i] = geodesic_acc.get(i) + self_force_acc.get(i);
+    get(*acceleration_terms)[i] = acc.get(i);
   }
 }
 }  // namespace CurvedScalarWave::Worldtube

src/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/IterateAccelerationTerms.hpp

@@ -17,7 +17,11 @@ namespace CurvedScalarWave::Worldtube {
 
 /*!
  * \brief Computes the next iteration of the acceleration due to scalar self
- * force from the current iteration of the regular field.
+ * force from the current iteration of the regular field, as well as the
+ * quantities required to compute the acceleration terms of the puncture field.
+ *
+ * \details Analytic expressions for the computed terms are given in Section V.B
+ * of \cite Wittek:2024gxn.
  */
 struct IterateAccelerationTerms {
   static constexpr size_t Dim = 3;
@@ -30,10 +34,12 @@ struct IterateAccelerationTerms {
       Stf::Tags::StfTensor<::Tags::dt<Tags::PsiWorldtube>, 0, Dim,
                            Frame::Inertial>,
       Stf::Tags::StfTensor<Tags::PsiWorldtube, 1, Dim, Frame::Inertial>,
+      Stf::Tags::StfTensor<::Tags::dt<Tags::PsiWorldtube>, 1, Dim,
+                           Frame::Inertial>,
       Tags::Charge, Tags::Mass, ::Tags::Time, Tags::SelfForceTurnOnTime,
-      Tags::SelfForceTurnOnInterval>;
+      Tags::SelfForceTurnOnInterval, Tags::CurrentIteration>;
   static void apply(
-      const gsl::not_null<Scalar<DataVector>*> acceleration_terms,
+      gsl::not_null<Scalar<DataVector>*> acceleration_terms,
       const std::array<tnsr::I<double, Dim>, 2>& pos_vel,
       const tuples::TaggedTuple<
           gr::Tags::SpacetimeMetric<double, Dim>,
@@ -43,10 +49,11 @@ struct IterateAccelerationTerms {
           Tags::TimeDilationFactor>& background,
       const tnsr::I<double, Dim, Frame::Inertial>& geodesic_acc,
       const Scalar<double>& psi_monopole, const Scalar<double>& dt_psi_monopole,
-      const tnsr::i<double, Dim, Frame::Inertial>& psi_dipole, double charge,
+      const tnsr::i<double, Dim, Frame::Inertial>& psi_dipole,
+      const tnsr::i<double, Dim, Frame::Inertial>& dt_psi_dipole, double charge,
       std::optional<double> mass, double time,
       std::optional<double> turn_on_time,
-      std::optional<double> turn_on_interval);
+      std::optional<double> turn_on_interval, size_t iteration);
 };
 
 }  // namespace CurvedScalarWave::Worldtube

tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt

@@ -18,6 +18,7 @@ set(LIBRARY_SOURCES
   SingletonActions/Test_ChangeSlabSize.cpp
   SingletonActions/Test_InitializeElementFacesGridCoordinates.cpp
   SingletonActions/Test_InitializeEvolvedVariables.cpp
+  SingletonActions/Test_IterateAccelerationTerms.cpp
   SingletonActions/Test_ObserveWorldtubeSolution.cpp
   SingletonActions/Test_UpdateAcceleration.cpp
   SingletonActions/Test_UpdateFunctionsOfTime.cpp