Enable time dependent translation with transition for Sphere domain

src/Domain/CoordinateMaps/TimeDependent/Translation.cpp

@@ -359,7 +359,8 @@ std::array<tt::remove_cvref_wrap_t<T>, Dim> Translation<Dim>::piecewise_helper(
       for (size_t i = 0; i < Dim; i++) {
         get_element(gsl::at(result, i), k) += gsl::at(func_or_deriv_of_time, i);
       }
-      ASSERT(max(magnitude(result)) <= outer_radius_.value(),
+      const double eps = std::numeric_limits<double>::epsilon() * 100.0;
+      ASSERT(max(magnitude(result)) - eps <= outer_radius_.value(),
              "Coordinates translated outside outer radius, This should not "
              "happen");
     } else if (get_element(radius, k) > inner_radius_.value() and
@@ -373,7 +374,8 @@ std::array<tt::remove_cvref_wrap_t<T>, Dim> Translation<Dim>::piecewise_helper(
         get_element(gsl::at(result, i), k) +=
             gsl::at(func_or_deriv_of_time, i) * radial_falloff_factor;
       }
-      ASSERT(max(magnitude(result)) <= outer_radius_.value(),
+      const double eps = std::numeric_limits<double>::epsilon() * 100.0;
+      ASSERT(max(magnitude(result)) - eps <= outer_radius_.value(),
              "Coordinates translated outside outer radius, This should not "
              "happen");
     }

src/Domain/Creators/Sphere.cpp

@@ -267,11 +267,23 @@ 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. Require at least
+      // one radial partition, as the transition requires spherical shape.
+
+      if (radial_partitioning_.empty()) {
+        PARSE_ERROR(context,
+                    "The hard-coded translation map requires at least two "
+                    "shells. Use at least one radial partition.");
+      }
+
       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
+                      std::pair<double, double>{radial_partitioning_.back(),
+                                                outer_radius_});
     }
   }
 }
@@ -374,9 +386,14 @@ Domain<3> Sphere::create_domain() const {
               time_dependent_options_.value());
 
       // First shell gets the distorted maps.
+      // Last shell gets translation with transition region
       for (size_t block_id = 0; block_id < num_blocks_; block_id++) {
         const bool include_distorted_map_in_first_shell =
             block_id < num_blocks_per_shell_;
+        // False if block_id is in the last shell
+        const bool use_rigid_translation =
+            block_id + num_blocks_per_shell_ + (fill_interior_ ? 1 : 0) <
+            num_blocks_;
         block_maps_grid_to_distorted[block_id] =
             hard_coded_options.grid_to_distorted_map(
                 include_distorted_map_in_first_shell);
@@ -385,7 +402,7 @@ Domain<3> Sphere::create_domain() const {
                 include_distorted_map_in_first_shell);
         block_maps_grid_to_inertial[block_id] =
             hard_coded_options.grid_to_inertial_map(
-                include_distorted_map_in_first_shell);
+                include_distorted_map_in_first_shell, use_rigid_translation);
       }
     } else {
       const auto& time_dependence = std::get<std::unique_ptr<

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 TranslationMapOptions& translation_map_options)
     : 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_translation_values_(translation_map_options.initial_values) {}
 
 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,32 @@ TimeDependentMapOptions::create_functions_of_time(
                                  {0.0}}},
       expiration_times.at(size_name));
 
+  DataVector initial_translation_center_temp{3, 0.0};
+  DataVector initial_translation_velocity_temp{3, 0.0};
+  for (size_t i = 0; i < 3; i++) {
+    initial_translation_center_temp[i] =
+        gsl::at(initial_translation_values_.front(), i);
+    initial_translation_velocity_temp[i] =
+        gsl::at(initial_translation_values_.back(), i);
+  }
+
+  // TranslationMap FunctionOfTime
+  result[translation_name] =
+      std::make_unique<FunctionsOfTime::PiecewisePolynomial<2>>(
+          initial_time_,
+          std::array<DataVector, 3>{
+              {std::move(initial_translation_center_temp),
+               std::move(initial_translation_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,
+    std::pair<double, double> translation_transition_radii) {
   std::unique_ptr<domain::CoordinateMaps::ShapeMapTransitionFunctions::
                       ShapeMapTransitionFunction>
       transition_func =
@@ -115,17 +139,22 @@ 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};
+
+  rigid_translation_map_ = TranslationMap{translation_name};
+
+  linear_falloff_translation_map_ =
+      TranslationMap{translation_name, translation_transition_radii.first,
+                     translation_transition_radii.second};
 }
 
 // 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>(
+        rigid_translation_map_);
   } else {
     return nullptr;
   }
