Compute scalar self-force in UpdateAcceleration

src/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt

@@ -12,6 +12,7 @@ spectre_target_sources(
   PunctureField.cpp
   PunctureFieldOrder0.cpp
   PunctureFieldOrder1.cpp
+  SelfForce.cpp
   )
 
 spectre_target_headers(
@@ -22,6 +23,7 @@ spectre_target_headers(
   SingletonChare.hpp
   Tags.hpp
   PunctureField.hpp
+  SelfForce.hpp
   Worldtube.hpp
   )
 

src/Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.cpp

@@ -0,0 +1,65 @@
+
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.hpp"
+
+#include <cstddef>
+
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "Utilities/Gsl.hpp"
+
+namespace CurvedScalarWave::Worldtube {
+
+template <size_t Dim>
+void self_force_acceleration(
+    gsl::not_null<tnsr::I<double, Dim>*> self_force_acc,
+    const Scalar<double>& dt_psi_monopole,
+    const tnsr::i<double, Dim>& psi_dipole,
+    const tnsr::I<double, Dim>& particle_velocity, const double particle_charge,
+    const double particle_mass, const tnsr::AA<double, Dim>& inverse_metric,
+    const Scalar<double>& dilation_factor) {
+  const double factor =
+      particle_charge / particle_mass / square(get(dilation_factor));
+  for (size_t i = 0; i < Dim; ++i) {
+    self_force_acc->get(i) =
+        (inverse_metric.get(i + 1, 0) -
+         particle_velocity.get(i) * inverse_metric.get(0, 0)) *
+        get(dt_psi_monopole) * factor;
+    for (size_t j = 0; j < Dim; ++j) {
+      self_force_acc->get(i) +=
+          (inverse_metric.get(i + 1, j + 1) -
+           particle_velocity.get(i) * inverse_metric.get(0, j + 1)) *
+          psi_dipole.get(j) * factor;
+    }
+  }
+}
+
+template <size_t Dim>
+tnsr::I<double, Dim> self_force_acceleration(
+    const Scalar<double>& dt_psi_monopole,
+    const tnsr::i<double, Dim>& psi_dipole,
+    const tnsr::I<double, Dim>& particle_velocity, const double particle_charge,
+    const double particle_mass, const tnsr::AA<double, Dim>& inverse_metric,
+    const Scalar<double>& dilation_factor) {
+  tnsr::I<double, Dim> self_force_acc{};
+  self_force_acceleration(make_not_null(&self_force_acc), dt_psi_monopole,
+                          psi_dipole, particle_velocity, particle_charge,
+                          particle_mass, inverse_metric, dilation_factor);
+  return self_force_acc;
+}
+
+// Instantiations
+template void self_force_acceleration(
+    gsl::not_null<tnsr::I<double, 3>*> self_force_acc,
+    const Scalar<double>& dt_psi_monopole, const tnsr::i<double, 3>& psi_dipole,
+    const tnsr::I<double, 3>& particle_velocity, const double particle_charge,
+    const double particle_mass, const tnsr::AA<double, 3>& inverse_metric,
+    const Scalar<double>& dilation_factor);
+
+template tnsr::I<double, 3> self_force_acceleration(
+    const Scalar<double>& dt_psi_monopole, const tnsr::i<double, 3>& psi_dipole,
+    const tnsr::I<double, 3>& particle_velocity, const double particle_charge,
+    const double particle_mass, const tnsr::AA<double, 3>& inverse_metric,
+    const Scalar<double>& dilation_factor);
+}  // namespace CurvedScalarWave::Worldtube

src/Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.hpp

@@ -0,0 +1,51 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <cstddef>
+
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "Utilities/Gsl.hpp"
+
+namespace CurvedScalarWave::Worldtube {
+
+/// @{
+/*!
+ * \brief Computes the coordinate acceleration due to the scalar self-force onto
+ * the charge.
+ *
+ * \details It is given by
+ *
+ * \begin{equation}
+ * (u^0)^2 \ddot{x}^i_p  = \frac{q}{\mu}(g^{i
+ * \alpha} - \dot{x}^i_p g^{0 \alpha} ) \partial_\alpha \Psi^R
+ * \end{equation}
+ *
+ * where $\dot{x}^i_p$ is the position of the scalar charge, $\Psi^R$ is the
+ * regular field, $q$ is the particle's charge, $\mu$ is the particle's mass,
+ * $u^\alpha$ is the four-velocity and $g^{\alpha \beta}$ is the inverse
+ * spacetime metric in the inertial frame, evaluated at the position of the
+ * particle. An overdot denotes a derivative with respect to coordinate time.
+ * Greek indices are spacetime indices and Latin indices are purely spatial.
+ * Note that the coordinate geodesic acceleration is NOT included.
+ */
+template <size_t Dim>
+void self_force_acceleration(
+    gsl::not_null<tnsr::I<double, Dim>*> self_force_acc,
+    const Scalar<double>& dt_psi_monopole,
+    const tnsr::i<double, Dim>& psi_dipole,
+    const tnsr::I<double, Dim>& particle_velocity, double particle_charge,
+    double particle_mass, const tnsr::AA<double, Dim>& inverse_metric,
+    const Scalar<double>& dilation_factor);
+
+template <size_t Dim>
+tnsr::I<double, Dim> self_force_acceleration(
+    const Scalar<double>& dt_psi_monopole,
+    const tnsr::i<double, Dim>& psi_dipole,
+    const tnsr::I<double, Dim>& particle_velocity, double particle_charge,
+    double particle_mass, const tnsr::AA<double, Dim>& inverse_metric,
+    const Scalar<double>& dilation_factor);
+
+/// @}
+}  // namespace CurvedScalarWave::Worldtube

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

@@ -7,6 +7,7 @@
 
 #include "DataStructures/Tensor/Tensor.hpp"
 #include "DataStructures/Variables.hpp"
+#include "Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.hpp"
 #include "Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp"
 #include "Utilities/Gsl.hpp"
 
@@ -17,13 +18,28 @@ void UpdateAcceleration::apply(
                              ::Tags::dt<Tags::EvolvedVelocity<Dim>>>>*>
         dt_evolved_vars,
     const std::array<tnsr::I<double, Dim>, 2>& pos_vel,
-    const tnsr::I<double, Dim, Frame::Inertial>& geodesic_acc) {
+    const tuples::TaggedTuple<gr::Tags::SpacetimeMetric<double, Dim>,
+                              gr::Tags::InverseSpacetimeMetric<double, Dim>,
+                              Tags::TimeDilationFactor>& background,
+    const tnsr::I<double, Dim, Frame::Inertial>& geodesic_acc,
+    const Scalar<double>& dt_psi_monopole,
+    const tnsr::i<double, Dim, Frame::Inertial>& psi_dipole,
+    const double charge, const double mass, const bool apply_self_force) {
+  tnsr::I<double, Dim> self_force_acc(0.);
   const auto& particle_velocity = pos_vel.at(1);
+  if (apply_self_force) {
+    const auto& inverse_metric =
+        get<gr::Tags::InverseSpacetimeMetric<double, Dim>>(background);
+    const auto& dilation_factor = get<Tags::TimeDilationFactor>(background);
+    self_force_acceleration(make_not_null(&self_force_acc), dt_psi_monopole,
+                            psi_dipole, particle_velocity, charge, mass,
+                            inverse_metric, dilation_factor);
+  }
   for (size_t i = 0; i < Dim; ++i) {
     get<::Tags::dt<Tags::EvolvedPosition<Dim>>>(*dt_evolved_vars).get(i)[0] =
         particle_velocity.get(i);
     get<::Tags::dt<Tags::EvolvedVelocity<Dim>>>(*dt_evolved_vars).get(i)[0] =
-        geodesic_acc.get(i);
+        geodesic_acc.get(i) + self_force_acc.get(i);
   }
 }
 }  // namespace CurvedScalarWave::Worldtube

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

