Add turn on function to self force acceleration

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

@@ -25,25 +25,31 @@ void UpdateAcceleration::apply(
         gr::Tags::TraceSpacetimeChristoffelSecondKind<double, Dim>,
         Tags::TimeDilationFactor>& background,
     const tnsr::I<double, Dim, Frame::Inertial>& geodesic_acc,
-    const Scalar<double>& dt_psi_monopole,
+    const Scalar<double>& psi_monopole, const Scalar<double>& dt_psi_monopole,
     const tnsr::i<double, Dim, Frame::Inertial>& psi_dipole,
     const double charge, const std::optional<double> mass,
-    const size_t max_iterations) {
+    const size_t max_iterations, const double time,
+    const std::optional<double> turn_on_time,
+    const std::optional<double> turn_on_interval) {
   tnsr::I<double, Dim> self_force_acc(0.);
   const auto& particle_velocity = pos_vel.at(1);
-  if (max_iterations > 0) {
+  double roll_on = 0.;
+  if (max_iterations > 0 and time > turn_on_time.value()) {
+    roll_on =
+        turn_on_function(time - turn_on_time.value(), turn_on_interval.value());
     const auto& inverse_metric =
         get<gr::Tags::InverseSpacetimeMetric<double, Dim>>(background);
     const auto& dilation_factor = get<Tags::TimeDilationFactor>(background);
+    const double evolved_mass = mass.value() - charge * get(psi_monopole);
     self_force_acceleration(make_not_null(&self_force_acc), dt_psi_monopole,
-                            psi_dipole, particle_velocity, charge, mass.value(),
+                            psi_dipole, particle_velocity, charge, evolved_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) + self_force_acc.get(i);
+        geodesic_acc.get(i) + roll_on * self_force_acc.get(i);
   }
 }
 }  // namespace CurvedScalarWave::Worldtube

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

@@ -36,10 +36,12 @@ struct UpdateAcceleration {
   using argument_tags = tmpl::list<
       Tags::ParticlePositionVelocity<Dim>, Tags::BackgroundQuantities<Dim>,
       Tags::GeodesicAcceleration<Dim>,
+      Stf::Tags::StfTensor<Tags::PsiWorldtube, 0, Dim, Frame::Inertial>,
       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::MaxIterations>;
+      Tags::Charge, Tags::Mass, Tags::MaxIterations, ::Tags::Time,
+      Tags::SelfForceTurnOnTime, Tags::SelfForceTurnOnInterval>;
   static void apply(
       gsl::not_null<
           Variables<tmpl::list<::Tags::dt<Tags::EvolvedPosition<Dim>>,
@@ -53,9 +55,11 @@ struct UpdateAcceleration {
           gr::Tags::TraceSpacetimeChristoffelSecondKind<double, Dim>,
           Tags::TimeDilationFactor>& background,
       const tnsr::I<double, Dim, Frame::Inertial>& geodesic_acc,
-      const Scalar<double>& dt_psi_monopole,
+      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, size_t max_iterations);
+      std::optional<double> mass, size_t max_iterations, double time,
+      std::optional<double> turn_on_time,
+      std::optional<double> turn_on_interval);
 };
 
 }  // namespace CurvedScalarWave::Worldtube

tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_Iterations.cpp

