Add UpdateFunctionsOfTime action

src/Domain/Creators/BinaryCompactObject.cpp

@@ -64,11 +64,13 @@ create_grid_anchors(const std::array<double, 3>& center_a,
 }
 }  // namespace bco
 
-bool BinaryCompactObject::Object::is_excised() const {
+template <bool UseWorldtube>
+bool BinaryCompactObject<UseWorldtube>::Object::is_excised() const {
   return inner_boundary_condition.has_value();
 }
 
-BinaryCompactObject::BinaryCompactObject(
+template <bool UseWorldtube>
+BinaryCompactObject<UseWorldtube>::BinaryCompactObject(
     typename ObjectA::type object_A, typename ObjectB::type object_B,
     std::array<double, 2> center_of_mass_offset, const double envelope_radius,
     const double outer_radius,
@@ -336,7 +338,8 @@ BinaryCompactObject::BinaryCompactObject(
   }
 }
 
-Domain<3> BinaryCompactObject::create_domain() const {
+template <bool UseWorldtube>
+Domain<3> BinaryCompactObject<UseWorldtube>::create_domain() const {
   const double inner_sphericity_A = is_excised_a_ ? 1.0 : 0.0;
   const double inner_sphericity_B = is_excised_b_ ? 1.0 : 0.0;
 
@@ -728,9 +731,10 @@ Domain<3> BinaryCompactObject::create_domain() const {
   return domain;
 }
 
+template <bool UseWorldtube>
 std::vector<DirectionMap<
     3, std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>>>
-BinaryCompactObject::external_boundary_conditions() const {
+BinaryCompactObject<UseWorldtube>::external_boundary_conditions() const {
   if (outer_boundary_condition_ == nullptr) {
     return {};
   }
@@ -766,16 +770,20 @@ BinaryCompactObject::external_boundary_conditions() const {
   return boundary_conditions;
 }
 
+template <bool UseWorldtube>
 std::unordered_map<std::string,
                    std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>
-BinaryCompactObject::functions_of_time(
+BinaryCompactObject<UseWorldtube>::functions_of_time(
     const std::unordered_map<std::string, double>& initial_expiration_times)
     const {
   return time_dependent_options_.has_value()
-             ? time_dependent_options_->create_functions_of_time(
+             ? time_dependent_options_->create_functions_of_time<UseWorldtube>(
                    initial_expiration_times)
              : std::unordered_map<
                    std::string,
                    std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>{};
 }
+
+template class BinaryCompactObject<true>;
+template class BinaryCompactObject<false>;
 }  // namespace domain::creators

src/Domain/Creators/BinaryCompactObject.hpp

@@ -144,7 +144,12 @@ create_grid_anchors(const std::array<double, 3>& center_a,
  * `domain::creators::bco::TimeDependentMapOptions`. This class must pass a
  * template parameter of `false` to
  * `domain::creators::bco::TimeDependentMapOptions`.
+ *
+ * The `UseWorldtube` template parameter is set to false by default. When set to
+ * true, some of the functions of time will be `IntegratedFunctionOfTime` used
+ * to control the orbit of the worldtube.
  */
+template <bool UseWorldtube = false>
 class BinaryCompactObject : public DomainCreator<3> {
  private:
   // Time-independent maps

src/Domain/Creators/Factory3D.hpp

@@ -21,8 +21,8 @@ template <>
 struct domain_creators<3> {
   using type =
       tmpl::list<domain::creators::AlignedLattice<3>,
-                 domain::creators::BinaryCompactObject, domain::creators::Brick,
-                 domain::creators::Cylinder,
+                 domain::creators::BinaryCompactObject<false>,
+                 domain::creators::Brick, domain::creators::Cylinder,
                  domain::creators::CylindricalBinaryCompactObject,
                  domain::creators::FrustalCloak,
                  domain::creators::RotatedBricks, domain::creators::Sphere>;

src/Domain/Creators/TimeDependentOptions/BinaryCompactObject.cpp

@@ -19,6 +19,7 @@
 #include "Domain/Creators/TimeDependentOptions/ShapeMap.hpp"
 #include "Domain/FunctionsOfTime/FixedSpeedCubic.hpp"
 #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
+#include "Domain/FunctionsOfTime/IntegratedFunctionOfTime.hpp"
 #include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp"
 #include "Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp"
 #include "NumericalAlgorithms/SphericalHarmonics/Spherepack.hpp"
@@ -71,11 +72,99 @@ TimeDependentMapOptions<IsCylindrical>::TimeDependentMapOptions(
 }
 
 template <bool IsCylindrical>
+std::unordered_map<std::string,
+                   std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>
+TimeDependentMapOptions<IsCylindrical>::create_worldtube_functions_of_time()
+    const {
+  if (translation_map_options_.has_value()) {
+    ERROR("Translation map is not implemented for worldtube evolutions.");
+  }
+  std::unordered_map<std::string,
+                     std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>
+      result{};
+  // The functions of time need to be valid only for the very first time step,
+  // after that they need to be updated by the worldtube singleton.
+  const double initial_expiration_time = initial_time_ + 1e-10;
+  if (not expansion_map_options_.has_value()) {
+    ERROR("Initial values for the expansion map need to be provided.");
+  }
+  result[expansion_name] =
+      std::make_unique<FunctionsOfTime::IntegratedFunctionOfTime>(
+          initial_time_,
+          std::array<double, 2>{
+              {{gsl::at(expansion_map_options_.value().initial_values, 0)},
+               {gsl::at(expansion_map_options_.value().initial_values, 1)}}},
+          initial_expiration_time, false);
+  result[expansion_outer_boundary_name] =
+      std::make_unique<FunctionsOfTime::FixedSpeedCubic>(
+          1.0, initial_time_,
+          expansion_map_options_.value().outer_boundary_velocity,
+          expansion_map_options_.value().outer_boundary_decay_time);
+  if (not rotation_map_options_.has_value()) {
+    ERROR(
+        "Initial values for the rotation map need to be provided when using "
+        "the worldtube.");
+  }
+
+  result[rotation_name] =
+      std::make_unique<FunctionsOfTime::IntegratedFunctionOfTime>(
+          initial_time_,
+          std::array<double, 2>{
+              0.,
+              gsl::at(rotation_map_options_.value().initial_angular_velocity,
+                      2)},
+          initial_expiration_time, true);
+
+  // Size and Shape FunctionOfTime for objects A and B. Only spherical excision
+  // spheres are supported currently.
+  if (not(shape_options_A_.has_value() and shape_options_B_.has_value())) {
+    ERROR(
+        "Initial size for both excision spheres need to be provided when using "
+        "the worldtube.");
+  }
+  for (size_t i = 0; i < shape_names.size(); i++) {
+    const auto make_initial_size_values = [](const auto& lambda_options) {
+      return std::array<double, 2>{
+          {gsl::at(lambda_options.value().initial_size_values.value(), 0),
+           gsl::at(lambda_options.value().initial_size_values.value(), 1)}};
+    };
+    const std::array<double, 2> initial_size_values =
+        i == 0 ? make_initial_size_values(shape_options_A_)
+               : make_initial_size_values(shape_options_B_);
+    const size_t initial_l_max = 2;
+    const DataVector shape_zeros{
+        ylm::Spherepack::spectral_size(initial_l_max, initial_l_max), 0.0};
+
+    result[gsl::at(shape_names, i)] =
+        std::make_unique<FunctionsOfTime::PiecewisePolynomial<2>>(
+            initial_time_,
+            std::array<DataVector, 3>{shape_zeros, shape_zeros, shape_zeros},
+            std::numeric_limits<double>::infinity());
+    result[gsl::at(size_names, i)] =
+        std::make_unique<FunctionsOfTime::IntegratedFunctionOfTime>(
+            initial_time_, initial_size_values, initial_expiration_time, false);
+  }
+  return result;
+}
+
+template <bool IsCylindrical>
+template <bool UseWorldtube>
 std::unordered_map<std::string,
                    std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>
 TimeDependentMapOptions<IsCylindrical>::create_functions_of_time(
     const std::unordered_map<std::string, double>& initial_expiration_times)
     const {
+  if constexpr (UseWorldtube) {
+    static_assert(not IsCylindrical,
+                  "Cylindrical map not supported with worldtube");
+    if (not initial_expiration_times.empty()) {
+      ERROR(
+          "Initial expiration times were specified with worldtube functions of "
+          "time. This is not supported, as the worldtube singleton has to set "
+          "the expiration times each time step");
+    }
+    return create_worldtube_functions_of_time();
+  }
   std::unordered_map<std::string,
                      std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>
       result{};
@@ -507,4 +596,17 @@ GENERATE_INSTANTIATIONS(INSTANTIATE, (true, false),
 #undef OBJECT
 #undef ISCYL
 #undef INSTANTIATE
+
+template std::unordered_map<
+    std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>
+TimeDependentMapOptions<false>::create_functions_of_time<true>(
+    const std::unordered_map<std::string, double>&) const;
+template std::unordered_map<
+    std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>
+TimeDependentMapOptions<false>::create_functions_of_time<false>(
+    const std::unordered_map<std::string, double>&) const;
+template std::unordered_map<
+    std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>
+TimeDependentMapOptions<true>::create_functions_of_time<false>(
+    const std::unordered_map<std::string, double>&) const;
 }  // namespace domain::creators::bco

src/Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp

@@ -277,7 +277,17 @@ struct TimeDependentMapOptions {
    * - Translation: `PiecewisePolynomial<3>`
    * - SizeA/B: `PiecewisePolynomial<3>`
    * - ShapeA/B: `PiecewisePolynomial<2>`
+   *
+   *  When `UseWorldtube` is set to true, they are
+   *
+   * - Expansion: `IntegratedFunctionOfTime`
+   * - ExpansionOuterBoundary: `FixedSpeedCubic`
+   * - Rotation: `IntegratedFunctionOfTime`
+   * - Translation: None
+   * - SizeA/B: IntegratedFunctionOfTime
+   * - ShapeA/B: `PiecewisePolynomial<2>`
    */
+  template <bool UseWorldtube = false>
   std::unordered_map<std::string,
                      std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>
   create_functions_of_time(const std::unordered_map<std::string, double>&
@@ -398,5 +408,10 @@ struct TimeDependentMapOptions {
       tmpl::conditional_t<IsCylindrical, std::array<std::optional<Shape>, 2>,
                           std::array<std::array<std::optional<Shape>, 6>, 2>>;
   ShapeMapType shape_maps_{};
+
+  // helper function that creates the functions of time used by the worldtube
+  std::unordered_map<std::string,
+                     std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>
+  create_worldtube_functions_of_time() const;
 };
 }  // namespace domain::creators::bco

src/Evolution/Executables/CurvedScalarWave/EvolveWorldtubeCurvedScalarWave.hpp

@@ -7,7 +7,7 @@
 #include <cstdint>
 #include <vector>
 
-#include "Domain/Creators/Factory3D.hpp"
+#include "Domain/Creators/BinaryCompactObject.hpp"
 #include "Domain/Creators/RegisterDerivedWithCharm.hpp"
 #include "Domain/Creators/TimeDependence/RegisterDerivedWithCharm.hpp"
 #include "Domain/ElementDistribution.hpp"
@@ -208,7 +208,8 @@ struct EvolutionMetavars {
                     standard_boundary_conditions<volume_dim>,
                 CurvedScalarWave::BoundaryConditions::Worldtube<volume_dim>>>>,
         tmpl::pair<DenseTrigger, DenseTriggers::standard_dense_triggers>,
-        tmpl::pair<DomainCreator<volume_dim>, domain_creators<volume_dim>>,
+        tmpl::pair<DomainCreator<volume_dim>,
+                   tmpl::list<domain::creators::BinaryCompactObject<false>>>,
         tmpl::pair<
             Event,
             tmpl::flatten<tmpl::list<

src/Evolution/Executables/GeneralizedHarmonic/EvolveGhBinaryBlackHole.hpp

@@ -451,7 +451,7 @@ struct EvolutionMetavars {
                        DenseTriggers::standard_dense_triggers>>>,
         tmpl::pair<
             DomainCreator<volume_dim>,
-            tmpl::list<::domain::creators::BinaryCompactObject,
+            tmpl::list<::domain::creators::BinaryCompactObject<false>,
                        ::domain::creators::CylindricalBinaryCompactObject>>,
         tmpl::pair<
             Event,

src/Evolution/Executables/GrMhd/GhValenciaDivClean/GhValenciaDivCleanBase.hpp

@@ -586,13 +586,14 @@ struct GhValenciaDivCleanTemplateBase<
                    tmpl::append<DenseTriggers::standard_dense_triggers,
                                 control_system::control_system_triggers<
                                     control_systems>>>,
-        tmpl::pair<DomainCreator<volume_dim>,
-                   // Currently control systems can only be used with BCO
-                   // domains
-                   tmpl::conditional_t<
-                       use_control_systems,
-                       tmpl::list<::domain::creators::BinaryCompactObject>,
-                       domain_creators<volume_dim>>>,
+        tmpl::pair<
+            DomainCreator<volume_dim>,
+            // Currently control systems can only be used with BCO
+            // domains
+            tmpl::conditional_t<
+                use_control_systems,
+                tmpl::list<::domain::creators::BinaryCompactObject<false>>,
+                domain_creators<volume_dim>>>,
         tmpl::pair<Event,
                    tmpl::flatten<tmpl::list<
                        Events::Completion,

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

@@ -15,6 +15,7 @@ spectre_target_sources(
   PunctureField.cpp
   PunctureFieldOrder0.cpp
   PunctureFieldOrder1.cpp
+  RadiusFunctions.cpp
   SelfForce.cpp
   )
 
@@ -27,13 +28,15 @@ spectre_target_headers(
   SingletonChare.hpp
   Tags.hpp
   PunctureField.hpp
+  RadiusFunctions.hpp
   SelfForce.hpp
   Worldtube.hpp
   )
 
 target_link_libraries(
   ${LIBRARY}
   PUBLIC
+  ControlSystem
   DataStructures
   Domain
   GeneralRelativity

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

@@ -0,0 +1,57 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Evolution/Systems/CurvedScalarWave/Worldtube/RadiusFunctions.hpp"
+
+#include <cmath>
+#include <cstddef>
+
+#include "Utilities/ErrorHandling/Error.hpp"
+
+namespace CurvedScalarWave::Worldtube {
+
+double smooth_broken_power_law(double orbit_radius, double alpha, double amp,
+                               double rb, double delta) {
+  if (alpha < 0. or alpha > 4.) {
+    ERROR(
+        "Only exponents between 0 and 4 have been tested for the broken law. "
+        "Values outside this range are likely not sensible.");
+  }
+  if (delta < 1e-3) {
+    ERROR(
+        "A value of delta less than 0.001 is disabled as it can lead to "
+        "prohibitively large floating precision errors.");
+  }
+  const double r_by_rb = orbit_radius / rb;
+  return amp * pow(r_by_rb, alpha) *
+         pow(0.5 * (1. + pow(r_by_rb, 1. / delta)), -alpha * delta);
+}
+
+double smooth_broken_power_law_derivative(double orbit_radius, double alpha,
+                                          double amp, double rb, double delta) {
+  if (alpha < 0. or alpha > 4.) {
+    ERROR(
+        "Only exponents between 0 and 4 have been tested for the broken law. "
+        "Values outside this range are likely not sensible.");
+  }
+  if (delta < 1e-3) {
+    ERROR(
+        "A value of delta less than 0.001 is disabled as it can lead to "
+        "prohibitively large floating precision errors.");
+  }
+
+  // During testing it was found that for large values of rb, the equation below
+  // can FPE even though the derivative is clearly zero up to double precision.
+  // To avoid this, we introduce this shortcut.
+  if (orbit_radius > rb + 1000. * delta) {
+    return 0.;
+  }
+  const double r_by_rb = orbit_radius / rb;
+  const double temp1 = pow(r_by_rb, 1. / delta - 1.);
+  const double temp2 = 0.5 * temp1 * r_by_rb + 0.5;
+
+  return amp * alpha * pow(r_by_rb, alpha - 1) / rb *
+         pow(temp2, -alpha * delta - 1.) * (temp2 - 0.5 * temp1 * r_by_rb);
+}
+
+}  // namespace CurvedScalarWave::Worldtube