@@ -10,15 +10,21 @@
 #include "DataStructures/Variables.hpp"
 #include "DataStructures/VariablesTag.hpp"
 #include "Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp"
+#include "NumericalAlgorithms/SphericalHarmonics/Tags.hpp"
 #include "Parallel/AlgorithmExecution.hpp"
 #include "Parallel/GlobalCache.hpp"
+#include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
 #include "Utilities/Gsl.hpp"
+#include "Utilities/TaggedTuple.hpp"
 
 namespace CurvedScalarWave::Worldtube {
 
 /*!
- * \brief Computes the geodesic acceleration of the particle, see
- * `Tags::GeodesicAccelerationCompute`. This mutator is run on the worldtube
+ * \brief Computes the final acceleration of the particle at this time step.
+ * \details If `apply_self_force` is `false`, the acceleration will simply be
+ * geodesic, see `gr::geodesic_acceleration`.  Otherwise, the acceleration due
+ * to the scalar self-force is additionally applied to it, see
+ * `self_force_acceleration`. This mutator is run on the worldtube
  * singleton chare.
  */
 struct UpdateAcceleration {
@@ -27,14 +33,26 @@ struct UpdateAcceleration {
       tmpl::list<Tags::EvolvedPosition<Dim>, Tags::EvolvedVelocity<Dim>>>;
   using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>;
   using return_tags = tmpl::list<dt_variables_tag>;
-  using argument_tags = tmpl::list<Tags::ParticlePositionVelocity<Dim>,
-                                   Tags::GeodesicAcceleration<Dim>>;
+  using argument_tags = tmpl::list<
+      Tags::ParticlePositionVelocity<Dim>, Tags::BackgroundQuantities<Dim>,
+      Tags::GeodesicAcceleration<Dim>,
+      Stf::Tags::StfTensor<::Tags::dt<Tags::PsiWorldtube>, 0, Dim,
+                           Frame::Inertial>,
+      Stf::Tags::StfTensor<Tags::PsiWorldtube, 1, Dim, Frame::Inertial>,
+      Tags::Charge, Tags::Mass, Tags::ApplySelfForce>;
   static void apply(
       gsl::not_null<
           Variables<tmpl::list<::Tags::dt<Tags::EvolvedPosition<Dim>>,
                                ::Tags::dt<Tags::EvolvedVelocity<Dim>>>>*>
           dt_evolved_vars,
       const std::array<tnsr::I<double, Dim>, 2>& pos_vel,
-      const tnsr::I<double, Dim, Frame::Inertial>& geodesic_acc);
+      const tuples::TaggedTuple<gr::Tags::SpacetimeMetric<double, Dim>,
+                                gr::Tags::InverseSpacetimeMetric<double, Dim>,
+                                Tags::TimeDilationFactor>& background,
+      const tnsr::I<double, Dim, Frame::Inertial>& geodesic_acc,
+      const Scalar<double>& dt_psi_monopole,
+      const tnsr::i<double, Dim, Frame::Inertial>& psi_dipole, double charge,
+      double mass, bool apply_self_force);
 };
+
 }  // namespace CurvedScalarWave::Worldtube

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