@@ -121,7 +121,8 @@ struct MockMetavariables {
       domain::Tags::Domain<Dim>,
       CurvedScalarWave::Tags::BackgroundSpacetime<gr::Solutions::KerrSchild>,
       Tags::ExcisionSphere<Dim>, Tags::WorldtubeRadius, Tags::ExpansionOrder,
-      Tags::MaxIterations, Tags::Charge, Tags::Mass>;
+      Tags::MaxIterations, Tags::Charge, Tags::Mass, ::Tags::Time,
+      Tags::SelfForceTurnOnTime, Tags::SelfForceTurnOnInterval>;
 };
 
 void test_iterations(const size_t max_iterations) {
@@ -163,19 +164,25 @@ void test_iterations(const size_t max_iterations) {
 
     const auto& initial_refinements = shell.initial_refinement_levels();
     const auto& initial_extents = shell.initial_extents();
+    // unused but the tag is needed to compile
+    const double time = std::numeric_limits<double>::signaling_NaN();
     tuples::TaggedTuple<
         domain::Tags::Domain<Dim>,
         CurvedScalarWave::Tags::BackgroundSpacetime<gr::Solutions::KerrSchild>,
         Tags::ExcisionSphere<Dim>, Tags::WorldtubeRadius, Tags::ExpansionOrder,
-        Tags::MaxIterations, Tags::Charge, Tags::Mass>
+        Tags::MaxIterations, Tags::Charge, Tags::Mass, ::Tags::Time,
+        Tags::SelfForceTurnOnTime, Tags::SelfForceTurnOnInterval>
         tuple_of_opts{shell.create_domain(),
                       kerr_schild,
                       excision_sphere,
                       excision_sphere.radius(),
                       expansion_order,
                       max_iterations,
                       charge,
-                      std::make_optional(mass)};
+                      std::make_optional(mass),
+                      time,
+                      std::nullopt,
+                      std::nullopt};
     ActionTesting::MockRuntimeSystem<metavars> runner{std::move(tuple_of_opts)};
     const auto element_ids = initial_element_ids(initial_refinements);
     const auto& blocks = shell_domain.blocks();

tests/Unit/Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/Test_SendToWorldtube.cpp

@@ -123,7 +123,8 @@ struct MockMetavariables {
   using const_global_cache_tags =
       tmpl::list<domain::Tags::Domain<Dim>, Tags::ExcisionSphere<Dim>,
                  Tags::WorldtubeRadius, Tags::ExpansionOrder,
-                 Tags::MaxIterations, Tags::Charge, Tags::Mass>;
+                 Tags::MaxIterations, Tags::Charge, Tags::Mass, ::Tags::Time,
+                 Tags::SelfForceTurnOnTime, Tags::SelfForceTurnOnInterval>;
 };
 
 // This test checks that `SendToWorldtube` integrates the regular field on the
@@ -142,6 +143,8 @@ SPECTRE_TEST_CASE("Unit.CurvedScalarWave.Worldtube.SendToWorldtube", "[Unit]") {
   const size_t initial_extent = 10;
   const size_t face_size = initial_extent * initial_extent;
   const auto quadrature = Spectral::Quadrature::GaussLobatto;
+  // unused but the tag is needed to compile
+  const double time = std::numeric_limits<double>::signaling_NaN();
   // we create several differently refined shells so a different number of
   // elements sends data
   for (const auto& [expansion_order, initial_refinement, worldtube_radius] :
@@ -166,13 +169,18 @@ SPECTRE_TEST_CASE("Unit.CurvedScalarWave.Worldtube.SendToWorldtube", "[Unit]") {
     // self force and therefore iterative scheme is turned off
     tuples::TaggedTuple<domain::Tags::Domain<Dim>, Tags::ExcisionSphere<Dim>,
                         Tags::WorldtubeRadius, Tags::ExpansionOrder,
-                        Tags::MaxIterations, Tags::Charge, Tags::Mass>
+                        Tags::MaxIterations, Tags::Charge, Tags::Mass,
+                        ::Tags::Time, Tags::SelfForceTurnOnTime,
+                        Tags::SelfForceTurnOnInterval>
         tuple_of_opts{shell.create_domain(),
                       excision_sphere,
                       excision_sphere.radius(),
                       expansion_order,
                       static_cast<size_t>(0),
                       0.1,
+                      std::nullopt,
+                      time,
+                      std::nullopt,
                       std::nullopt};
     ActionTesting::MockRuntimeSystem<metavars> runner{std::move(tuple_of_opts)};
     const auto element_ids = initial_element_ids(initial_refinements);

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

@@ -21,6 +21,7 @@
 #include "Helpers/DataStructures/MakeWithRandomValues.hpp"
 #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
 #include "ParallelAlgorithms/Actions/MutateApply.hpp"
+#include "Time/Tags/Time.hpp"
 #include "Utilities/TMPL.hpp"
 
 namespace CurvedScalarWave::Worldtube {
@@ -56,20 +57,33 @@ SPECTRE_TEST_CASE("Unit.Evolution.Systems.CSW.Worldtube.UpdateAcceleration",
       make_with_random_values<Scalar<double>>(make_not_null(&gen), dist, 1);
   const Tags::BackgroundQuantities<Dim>::type background_quantities{
       metric, inverse_metric, christoffel, trace_christoffel, dilation};
+  const auto psi_monopole =
+      make_with_random_values<Scalar<double>>(make_not_null(&gen), dist, 1);
   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 double mass = 0.1;
+  const double charge = 0.2;
+
+  const double time = 105.;
+  const double turn_on_time = 100.;
+  const double turn_on_interval = 10.;
   const size_t max_iterations = 0;
   auto box = db::create<db::AddSimpleTags<
       dt_variables_tag, Tags::ParticlePositionVelocity<Dim>,
       Tags::BackgroundQuantities<Dim>, Tags::GeodesicAcceleration<Dim>,
+      Stf::Tags::StfTensor<Tags::PsiWorldtube, 0, Dim, Frame::Inertial>,
       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::MaxIterations>>(
+      Tags::Charge, Tags::Mass, Tags::MaxIterations, ::Tags::Time,
+      Tags::SelfForceTurnOnTime, Tags::SelfForceTurnOnInterval>>(
       std::move(dt_evolved_vars), pos_vel, background_quantities, geodesic_acc,
-      dt_psi_monopole, psi_dipole, 1., std::make_optional(1.), max_iterations);
+      psi_monopole, dt_psi_monopole, psi_dipole, charge,
+      std::make_optional(mass), max_iterations, time,
+      std::make_optional(turn_on_time), std::make_optional(turn_on_interval));
 
   db::mutate_apply<UpdateAcceleration>(make_not_null(&box));
   const auto& dt_vars = db::get<dt_variables_tag>(box);
@@ -87,13 +101,18 @@ SPECTRE_TEST_CASE("Unit.Evolution.Systems.CSW.Worldtube.UpdateAcceleration",
       make_not_null(&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);
+
+  const double roll_on =
+      turn_on_function(time - turn_on_time, turn_on_interval);
+  const double evolved_mass = mass - charge * get(psi_monopole);
+  const auto self_force_acc =
+      self_force_acceleration(dt_psi_monopole, psi_dipole, vel, charge,
+                              evolved_mass, inverse_metric, dilation);
   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) + self_force_acc.get(i));
+          geodesic_acc.get(i) + roll_on * self_force_acc.get(i));
   }
 }
 }  // namespace