Add translation to sphere domain time maps

src/Domain/Creators/Sphere.cpp

@@ -267,11 +267,21 @@ Sphere::Sphere(
       // Build the maps. We only apply the maps in the inner-most shell. The
       // inner radius is what's passed in, but the outer radius is the outer
       // radius of the inner-most shell so we have to see how many shells we
-      // have
+      // have.
+      // The translation map linear transition occurs only in the outer-most
+      // shell, such that the outer boundary is not translated. If there
+      // is no radial partition, a uniform translation is applied instead, as
+      // the transition requires spherical shape.
       std::get<sphere::TimeDependentMapOptions>(time_dependent_options_.value())
           .build_maps(std::array{0.0, 0.0, 0.0}, inner_radius_,
                       radial_partitioning_.empty() ? outer_radius_
-                                                   : radial_partitioning_[0]);
+                                                   : radial_partitioning_[0],
+                      // Translation inner radius
+                      radial_partitioning_.empty()
+                          ? outer_radius_
+                          : radial_partitioning_.back(),
+                      // Translation outer radius
+                      outer_radius_);
     }
   }
 }

src/Domain/Creators/SphereTimeDependentMaps.cpp

@@ -26,11 +26,14 @@
 #include "Utilities/ErrorHandling/Error.hpp"
 
 namespace domain::creators::sphere {
+
 TimeDependentMapOptions::TimeDependentMapOptions(
-    const double initial_time, const ShapeMapOptions& shape_map_options)
+    const double initial_time, const ShapeMapOptions& shape_map_options,
+    const std::array<double, 3>& initial_velocity)
     : initial_time_(initial_time),
       initial_l_max_(shape_map_options.l_max),
-      initial_shape_values_(shape_map_options.initial_values) {}
+      initial_shape_values_(shape_map_options.initial_values),
+      initial_velocity_(initial_velocity) {}
 
 std::unordered_map<std::string,
                    std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>
@@ -46,7 +49,8 @@ TimeDependentMapOptions::create_functions_of_time(
   // their initial expiration time to infinity (i.e. not expiring)
   std::unordered_map<std::string, double> expiration_times{
       {size_name, std::numeric_limits<double>::infinity()},
-      {shape_name, std::numeric_limits<double>::infinity()}};
+      {shape_name, std::numeric_limits<double>::infinity()},
+      {translation_name, std::numeric_limits<double>::infinity()}};
 
   // If we have control systems, overwrite these expiration times with the ones
   // supplied by the control system
@@ -101,12 +105,26 @@ TimeDependentMapOptions::create_functions_of_time(
                                  {0.0}}},
       expiration_times.at(size_name));
 
+  DataVector initial_velocity_temp{3, 0.0};
+  for (size_t i = 0; i < 3; i++) {
+    initial_velocity_temp[i] = gsl::at(initial_velocity_, i);
+  }
+
+  // TranslationMap FunctionOfTime
+  result[translation_name] =
+      std::make_unique<FunctionsOfTime::PiecewisePolynomial<2>>(
+          initial_time_,
+          std::array<DataVector, 3>{
+              {{3, 0.0}, initial_velocity_temp, {3, 0.0}}},
+          expiration_times.at(translation_name));
+
   return result;
 }
 