@@ -143,12 +172,15 @@ TimeDependentMapOptions::grid_to_distorted_map(
 
 TimeDependentMapOptions::MapType<Frame::Grid, Frame::Inertial>
 TimeDependentMapOptions::grid_to_inertial_map(
-    const bool include_distorted_map) const {
+    const bool include_distorted_map, const bool use_rigid_translation) const {
   if (include_distorted_map) {
-    return std::make_unique<GridToInertialComposition>(shape_map_);
+    return std::make_unique<GridToInertialComposition>(shape_map_,
+                                                       rigid_translation_map_);
+  } else if (use_rigid_translation) {
+    return std::make_unique<GridToInertialSimple>(rigid_translation_map_);
   } else {
-    return std::make_unique<
-        IdentityForComposition<Frame::Grid, Frame::Inertial>>(IdentityMap{});
+    return std::make_unique<GridToInertialSimple>(
+        linear_falloff_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 {
@@ -131,15 +139,35 @@ struct TimeDependentMapOptions {
     std::optional<std::variant<KerrSchildFromBoyerLindquist>> initial_values{};
   };
 
-  using options = tmpl::list<InitialTime, ShapeMapOptions>;
+  struct TranslationMapOptions {
+    using type = TranslationMapOptions;
+    static std::string name() { return "TranslationMap"; }
+    static constexpr Options::String help = {
+        "Options for a time-dependent translation map in that keeps the outer "
+        "boundary fixed."};
+
+    struct InitialValues {
+      using type = std::array<std::array<double, 3>, 2>;
+      static constexpr Options::String help = {
+          "Initial values for the translation map."};
+    };
+
+    using options = tmpl::list<InitialValues>;
+
+    std::array<std::array<double, 3>, 2> initial_values{};
+  };
+
+  using options =
+      tmpl::list<InitialTime, ShapeMapOptions, TranslationMapOptions>;
   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 TranslationMapOptions& translation_map_options);
 
   /*!
    * \brief Create the function of time map using the options that were
@@ -149,6 +177,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 +191,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, double inner_radius,
+      double outer_radius,
+      std::pair<double, double> translation_transition_radii);
 
   /*!
    * \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
@@ -195,16 +227,20 @@ struct TimeDependentMapOptions {
    * only be an identity map.
    */
   MapType<Frame::Grid, Frame::Inertial> grid_to_inertial_map(
-      bool include_distorted_map) const;
+      bool include_distorted_map, bool use_rigid_translation) const;
 
   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 rigid_translation_map_{};
+  TranslationMap linear_falloff_translation_map_{};
   std::optional<std::variant<KerrSchildFromBoyerLindquist>>
       initial_shape_values_{};
+  std::array<std::array<double, 3>, 2> initial_translation_values_{};
 };
 }  // namespace domain::creators::sphere

tests/InputFiles/Xcts/KerrSchild.yaml

@@ -4,7 +4,7 @@
 Executable: SolveXcts
 Testing:
   Check: parse;execute_check_output
-  Timeout: 40
+  Timeout: 50
   Priority: High
 ExpectedOutput:
   - KerrSchildReductions.h5
@@ -26,7 +26,7 @@ Background: &solution
   KerrSchild:
     Mass: &KerrMass 1.
     Spin: &KerrSpin [0., 0., 0.6]
-    Center: [0., 0., 0.]
+    Center: &Center [0., 0., 0.]
 
 InitialGuess: Flatness
 
@@ -48,15 +48,17 @@ DomainCreator:
     UseEquiangularMap: True
     EquatorialCompression: None
     WhichWedges: All
-    RadialPartitioning: []
-    RadialDistribution: [Logarithmic]
+    RadialPartitioning: [3.0]
+    RadialDistribution: [Linear, Logarithmic]
     TimeDependentMaps:
       InitialTime: 0.
       ShapeMap:
         LMax: 20
         InitialValues:
           Mass: *KerrMass
           Spin: *KerrSpin
+      TranslationMap:
+        InitialValues: [*Center, [0., 0., 0.]]
     OuterBoundaryCondition:
       AnalyticSolution:
         Solution: *solution