Merge pull request #5847 from nikwit/iteration-tags

Add tags for worldtube iterations

src/Evolution/Executables/CurvedScalarWave/EvolveWorldtubeCurvedScalarWave.hpp

@@ -267,7 +267,7 @@ struct EvolutionMetavars {
       CurvedScalarWave::Worldtube::Tags::ExpansionOrder,
       CurvedScalarWave::Worldtube::Tags::Charge,
       CurvedScalarWave::Worldtube::Tags::Mass,
-      CurvedScalarWave::Worldtube::Tags::ApplySelfForce,
+      CurvedScalarWave::Worldtube::Tags::MaxIterations,
       CurvedScalarWave::Worldtube::Tags::ObserveCoefficientsTrigger>;
 
   using dg_registration_list =

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

@@ -12,6 +12,7 @@ spectre_target_headers(
   INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src
   HEADERS
   InitializeConstraintGammas.hpp
+  InitializeCurrentIteration.hpp
   ReceiveWorldtubeData.hpp
   SendToWorldtube.hpp
 )

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

@@ -0,0 +1,30 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <cstddef>
+
+#include "DataStructures/DataBox/Protocols/Mutator.hpp"
+#include "Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp"
+#include "Utilities/ProtocolHelpers.hpp"
+
+namespace CurvedScalarWave::Worldtube::Initialization {
+
+/*!
+ * \brief Sets the initial value of `CurrentIteration to 0.
+ */
+struct InitializeCurrentIteration : tt::ConformsTo<db::protocols::Mutator> {
+  using return_tags =
+      tmpl::list<CurvedScalarWave::Worldtube::Tags::CurrentIteration>;
+  using argument_tags = tmpl::list<>;
+  using simple_tags = return_tags;
+  using compute_tags = tmpl::list<>;
+  using simple_tags_from_options = tmpl::list<>;
+  using const_global_cache_tags = tmpl::list<>;
+  using mutable_global_cache_tags = tmpl::list<>;
+  static void apply(const gsl::not_null<size_t*> current_iteration) {
+    *current_iteration = 0;
+  }
+};
+}  // namespace CurvedScalarWave::Worldtube::Initialization

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

@@ -2,28 +2,24 @@
 // See LICENSE.txt for details.
 
 #pragma once
-
 #include <cstddef>
-
 #include "DataStructures/DataBox/Prefixes.hpp"
 #include "DataStructures/Variables.hpp"
 #include "DataStructures/VariablesTag.hpp"
 #include "Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp"
 #include "Time/Tags/HistoryEvolvedVariables.hpp"
 #include "Time/TimeSteppers/TimeStepper.hpp"
 #include "Utilities/Gsl.hpp"
-
 /// \cond
 namespace Tags {
 template <typename StepperInterface>
 struct TimeStepper;
 }  // namespace Tags
 /// \endcond
-
 namespace CurvedScalarWave::Worldtube::Initialization {
 /*!
  * \brief Initializes the time stepper and evolved variables used by the
- * worldtube system.
+ * worldtube system. Also sets `Tags::CurrentIteration` to 0.
  *
  * \details Sets the initial position and velocity of the particle to the values
  * specified in the input file. The time stepper history is set analogous to the
@@ -36,7 +32,7 @@ struct InitializeEvolvedVariables {
   using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>;
 
   using simple_tags =
-      tmpl::list<variables_tag, dt_variables_tag,
+      tmpl::list<variables_tag, dt_variables_tag, Tags::CurrentIteration,
                  ::Tags::HistoryEvolvedVariables<variables_tag>>;
   using return_tags = simple_tags;
 
@@ -54,11 +50,12 @@ struct InitializeEvolvedVariables {
           Variables<tmpl::list<::Tags::dt<Tags::EvolvedPosition<Dim>>,
                                ::Tags::dt<Tags::EvolvedVelocity<Dim>>>>*>
           dt_evolved_vars,
+      const gsl::not_null<size_t*> current_iteration,
       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) {
+    *current_iteration = 0;
     const size_t starting_order =
         time_stepper.number_of_past_steps() == 0 ? time_stepper.order() : 1;
     *time_stepper_history =

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

@@ -24,15 +24,16 @@ void UpdateAcceleration::apply(
     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) {
+    const double charge, const std::optional<double> mass,
+    const size_t max_iterations) {
   tnsr::I<double, Dim> self_force_acc(0.);
   const auto& particle_velocity = pos_vel.at(1);
-  if (apply_self_force) {
+  if (max_iterations > 0) {
     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,
+                            psi_dipole, particle_velocity, charge, mass.value(),
                             inverse_metric, dilation_factor);
   }
   for (size_t i = 0; i < Dim; ++i) {

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

@@ -21,7 +21,7 @@ namespace CurvedScalarWave::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
+ * \details If `max_iterations` is 0, 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
@@ -39,7 +39,7 @@ struct UpdateAcceleration {
       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>;
+      Tags::Charge, Tags::Mass, Tags::MaxIterations>;
   static void apply(
       gsl::not_null<
           Variables<tmpl::list<::Tags::dt<Tags::EvolvedPosition<Dim>>,
@@ -52,7 +52,7 @@ struct UpdateAcceleration {
       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);
+      std::optional<double> mass, size_t max_iterations);
 };
 
 }  // namespace CurvedScalarWave::Worldtube

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

@@ -25,6 +25,7 @@
 #include "Evolution/Systems/CurvedScalarWave/BackgroundSpacetime.hpp"
 #include "Evolution/Systems/CurvedScalarWave/Tags.hpp"
 #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
+#include "Options/Auto.hpp"
 #include "Options/String.hpp"
 #include "ParallelAlgorithms/EventsAndTriggers/Trigger.hpp"
 #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp"
@@ -66,34 +67,40 @@ struct Charge {
 };
 
 /*!
- * \brief Options for the scalar self-force
+ * \brief Options for the scalar self-force. Select `None` for a purely geodesic
+ * evolution
  */
 struct SelfForceOptions {
-  static constexpr Options::String help = {"Options for the scalar self-force"};
+  static constexpr Options::String help = {
+      "Options for the scalar self-force. Select `None` for a purely geodesic "
+      "evolution"};
   using group = Worldtube;
-};
+  using type = Options::Auto<SelfForceOptions, Options::AutoLabel::None>;
 
-/*!
- * \brief The mass of the scalar particle in units of the black hole mass M.
- */
-struct Mass {
-  using type = double;
-  static constexpr Options::String help{
-      "The mass of the scalar particlein units of the black hole mass M."};
-  static double lower_bound() { return 0.; }
-  using group = SelfForceOptions;
-};
+  struct Mass {
+    using type = double;
+    static constexpr Options::String help{
+        "The mass of the scalar particle in units of the black hole mass M."};
+    static double lower_bound() { return 0.; }
+  };
 
-/*!
- * \brief Whether to apply the acceleration due to the scalar self-force to the
- * particle
- */
-struct ApplySelfForce {
-  using type = bool;
-  static constexpr Options::String help{
-      "Whether to apply the acceleration due to the scalar self-force to the "
-      "particle"};
-  using group = SelfForceOptions;
+  struct Iterations {
+    using type = size_t;
+    static constexpr Options::String help{
+        "The number of iterations used to compute the particle acceleration. "
+        "Must be at least 1 as 0 iterations corresponds to the geodesic "
+        "acceleration."};
+    static size_t lower_bound() { return 1; }
+  };
+  void pup(PUP::er& p) {
+    p | mass;
+    p | iterations;
+  }
+
+  using options = tmpl::list<Mass, Iterations>;
+
+  double mass{};
+  size_t iterations{};
 };
 
 /*!
@@ -246,26 +253,41 @@ struct Charge : db::SimpleTag {
 };
 
 /*!
- * \brief The mass of the particle.
+ * \brief The mass of the scalar charge. Only has a value if the scalar self
+ * force is applied.
  */
 struct Mass : db::SimpleTag {
-  using type = double;
-  using option_tags = tmpl::list<OptionTags::Mass>;
+  using type = std::optional<double>;
+  using option_tags = tmpl::list<OptionTags::SelfForceOptions>;
   static constexpr bool pass_metavariables = false;
-  static double create_from_options(const double mass) { return mass; };
+  static std::optional<double> create_from_options(
+      const std::optional<OptionTags::SelfForceOptions>& self_force_options) {
+    return self_force_options.has_value()
+               ? std::make_optional(self_force_options->mass)
+               : std::nullopt;
+  }
 };
 
 /*!
- * \brief Whether to apply the acceleration due to the scalar self-force to the
- * particle
+ * \brief The maximum number of iterations that will be applied to the
+ * acceleration of the particle.
  */
-struct ApplySelfForce : db::SimpleTag {
-  using type = bool;
-  using option_tags = tmpl::list<OptionTags::ApplySelfForce>;
+struct MaxIterations : db::SimpleTag {
+  using type = size_t;
+  using option_tags = tmpl::list<OptionTags::SelfForceOptions>;
   static constexpr bool pass_metavariables = false;
-  static bool create_from_options(const bool apply_self_force) {
-    return apply_self_force;
-  };
+  static size_t create_from_options(
+      const std::optional<OptionTags::SelfForceOptions>& self_force_options) {
+    return self_force_options.has_value() ? self_force_options->iterations : 0;
+  }
+};
+
+/*!
+ * \brief The current number of iterations that has been applied to the
+ * acceleration of the particle.
+ */
+struct CurrentIteration : db::SimpleTag {
+  using type = size_t;
 };
 
 /*!

tests/InputFiles/CurvedScalarWave/WorldtubeKerrSchild.yaml

@@ -135,9 +135,7 @@ Worldtube:
   ExcisionSphere: ExcisionSphereA
   ExpansionOrder: 1
   Charge: 0.5
-  SelfForceOptions:
-    ApplySelfForce: false
-    Mass: 0.5
+  SelfForceOptions: None
   ObserveCoefficientsTrigger:
     Slabs:
       EvenlySpaced:

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

@@ -11,6 +11,7 @@ set(LIBRARY_SOURCES
   ElementActions/Test_SendToWorldtube.cpp
   ElementActions/Test_ReceiveWorldtubeData.cpp
   ElementActions/Test_InitializeConstraintGammas.cpp
+  ElementActions/Test_InitializeCurrentIteration.cpp
   SingletonActions/Test_ChangeSlabSize.cpp
   SingletonActions/Test_InitializeElementFacesGridCoordinates.cpp
   SingletonActions/Test_InitializeEvolvedVariables.cpp

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

@@ -0,0 +1,25 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include <array>
+#include <memory>
+
+#include "DataStructures/DataBox/DataBox.hpp"
+#include "Evolution/Systems/CurvedScalarWave/Worldtube/ElementActions/InitializeCurrentIteration.hpp"
+#include "Evolution/Systems/CurvedScalarWave/Worldtube/Tags.hpp"
+#include "Framework/TestingFramework.hpp"
+
+namespace CurvedScalarWave::Worldtube {
+namespace {
+SPECTRE_TEST_CASE(
+    "Unit.Evolution.Systems.CSW.Worldtube.InitializeCurrentIteration",
+    "[Unit][Evolution]") {
+  const size_t current_iteration = 77;
+  auto box =
+      db::create<db::AddSimpleTags<Tags::CurrentIteration>>(current_iteration);
+  db::mutate_apply<Initialization::InitializeCurrentIteration>(
+      make_not_null(&box));
+  CHECK(get<Tags::CurrentIteration>(box) == 0);
+}
+}  // namespace
+}  // namespace CurvedScalarWave::Worldtube

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

@@ -1,22 +1,17 @@
 // Distributed under the MIT License.
 // See LICENSE.txt for details.
-
-#include "Framework/TestingFramework.hpp"
-
 #include <array>
 #include <memory>
-
 #include "DataStructures/DataBox/DataBox.hpp"
 #include "DataStructures/Tensor/Tensor.hpp"
 #include "Evolution/Systems/CurvedScalarWave/Worldtube/SingletonActions/InitializeEvolvedVariables.hpp"
+#include "Framework/TestingFramework.hpp"
 #include "Time/Tags/HistoryEvolvedVariables.hpp"
 #include "Time/Tags/TimeStepper.hpp"
 #include "Time/TimeSteppers/AdamsBashforth.hpp"
 #include "Time/TimeSteppers/TimeStepper.hpp"
-
 namespace CurvedScalarWave::Worldtube {
 namespace {
-
 SPECTRE_TEST_CASE(
     "Unit.Evolution.Systems.CSW.Worldtube.InitializeEvolvedVariables",
     "[Unit][Evolution]") {
@@ -25,21 +20,23 @@ SPECTRE_TEST_CASE(
   using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>;
   const tnsr::I<double, 3> initial_pos{{1., 2., 3.}};
   const tnsr::I<double, 3> initial_vel{{4., 5., 6.}};
-  auto box = db::create<
-      db::AddSimpleTags<variables_tag, dt_variables_tag,
-                        ::Tags::HistoryEvolvedVariables<variables_tag>,
-                        ::Tags::ConcreteTimeStepper<TimeStepper>,
-                        Tags::InitialPositionAndVelocity>,
-      time_stepper_ref_tags<TimeStepper>>(
-      variables_tag::type{}, dt_variables_tag::type{},
-      TimeSteppers::History<variables_tag::type>{},
-      static_cast<std::unique_ptr<TimeStepper>>(
-          std::make_unique<TimeSteppers::AdamsBashforth>(4)),
-      std::array<tnsr::I<double, 3>, 2>{{initial_pos, initial_vel}});
+  const size_t current_iteration = 77;
+  auto box =
+      db::create<db::AddSimpleTags<
+                     variables_tag, dt_variables_tag,
+                     ::Tags::HistoryEvolvedVariables<variables_tag>,
+                     ::Tags::ConcreteTimeStepper<TimeStepper>,
+                     Tags::InitialPositionAndVelocity, Tags::CurrentIteration>,
+                 time_stepper_ref_tags<TimeStepper>>(
+          variables_tag::type{}, dt_variables_tag::type{},
+          TimeSteppers::History<variables_tag::type>{},
+          static_cast<std::unique_ptr<TimeStepper>>(
+              std::make_unique<TimeSteppers::AdamsBashforth>(4)),
+          std::array<tnsr::I<double, 3>, 2>{{initial_pos, initial_vel}},
+          current_iteration);
 
   db::mutate_apply<Initialization::InitializeEvolvedVariables>(
       make_not_null(&box));
-
   const auto vars = db::get<variables_tag>(box);
   for (size_t i = 0; i < 3; ++i) {
     CHECK(get<Tags::EvolvedPosition<3>>(vars).get(i)[0] == initial_pos.get(i));
@@ -49,6 +46,7 @@ SPECTRE_TEST_CASE(
         dt_variables_tag::type(size_t(1), 0.));
   CHECK(db::get<::Tags::HistoryEvolvedVariables<variables_tag>>(box) ==
         TimeSteppers::History<variables_tag::type>(1));
+  CHECK(get<Tags::CurrentIteration>(box) == 0);
 }
 }  // namespace
 }  // namespace CurvedScalarWave::Worldtube