-void TimeDependentMapOptions::build_maps(const std::array<double, 3>& center,
-                                         const double inner_radius,
-                                         const double outer_radius) {
+void TimeDependentMapOptions::build_maps(
+    const std::array<double, 3>& center, const double inner_radius,
+    const double outer_radius, const double translation_inner_radius,
+    const double translation_outer_radius) {
   std::unique_ptr<domain::CoordinateMaps::ShapeMapTransitionFunctions::
                       ShapeMapTransitionFunction>
       transition_func =
@@ -115,17 +133,24 @@ void TimeDependentMapOptions::build_maps(const std::array<double, 3>& center,
   shape_map_ = ShapeMap{center,         initial_l_max_,
                         initial_l_max_, std::move(transition_func),
                         shape_name,     size_name};
+
+  if (translation_inner_radius == translation_outer_radius) {
+    // Uniform translation
+    translation_map_ = TranslationMap{translation_name};
+  } else {
+    // Translation with transition
+    translation_map_ = TranslationMap{
+        translation_name, translation_inner_radius, translation_outer_radius};
+  }
 }
 
 // If you edit any of the functions below, be sure to update the documentation
 // in the Sphere domain creator as well as this class' documentation.
 TimeDependentMapOptions::MapType<Frame::Distorted, Frame::Inertial>
 TimeDependentMapOptions::distorted_to_inertial_map(
-    const bool include_distorted_map) {
+    const bool include_distorted_map) const {
   if (include_distorted_map) {
-    return std::make_unique<
-        IdentityForComposition<Frame::Distorted, Frame::Inertial>>(
-        IdentityMap{});
+    return std::make_unique<DistortedToInertialComposition>(translation_map_);
   } else {
     return nullptr;
   }
@@ -145,10 +170,10 @@ TimeDependentMapOptions::MapType<Frame::Grid, Frame::Inertial>
 TimeDependentMapOptions::grid_to_inertial_map(
     const bool include_distorted_map) const {
   if (include_distorted_map) {
-    return std::make_unique<GridToInertialComposition>(shape_map_);
+    return std::make_unique<GridToInertialComposition>(shape_map_,
+                                                       translation_map_);
   } else {
-    return std::make_unique<
-        IdentityForComposition<Frame::Grid, Frame::Inertial>>(IdentityMap{});
+    return std::make_unique<GridToInertialSimple>(translation_map_);
   }
 }
 }  // namespace domain::creators::sphere

src/Domain/Creators/SphereTimeDependentMaps.hpp

@@ -14,6 +14,7 @@
 #include "Domain/CoordinateMaps/CoordinateMap.hpp"
 #include "Domain/CoordinateMaps/Identity.hpp"
 #include "Domain/CoordinateMaps/TimeDependent/Shape.hpp"
+#include "Domain/CoordinateMaps/TimeDependent/Translation.hpp"
 #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
 #include "Options/Auto.hpp"
 #include "Options/String.hpp"
@@ -81,21 +82,28 @@ struct TimeDependentMapOptions {
   using IdentityMap = domain::CoordinateMaps::Identity<3>;
   // Time-dependent maps
   using ShapeMap = domain::CoordinateMaps::TimeDependent::Shape;
+  using TranslationMap = domain::CoordinateMaps::TimeDependent::Translation<3>;
 
   template <typename SourceFrame, typename TargetFrame>
   using IdentityForComposition =
       domain::CoordinateMap<SourceFrame, TargetFrame, IdentityMap>;
   using GridToDistortedComposition =
       domain::CoordinateMap<Frame::Grid, Frame::Distorted, ShapeMap>;
   using GridToInertialComposition =
-      domain::CoordinateMap<Frame::Grid, Frame::Inertial, ShapeMap>;
+      domain::CoordinateMap<Frame::Grid, Frame::Inertial, ShapeMap,
+                            TranslationMap>;
+  using GridToInertialSimple =
+      domain::CoordinateMap<Frame::Grid, Frame::Inertial, TranslationMap>;
+  using DistortedToInertialComposition =
+      domain::CoordinateMap<Frame::Distorted, Frame::Inertial, TranslationMap>;
 
  public:
   using maps_list =
       tmpl::list<IdentityForComposition<Frame::Grid, Frame::Inertial>,
                  IdentityForComposition<Frame::Grid, Frame::Distorted>,
                  IdentityForComposition<Frame::Distorted, Frame::Inertial>,
-                 GridToDistortedComposition, GridToInertialComposition>;
+                 GridToDistortedComposition, GridToInertialSimple,
+                 GridToInertialComposition, DistortedToInertialComposition>;
 
   /// \brief The initial time of the functions of time.
   struct InitialTime {
@@ -104,6 +112,13 @@ struct TimeDependentMapOptions {
         "The initial time of the functions of time"};
   };
 
+  struct TranslationMapOptions {
+    static std::string name() { return "TranslationMap"; }
+    static constexpr Options::String help = {
+        "Options for a time-dependent translation map in the inner-most shell "
+        "of the domain."};
+  };
+
   struct ShapeMapOptions {
     using type = ShapeMapOptions;
     static std::string name() { return "ShapeMap"; }
@@ -131,15 +146,24 @@ struct TimeDependentMapOptions {
     std::optional<std::variant<KerrSchildFromBoyerLindquist>> initial_values{};
   };
 
-  using options = tmpl::list<InitialTime, ShapeMapOptions>;
+  struct InitialVelocity {
+    static std::string name() { return "InitialVelocity"; }
+    using type = std::array<double, 3>;
+    static constexpr Options::String help = {
+        "Velocity of the translation map."};
+    using group = TranslationMapOptions;
+  };
+
+  using options = tmpl::list<InitialTime, ShapeMapOptions, InitialVelocity>;
   static constexpr Options::String help{
       "The options for all the hard-coded time dependent maps in the Sphere "
       "domain."};
 
   TimeDependentMapOptions() = default;
 
   TimeDependentMapOptions(double initial_time,
-                          const ShapeMapOptions& shape_map_options);
+                          const ShapeMapOptions& shape_map_options,
+                          const std::array<double, 3>& initial_velocity);
 
   /*!
    * \brief Create the function of time map using the options that were
@@ -149,6 +173,7 @@ struct TimeDependentMapOptions {
    *
    * - Size: `PiecewisePolynomial<3>`
    * - Shape: `PiecewisePolynomial<2>`
+   * - Translation: `PiecewisePolynomial<2>`
    */
   std::unordered_map<std::string,
                      std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>
@@ -162,19 +187,22 @@ struct TimeDependentMapOptions {
    * Currently, this constructs a:
    *
    * - Shape: `Shape` (with a size function of time)
+   * - Translation: `Translation`
    */
-  void build_maps(const std::array<double, 3>& center, double inner_radius,
-                  double outer_radius);
+  void build_maps(const std::array<double, 3>& center,
+                  const double inner_radius, const double outer_radius,
+                  const double translation_inner_radius,
+                  const double translation_outer_radius);
 
   /*!
    * \brief This will construct the map from `Frame::Distorted` to
    * `Frame::Inertial`.
    *
-   * If the argument `include_distorted_map` is true, then this will be an
-   * identity map. If it is false, then this returns `nullptr`.
+   * If the argument `include_distorted_map` is true, then this will be a
+   * translation map. If it is false, then this returns `nullptr`.
    */
-  static MapType<Frame::Distorted, Frame::Inertial> distorted_to_inertial_map(
-      bool include_distorted_map);
+  MapType<Frame::Distorted, Frame::Inertial> distorted_to_inertial_map(
+      bool include_distorted_map) const;
 
   /*!
    * \brief This will construct the map from `Frame::Grid` to
@@ -199,12 +227,16 @@ struct TimeDependentMapOptions {
 
   inline static const std::string size_name{"Size"};
   inline static const std::string shape_name{"Shape"};
+  inline static const std::string translation_name{"Translation"};
 
  private:
   double initial_time_{std::numeric_limits<double>::signaling_NaN()};
   size_t initial_l_max_{};
   ShapeMap shape_map_{};
+  TranslationMap translation_map_{};
   std::optional<std::variant<KerrSchildFromBoyerLindquist>>
       initial_shape_values_{};
+  std::array<double, 3> initial_velocity_{
+      std::numeric_limits<std::array<double, 3>>::signaling_NaN()};
 };
 }  // namespace domain::creators::sphere

tests/Unit/Domain/Creators/Test_Sphere.cpp

@@ -139,6 +139,8 @@ std::string option_string(
                               "    ShapeMap:\n"
                               "      LMax: 10\n"
                               "      InitialValues: Spherical\n"
+                              "    TranslationMap:\n"
+                              "      InitialVelocity: [-0.3, 0.5, -0.7]\n"
                             : "  TimeDependentMaps:\n"
                               "    UniformTranslation:\n"
                               "      InitialTime: 1.0\n"
@@ -535,11 +537,14 @@ void test_parse_errors() {
           "an outflow-type boundary condition, you must use that."));
   using TDMO = domain::creators::sphere::TimeDependentMapOptions;
   CHECK_THROWS_WITH(
-      creators::Sphere(
-          inner_radius, outer_radius, inner_cube, refinement, initial_extents,
-          use_equiangular_map, equatorial_compression, radial_partitioning,
-          radial_distribution, which_wedges,
-          TDMO{1.0, TDMO::ShapeMapOptions{5, std::nullopt}}, nullptr),
+      creators::Sphere(inner_radius, outer_radius, inner_cube, refinement,
+                       initial_extents, use_equiangular_map,
+                       equatorial_compression, radial_partitioning,
+                       radial_distribution, which_wedges,
+                       TDMO{1.0, TDMO::ShapeMapOptions{5, std::nullopt},
+                            // Translation options: InitialVelocity
+                            std::array<double, 3>{-0.3, 0.5, -0.7}},
+                       nullptr),
       Catch::Matchers::ContainsSubstring(
           "Currently cannot use hard-coded time dependent maps with an inner "
           "cube. Use a TimeDependence instead."));
@@ -630,8 +635,11 @@ void test_sphere() {
     if (time_dependent) {
       if (use_hard_coded_time_dep_options) {
         time_dependent_options = creators::sphere::TimeDependentMapOptions(
-            1.0, creators::sphere::TimeDependentMapOptions::ShapeMapOptions{
-                     l_max, std::nullopt});
+            1.0,
+            creators::sphere::TimeDependentMapOptions::ShapeMapOptions{
+                l_max, std::nullopt},
+            // Translation options: InitialVelocity
+            std::array<double, 3>{-0.3, 0.5, -0.7});
       } else {
         time_dependent_options = std::make_unique<
             domain::creators::time_dependence::UniformTranslation<3, 0>>(

tests/Unit/Domain/Creators/Test_SphereTimeDependentMaps.cpp

@@ -34,7 +34,9 @@ std::string create_option_string(const bool use_non_zero_shape) {
          (use_non_zero_shape ? "  InitialValues:\n"
                                "    Mass: 1.0\n"
                                "    Spin: [0.0, 0.0, 0.0]\n"s
-                             : "  InitialValues: Spherical\n"s);
+                             : "  InitialValues: Spherical\n"s) +
+         "TranslationMap:\n"
+         "  InitialVelocity: [-0.3, 0.5, -0.7]\n"s;
 }
 void test(const bool use_non_zero_shape) {
   CAPTURE(use_non_zero_shape);
@@ -48,11 +50,14 @@ void test(const bool use_non_zero_shape) {
   std::unordered_map<std::string, double> expiration_times{
       {TimeDependentMapOptions::shape_name, 15.5},
       {TimeDependentMapOptions::size_name,
+       std::numeric_limits<double>::infinity()},
+      {TimeDependentMapOptions::translation_name,
        std::numeric_limits<double>::infinity()}};
 
   // These are hard-coded so this is just a regression test
   CHECK(TimeDependentMapOptions::size_name == "Size"s);
   CHECK(TimeDependentMapOptions::shape_name == "Shape"s);
+  CHECK(TimeDependentMapOptions::translation_name == "Translation"s);
 
   using PP2 = domain::FunctionsOfTime::PiecewisePolynomial<2>;
   using PP3 = domain::FunctionsOfTime::PiecewisePolynomial<3>;
@@ -70,9 +75,17 @@ void test(const bool use_non_zero_shape) {
       std::array<DataVector, 3>{shape_zeros, shape_zeros, shape_zeros},
       expiration_times.at(TimeDependentMapOptions::shape_name)};
 
+  DataVector velocity_temp{{-0.3, 0.5, -0.7}};
+  PP2 translation_non_zero{
+      initial_time,
+      std::array<DataVector, 3>{{{3, 0.0}, velocity_temp, {3, 0.0}}},
+      expiration_times.at(TimeDependentMapOptions::translation_name)};
+
   const std::array<double, 3> center{-5.0, -0.01, -0.02};
   const double inner_radius = 0.5;
   const double outer_radius = 2.1;
+  const double translation_inner_radius = 3.1;
+  const double translation_outer_radius = 3.9;
 
   const auto functions_of_time =
       time_dep_options.create_functions_of_time(inner_radius, expiration_times);
@@ -86,8 +99,14 @@ void test(const bool use_non_zero_shape) {
             *functions_of_time.at(TimeDependentMapOptions::shape_name).get()) ==
         shape_all_zero);
 
+  CHECK(dynamic_cast<PP2&>(
+            *functions_of_time.at(TimeDependentMapOptions::translation_name)
+                 .get()) == translation_non_zero);
+
   for (const bool include_distorted : make_array(true, false)) {
-    time_dep_options.build_maps(center, inner_radius, outer_radius);
+    time_dep_options.build_maps(center, inner_radius, outer_radius,
+                                translation_inner_radius,
+                                translation_outer_radius);
 
     const auto grid_to_distorted_map =
         time_dep_options.grid_to_distorted_map(include_distorted);
@@ -112,9 +131,13 @@ void test(const bool use_non_zero_shape) {
       }
     };
 
-    check_map(grid_to_inertial_map, false, not include_distorted);
+    // There is no null pointer in the grid to inertial map
+    check_map(grid_to_inertial_map, false, false);
+
     check_map(grid_to_distorted_map, not include_distorted, false);
-    check_map(distorted_to_inertial_map, not include_distorted, true);
+
+    // If no shape distortion, there is only the translation
+    check_map(distorted_to_inertial_map, not include_distorted, false);
   }
 }