@@ -5,6 +5,7 @@ set(LIBRARY "Test_ScalarWaveWorldtube")
 
 set(LIBRARY_SOURCES
   Test_PunctureField.cpp
+  Test_SelfForce.cpp
   Test_Tags.cpp
   ElementActions/Test_SendToWorldtube.cpp
   ElementActions/Test_ReceiveWorldtubeData.cpp

tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.py

@@ -0,0 +1,23 @@
+# Distributed under the MIT License.
+# See LICENSE.txt for details.
+
+import numpy as np
+
+
+def self_force_acceleration(
+    dt_psi_monopole,
+    psi_dipole,
+    vel,
+    charge,
+    mass,
+    inverse_metric,
+    dilation,
+):
+    # Prepend extra value so dimensions work out for einsum.
+    # The 0th component does not affect the final result
+    four_vel = np.concatenate((np.empty(1), vel), axis=0)
+    d_psi = np.concatenate(([dt_psi_monopole], psi_dipole), axis=0)
+    self_force_acc = np.einsum("ij,j", inverse_metric, d_psi)
+    self_force_acc -= np.einsum("i,j,j", four_vel, inverse_metric[0], d_psi)
+    self_force_acc *= charge / mass / dilation**2
+    return self_force_acc[1:]

tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/Test_UpdateAcceleration.cpp

@@ -12,6 +12,7 @@
 #include "DataStructures/DataBox/PrefixHelpers.hpp"
 #include "DataStructures/Variables.hpp"
 #include "DataStructures/VariablesTag.hpp"
+#include "Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.hpp"
 #include "Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/UpdateAcceleration.hpp"
 #include "Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp"
 #include "Framework/TestCreation.hpp"
@@ -34,28 +35,61 @@ SPECTRE_TEST_CASE("Unit.Evolution.Systems.CSW.Worldtube.UpdateAcceleration",
   MAKE_GENERATOR(gen);
   std::uniform_real_distribution<> dist(-1., 1.);
   const DataVector used_for_size(1);
-  const auto pos = make_with_random_values<tnsr::I<double, 3>>(
-      make_not_null(&gen), dist, used_for_size);
-  const auto vel = make_with_random_values<tnsr::I<double, 3>>(
-      make_not_null(&gen), dist, used_for_size);
   auto dt_evolved_vars = make_with_random_values<dt_variables_tag::type>(
       make_not_null(&gen), dist, used_for_size);
-  const auto acceleration =
-      make_with_random_values<tnsr::I<double, 3>>(make_not_null(&gen), dist, 1);
-  auto box = db::create<
-      db::AddSimpleTags<Tags::ParticlePositionVelocity<Dim>, dt_variables_tag,
-                        Tags::GeodesicAcceleration<Dim>>>(
-      std::array<tnsr::I<double, 3>, 2>{{pos, vel}}, std::move(dt_evolved_vars),
-      acceleration);
+  const auto geodesic_acc = make_with_random_values<tnsr::I<double, Dim>>(
+      make_not_null(&gen), dist, 1);
+  const auto vel = make_with_random_values<tnsr::I<double, Dim>>(
+      make_not_null(&gen), dist, 1);
+  const auto pos = make_with_random_values<tnsr::I<double, Dim>>(
+      make_not_null(&gen), dist, 1);
+  const std::array<tnsr::I<double, Dim>, 2> pos_vel{{pos, vel}};
+  const auto metric = make_with_random_values<tnsr::aa<double, Dim>>(
+      make_not_null(&gen), dist, 1);
+  const auto inverse_metric = make_with_random_values<tnsr::AA<double, Dim>>(
+      make_not_null(&gen), dist, 1);
+  const auto dilation =
+      make_with_random_values<Scalar<double>>(make_not_null(&gen), dist, 1);
+  const Tags::BackgroundQuantities<Dim>::type background_quantities{
+      metric, inverse_metric, dilation};
+  const auto dt_psi_monopole =
+      make_with_random_values<Scalar<double>>(make_not_null(&gen), dist, 1);
+  const auto psi_dipole = make_with_random_values<tnsr::i<double, Dim>>(
+      make_not_null(&gen), dist, 1);
+  const bool apply_self_force = false;
+  auto box = db::create<db::AddSimpleTags<
+      dt_variables_tag, Tags::ParticlePositionVelocity<Dim>,
+      Tags::BackgroundQuantities<Dim>, Tags::GeodesicAcceleration<Dim>,
+      Stf::Tags::StfTensor<::Tags::dt<Tags::PsiWorldtube>, 0, Dim,
+                           Frame::Inertial>,
+      Stf::Tags::StfTensor<Tags::PsiWorldtube, 1, Dim, Frame::Inertial>,
+      Tags::Charge, Tags::Mass, Tags::ApplySelfForce>>(
+      std::move(dt_evolved_vars), pos_vel, background_quantities, geodesic_acc,
+      dt_psi_monopole, psi_dipole, 1., 1., apply_self_force);
 
   db::mutate_apply<UpdateAcceleration>(make_not_null(&box));
+  const auto& dt_vars = db::get<dt_variables_tag>(box);
+  for (size_t i = 0; i < Dim; ++i) {
+    CHECK(get<::Tags::dt<Tags::EvolvedPosition<Dim>>>(dt_vars).get(i)[0] ==
+          vel.get(i));
+    CHECK(get<::Tags::dt<Tags::EvolvedVelocity<Dim>>>(dt_vars).get(i)[0] ==
+          geodesic_acc.get(i));
+  }
+
+  db::mutate<Tags::ApplySelfForce>(
+      [](const gsl::not_null<bool*> apply_self_force_arg) {
+        *apply_self_force_arg = true;
+      },
+      make_not_null(&box));
 
-  const auto& dt_vars_after_mutate = db::get<dt_variables_tag>(box);
+  db::mutate_apply<UpdateAcceleration>(make_not_null(&box));
+  const auto self_force_acc = self_force_acceleration(
+      dt_psi_monopole, psi_dipole, vel, 1., 1., inverse_metric, dilation);
   for (size_t i = 0; i < Dim; ++i) {
-    CHECK(get<::Tags::dt<Tags::EvolvedPosition<Dim>>>(dt_vars_after_mutate)
-              .get(i)[0] == vel.get(i));
-    CHECK(get<::Tags::dt<Tags::EvolvedVelocity<Dim>>>(dt_vars_after_mutate)
-              .get(i)[0] == acceleration.get(i));
+    CHECK(get<::Tags::dt<Tags::EvolvedPosition<Dim>>>(dt_vars).get(i)[0] ==
+          vel.get(i));
+    CHECK(get<::Tags::dt<Tags::EvolvedVelocity<Dim>>>(dt_vars).get(i)[0] ==
+          geodesic_acc.get(i) + self_force_acc.get(i));
   }
 }
 }  // namespace

tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/Test_SelfForce.cpp

@@ -0,0 +1,35 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Framework/TestingFramework.hpp"
+
+#include <random>
+
+#include "Evolution/Systems/CurvedScalarWave/Worldtube/SelfForce.hpp"
+#include "Framework/CheckWithRandomValues.hpp"
+#include "Framework/SetupLocalPythonEnvironment.hpp"
+#include "Framework/TestCreation.hpp"
+#include "Framework/TestHelpers.hpp"
+#include "Utilities/TMPL.hpp"
+
+namespace CurvedScalarWave::Worldtube {
+namespace {
+
+void test_self_force_acceleration() {
+  pypp::SetupLocalPythonEnvironment local_python_env{
+      "Evolution/Systems/CurvedScalarWave/Worldtube"};
+  pypp::check_with_random_values<1>(
+      static_cast<tnsr::I<double, 3> (*)(
+          const Scalar<double>&, const tnsr::i<double, 3>&,
+          const tnsr::I<double, 3>&, const double, const double,
+          const tnsr::AA<double, 3>&, const Scalar<double>&)>(
+          self_force_acceleration<3>),
+      "SelfForce", "self_force_acceleration", {{{-2.0, 2.0}}}, 1);
+}
+
+SPECTRE_TEST_CASE("Unit.Evolution.Systems.CSW.Worldtube.SelfForce",
+                  "[Unit][Evolution]") {
+  test_self_force_acceleration();
+}
+}  // namespace
+}  // namespace CurvedScalarWave::Worldtube