fixup

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

@@ -36,7 +36,6 @@ spectre_target_headers(
 target_link_libraries(
   ${LIBRARY}
   PUBLIC
-  ControlSystem
   DataStructures
   Domain
   GeneralRelativity

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

@@ -33,17 +33,17 @@ struct InitializeEvolvedVariables {
 
   using simple_tags =
       tmpl::list<variables_tag, dt_variables_tag, Tags::CurrentIteration,
-                 Tags::ExpirationTime,
+                 Tags::ExpirationTime, Tags::WorldtubeRadius,
                  ::Tags::HistoryEvolvedVariables<variables_tag>>;
   using return_tags = simple_tags;
 
   using compute_tags = tmpl::list<>;
   using const_global_cache_tags = tmpl::list<>;
   using mutable_global_cache_tags = tmpl::list<>;
   using simple_tags_from_options = tmpl::list<Tags::InitialPositionAndVelocity>;
-  using argument_tags =
-      tmpl::list<::Tags::TimeStepper<TimeStepper>,
-                 Tags::InitialPositionAndVelocity, ::Tags::Time>;
+  using argument_tags = tmpl::list<::Tags::TimeStepper<TimeStepper>,
+                                   Tags::InitialPositionAndVelocity,
+                                   ::Tags::Time, Tags::ExcisionSphere<Dim>>;
   static void apply(
       const gsl::not_null<Variables<
           tmpl::list<Tags::EvolvedPosition<Dim>, Tags::EvolvedVelocity<Dim>>>*>
@@ -54,13 +54,18 @@ struct InitializeEvolvedVariables {
           dt_evolved_vars,
       const gsl::not_null<size_t*> current_iteration,
       const gsl::not_null<double*> expiration_time,
+      const gsl::not_null<double*> worldtube_radius,
       const gsl::not_null<::Tags::HistoryEvolvedVariables<variables_tag>::type*>
           time_stepper_history,
       const TimeStepper& time_stepper,
       const std::array<tnsr::I<double, Dim>, 2>& initial_pos_and_vel,
-      const double initial_time) {
+      const double initial_time, const ExcisionSphere<Dim> excision_sphere) {
     *current_iteration = 0;
+
+    // the functions of time should be ready during the first step but not the
+    // second. We choose the arbitrary value of 1e-10 here to ensure this.
     *expiration_time = initial_time + 1e-10;
+    *worldtube_radius = excision_sphere.radius();
 
     const size_t starting_order =
         time_stepper.number_of_past_steps() == 0 ? time_stepper.order() : 1;

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

@@ -102,7 +102,7 @@ struct ReceiveElementData {
     for (const auto& [_, element_ylm_coefs] : inbox.at(time_step_id)) {
       external_ylm_coefs += element_ylm_coefs;
     }
-    const double wt_radius = db::get<Tags::ExcisionSphere<Dim>>(box).radius();
+    const double wt_radius = db::get<Tags::WorldtubeRadius>(box);
     external_ylm_coefs /= wt_radius * wt_radius;
 
     DataVector& psi_ylm_coefs =

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

@@ -26,6 +26,11 @@
 #include "Utilities/Gsl.hpp"
 #include "Utilities/TMPL.hpp"
 
+// forward declare
+namespace control_system {
+struct UpdateMultipleFunctionsOfTime;
+}
+
 namespace CurvedScalarWave::Worldtube::Actions {
 
 /*!
@@ -61,15 +66,10 @@ struct UpdateFunctionsOfTime {
     const double grid_radius_particle =
         get(magnitude(excision_sphere.center()));
 
-    DataVector angular_update(2);
-    DataVector expansion_update(2);
-    DataVector size_a_update(2);
-    DataVector size_b_update(2);
-
-    angular_update[0] = atan2(y, x);
-    angular_update[1] = (x * ydot - y * xdot) / square(r);
-    expansion_update[0] = r / grid_radius_particle;
-    expansion_update[1] = radial_vel / grid_radius_particle;
+    const DataVector angular_update{atan2(y, x),
+                                    (x * ydot - y * xdot) / square(r)};
+    const DataVector expansion_update{r / grid_radius_particle,
+                                      radial_vel / grid_radius_particle};
 
     const auto& wt_radius_params =
         db::get<Tags::WorldtubeRadiusParameters>(box);
@@ -101,21 +101,20 @@ struct UpdateFunctionsOfTime {
     const double sqrt_4_pi = sqrt(4. * M_PI);
     // The expansion map has already stretched the excision spheres and we need
     // to account for that.
-    size_a_update[0] =
-        sqrt_4_pi * (wt_radius_grid - wt_radius_inertial / expansion_update[0]);
-    size_a_update[1] = sqrt_4_pi *
-                       (wt_radius_inertial * expansion_update[1] -
-                        wt_radius_time_derivative * expansion_update[0]) /
-                       square(expansion_update[0]);
+    const DataVector size_a_update{
+        sqrt_4_pi * (wt_radius_grid - wt_radius_inertial / expansion_update[0]),
+        sqrt_4_pi *
+            (wt_radius_inertial * expansion_update[1] -
+             wt_radius_time_derivative * expansion_update[0]) /
+            square(expansion_update[0])};
 
-    size_b_update[0] =
+    const DataVector size_b_update{
         sqrt_4_pi * (bh_excision_radius_grid -
-                     bh_excision_radius_inertial / expansion_update[0]);
-    size_b_update[1] =
+                     bh_excision_radius_inertial / expansion_update[0]),
         sqrt_4_pi *
-        (bh_excision_radius_inertial * expansion_update[1] -
-         bh_excision_radius_time_derivative * expansion_update[0]) /
-        square(expansion_update[0]);
+            (bh_excision_radius_inertial * expansion_update[1] -
+             bh_excision_radius_time_derivative * expansion_update[0]) /
+            square(expansion_update[0])};
 
     // we set the new expiration time to 1/100th of a time step next to the
     // current time. This is small enough that it can handle rapid time step
@@ -138,12 +137,11 @@ struct UpdateFunctionsOfTime {
     std::unordered_map<std::string, std::pair<DataVector, double>>
         all_updates{};
     all_updates["Rotation"] =
-        std::make_pair(std::move(angular_update), new_expiration_time);
+        std::make_pair(angular_update, new_expiration_time);
     all_updates["Expansion"] =
-        std::make_pair(std::move(expansion_update), new_expiration_time);
+        std::make_pair(expansion_update, new_expiration_time);
     all_updates["SizeA"] = std::make_pair(size_a_update, new_expiration_time);
-    all_updates["SizeB"] =
-        std::make_pair(std::move(size_b_update), new_expiration_time);
+    all_updates["SizeB"] = std::make_pair(size_b_update, new_expiration_time);
 
     Parallel::mutate<::domain::Tags::FunctionsOfTime,
                      control_system::UpdateMultipleFunctionsOfTime>(

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

@@ -28,6 +28,7 @@ add_test_library(${LIBRARY} "${LIBRARY_SOURCES}")
 target_link_libraries(
   ${LIBRARY}
   PRIVATE
+  ControlSystem
   CurvedScalarWave
   CurvedScalarWaveHelpers
   DataStructures

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

@@ -120,8 +120,8 @@ struct MockMetavariables {
   using const_global_cache_tags = tmpl::list<
       domain::Tags::Domain<Dim>,
       CurvedScalarWave::Tags::BackgroundSpacetime<gr::Solutions::KerrSchild>,
-      Tags::ExcisionSphere<Dim>, Tags::ExpansionOrder, Tags::MaxIterations,
-      Tags::Charge, Tags::Mass>;
+      Tags::ExcisionSphere<Dim>, Tags::WorldtubeRadius, Tags::ExpansionOrder,
+      Tags::MaxIterations, Tags::Charge, Tags::Mass>;
 };
 
 void test_iterations(const size_t max_iterations) {
@@ -166,10 +166,15 @@ void test_iterations(const size_t max_iterations) {
     tuples::TaggedTuple<
         domain::Tags::Domain<Dim>,
         CurvedScalarWave::Tags::BackgroundSpacetime<gr::Solutions::KerrSchild>,
-        Tags::ExcisionSphere<Dim>, Tags::ExpansionOrder, Tags::MaxIterations,
-        Tags::Charge, Tags::Mass>
-        tuple_of_opts{shell.create_domain(),   kerr_schild,    excision_sphere,
-                      expansion_order,         max_iterations, charge,
+        Tags::ExcisionSphere<Dim>, Tags::WorldtubeRadius, Tags::ExpansionOrder,
+        Tags::MaxIterations, Tags::Charge, Tags::Mass>
+        tuple_of_opts{shell.create_domain(),
+                      kerr_schild,
+                      excision_sphere,
+                      excision_sphere.radius(),
+                      expansion_order,
+                      max_iterations,
+                      charge,
                       std::make_optional(mass)};
     ActionTesting::MockRuntimeSystem<metavars> runner{std::move(tuple_of_opts)};
     const auto element_ids = initial_element_ids(initial_refinements);

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

@@ -122,8 +122,8 @@ struct MockMetavariables {
   using dg_element_array = MockElementArray<MockMetavariables>;
   using const_global_cache_tags =
       tmpl::list<domain::Tags::Domain<Dim>, Tags::ExcisionSphere<Dim>,
-                 Tags::ExpansionOrder, Tags::MaxIterations, Tags::Charge,
-                 Tags::Mass>;
+                 Tags::WorldtubeRadius, Tags::ExpansionOrder,
+                 Tags::MaxIterations, Tags::Charge, Tags::Mass>;
 };
 
 // This test checks that `SendToWorldtube` integrates the regular field on the
@@ -165,10 +165,11 @@ SPECTRE_TEST_CASE("Unit.CurvedScalarWave.Worldtube.SendToWorldtube", "[Unit]") {
     const auto& initial_extents = shell.initial_extents();
     // self force and therefore iterative scheme is turned off
     tuples::TaggedTuple<domain::Tags::Domain<Dim>, Tags::ExcisionSphere<Dim>,
-                        Tags::ExpansionOrder, Tags::MaxIterations, Tags::Charge,
-                        Tags::Mass>
+                        Tags::WorldtubeRadius, Tags::ExpansionOrder,
+                        Tags::MaxIterations, Tags::Charge, Tags::Mass>
         tuple_of_opts{shell.create_domain(),
                       excision_sphere,
+                      excision_sphere.radius(),
                       expansion_order,
                       static_cast<size_t>(0),
                       0.1,

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

@@ -70,6 +70,10 @@ struct MockMetavariables {
                  Tags::BlackHoleRadiusParameters>;
 };
 
+// this test takes a random position and velocity of the scalar charge and then
+// calls `UpdateFunctionsOfTime`. After this action, the functions of time
+// should be set in such a way that the worldtube excision sphere has the same
+// position, velocity as the particle with the appropriate worldtube radius.
 void test_particle_motion() {
   static constexpr size_t Dim = 3;
   using metavars = MockMetavariables<Dim>;
@@ -88,6 +92,7 @@ void test_particle_motion() {
   const auto bh_excision_sphere =
       domain.excision_spheres().at("ExcisionSphereB");
 
+  // these will be overwritten
   const std::array<double, 4> wt_radius_options{{0.1, 5., 1e-2, 2.3}};
   const std::array<double, 4> bh_radius_options{{0.1, 5., 1e-2, 0.123}};
 
@@ -156,13 +161,18 @@ void test_particle_motion() {
   // includes worldtube shape map
   const auto& grid_to_inertial_map_wt_boundary =
       wt_neighbor_block.moving_mesh_grid_to_inertial_map();
+
+  // checks the position and velocity of the worldtube excision sphere are now
+  // mapped according to the position and velocity of the particle
   const auto mapped_particle =
       grid_to_inertial_map_particle.coords_frame_velocity_jacobians(
           wt_excision_sphere.center(), 0.5, fots);
   CHECK_ITERABLE_APPROX(std::get<0>(mapped_particle), particle_pos_vel.at(0));
   CHECK_ITERABLE_APPROX(std::get<3>(mapped_particle), particle_pos_vel.at(1));
   const double orbit_radius = get(magnitude(particle_pos_vel[0]));
 
+  // check that the grid points on the worldtube boundary are mapped according
+  // to the prescribed worldtube radius
   auto wt_boundary_point = wt_excision_sphere.center();
   wt_boundary_point[0] += wt_excision_sphere.radius();
   const auto mapped_wt_boundary_data =