BinaryCompactObject domain: support initial Kerr horizon shape

src/Domain/Creators/BinaryCompactObject.cpp

@@ -564,12 +564,14 @@ Domain<3> BinaryCompactObject::create_domain() const {
     // Inside the excision sphere we add the grid to inertial map from the outer
     // shell. This allows the center of the excisions/horizons to be mapped
     // properly to the inertial frame.
-    if (is_excised_a_) {
+    if (is_excised_a_ and
+        grid_to_inertial_block_maps[number_of_blocks_ - 1] != nullptr) {
       domain.inject_time_dependent_map_for_excision_sphere(
           "ExcisionSphereA",
           grid_to_inertial_block_maps[number_of_blocks_ - 1]->get_clone());
     }
-    if (is_excised_b_) {
+    if (is_excised_b_ and
+        grid_to_inertial_block_maps[number_of_blocks_ - 1] != nullptr) {
       domain.inject_time_dependent_map_for_excision_sphere(
           "ExcisionSphereB",
           grid_to_inertial_block_maps[number_of_blocks_ - 1]->get_clone());
@@ -620,14 +622,18 @@ Domain<3> BinaryCompactObject::create_domain() const {
             time_dependent_options_
                 ->distorted_to_inertial_map<domain::ObjectLabel::B>(
                     block_for_distorted_frame);
-      } else {
+      } else if (grid_to_inertial_block_maps[number_of_blocks_ - 1] !=
+                 nullptr) {
         // No distorted frame
         grid_to_inertial_block_maps[block] =
             grid_to_inertial_block_maps[number_of_blocks_ - 1]->get_clone();
       }
     }
     // Finally, inject the time dependent maps into the corresponding blocks
     for (size_t block = 0; block < number_of_blocks_; ++block) {
+      if (grid_to_inertial_block_maps[block] == nullptr) {
+        continue;
+      }
       domain.inject_time_dependent_map_for_block(
           block, std::move(grid_to_inertial_block_maps[block]),
           std::move(grid_to_distorted_block_maps[block]),
@@ -682,6 +688,14 @@ BinaryCompactObject::functions_of_time(
     const {
   return time_dependent_options_.has_value()
              ? time_dependent_options_->create_functions_of_time(
+                   std::holds_alternative<Object>(object_A_)
+                       ? std::make_optional(
+                             std::get<Object>(object_A_).inner_radius)
+                       : std::nullopt,
+                   std::holds_alternative<Object>(object_B_)
+                       ? std::make_optional(
+                             std::get<Object>(object_B_).inner_radius)
+                       : std::nullopt,
                    initial_expiration_times)
              : std::unordered_map<
                    std::string,

src/Domain/Creators/BinaryCompactObjectHelpers.cpp

@@ -22,6 +22,7 @@
 #include "Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp"
 #include "NumericalAlgorithms/SphericalHarmonics/Spherepack.hpp"
 #include "Options/ParseError.hpp"
+#include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrHorizon.hpp"
 #include "Utilities/ErrorHandling/Assert.hpp"
 #include "Utilities/ErrorHandling/Error.hpp"
 #include "Utilities/GenerateInstantiations.hpp"
@@ -80,6 +81,8 @@ template <bool IsCylindrical>
 std::unordered_map<std::string,
                    std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>
 TimeDependentMapOptions<IsCylindrical>::create_functions_of_time(
+    const std::optional<double> inner_radius_A,
+    const std::optional<double> inner_radius_B,
     const std::unordered_map<std::string, double>& initial_expiration_times)
     const {
   std::unordered_map<std::string,
@@ -145,34 +148,75 @@ TimeDependentMapOptions<IsCylindrical>::create_functions_of_time(
   }
 
   // Size and Shape FunctionOfTime for objects A and B
-  for (size_t i = 0; i < shape_names.size(); i++) {
-    if (i == 0 ? shape_options_A_.has_value() : shape_options_B_.has_value()) {
-      const auto make_initial_size_values = [](const auto& lambda_options) {
-        return std::array<DataVector, 4>{
-            {{gsl::at(lambda_options.value().initial_size_values, 0)},
-             {gsl::at(lambda_options.value().initial_size_values, 1)},
-             {gsl::at(lambda_options.value().initial_size_values, 2)},
-             {0.0}}};
-      };
+  const auto build_shape_and_size_fot = [&result, &expiration_times, this](
+                                            const auto& shape_options,
+                                            const double inner_radius,
+                                            const std::string& shape_name,
+                                            const std::string& size_name) {
+    const DataVector shape_zeros{ylm::Spherepack::spectral_size(
+                                     shape_options.l_max, shape_options.l_max),
+                                 0.0};
+    DataVector shape_func{};
+    DataVector size_func{1, shape_options.initial_size_values[0]};
+    if (shape_options.initial_values.has_value()) {
+      if (std::holds_alternative<sphere::KerrSchildFromBoyerLindquist>(
+              shape_options.initial_values.value())) {
+        const ylm::Spherepack ylm{shape_options.l_max, shape_options.l_max};
+        const auto& mass_and_spin =
+            std::get<sphere::KerrSchildFromBoyerLindquist>(
+                shape_options.initial_values.value());
+        const DataVector radial_distortion =
+            1.0 - get(gr::Solutions::kerr_schild_radius_from_boyer_lindquist(
+                      inner_radius, ylm.theta_phi_points(), mass_and_spin.mass,
+                      mass_and_spin.spin)) /
+                      inner_radius;
+        shape_func = ylm.phys_to_spec(radial_distortion);
+        // Transform from SPHEREPACK to actual Ylm for size func
+        if (size_func[0] != 0.0) {
+          ERROR(
+              "Initial value for size map must be zero, because it is "
+              "overridden by the initial shape map values.");
+        }
+        size_func[0] = shape_func[0] * sqrt(0.5 * M_PI);
+        // Set l=0 for shape map to 0 because size is going to be used
+        shape_func[0] = 0.0;
+      }
+    } else {
+      shape_func = shape_zeros;
+    }
 
-      const size_t initial_l_max = i == 0 ? shape_options_A_.value().l_max
-                                          : shape_options_B_.value().l_max;
-      const std::array<DataVector, 4> initial_size_values =
-          i == 0 ? make_initial_size_values(shape_options_A_)
-                 : make_initial_size_values(shape_options_B_);
-      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},
-              expiration_times.at(gsl::at(shape_names, i)));
-      result[gsl::at(size_names, i)] =
-          std::make_unique<FunctionsOfTime::PiecewisePolynomial<3>>(
-              initial_time_, initial_size_values,
-              expiration_times.at(gsl::at(size_names, i)));
+    result[shape_name] =
+        std::make_unique<FunctionsOfTime::PiecewisePolynomial<2>>(
+            initial_time_,
+            std::array<DataVector, 3>{std::move(shape_func), shape_zeros,
+                                      shape_zeros},
+            expiration_times.at(shape_name));
+    result[size_name] =
+        std::make_unique<FunctionsOfTime::PiecewisePolynomial<3>>(
+            initial_time_,
+            std::array<DataVector, 4>{{std::move(size_func),
+                                       {shape_options.initial_size_values[1]},
+                                       {shape_options.initial_size_values[2]},
+                                       {0.0}}},
+            expiration_times.at(size_name));
+  };
+  if (shape_options_A_.has_value()) {
+    if (not inner_radius_A.has_value()) {
+      ERROR(
+          "A shape map was specified for object A, but no inner radius is "
+          "available. The object must be enclosed by a sphere.");
     }
+    build_shape_and_size_fot(shape_options_A_.value(), *inner_radius_A,
+                             gsl::at(shape_names, 0), gsl::at(size_names, 0));
+  }
+  if (shape_options_B_.has_value()) {
+    if (not inner_radius_A.has_value()) {
+      ERROR(
+          "A shape map was specified for object A, but no inner radius is "
+          "available. The object must be enclosed by a sphere.");
+    }
+    build_shape_and_size_fot(shape_options_B_.value(), *inner_radius_B,
+                             gsl::at(shape_names, 1), gsl::at(size_names, 1));
   }
 
   return result;
@@ -420,13 +464,7 @@ TimeDependentMapOptions<IsCylindrical>::grid_to_inertial_map(
     } else if (rotation_map_.has_value()) {
       return std::make_unique<detail::gi_map<Rotation>>(rotation_map_.value());
     } else {
-      ERROR(
-          "Requesting grid to inertial map without a distorted frame and "
-          "without a Rotation or Expansion map for object "
-          << Object
-          << ". This means there are no time dependent maps. If you don't want "
-             "time dependent maps, specify 'None' for TimeDependentMapOptions. "
-             "Otherwise specify at least one time dependent map.");
+      return nullptr;
     }
   }
 }

src/Domain/Creators/BinaryCompactObjectHelpers.hpp

@@ -16,6 +16,7 @@
 #include "Domain/CoordinateMaps/TimeDependent/Rotation.hpp"
 #include "Domain/CoordinateMaps/TimeDependent/Shape.hpp"
 #include "Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/ShapeMapTransitionFunction.hpp"
+#include "Domain/Creators/SphereTimeDependentMaps.hpp"
 #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
 #include "Domain/Structure/ObjectLabel.hpp"
 #include "Options/Auto.hpp"
@@ -68,6 +69,15 @@ struct ShapeMapOptions {
         "zero."};
   };
 
+  struct InitialValues {
+    using type =
+        Options::Auto<std::variant<sphere::KerrSchildFromBoyerLindquist>,
+                      sphere::Spherical>;
+    static constexpr Options::String help = {
+        "Initial Ylm coefficients for the shape map. Specify 'Spherical' for "
+        "all coefficients to be initialized to zero."};
+  };
+
   struct SizeInitialValues {
     using type = std::array<double, 3>;
     static constexpr Options::String help = {
@@ -82,13 +92,15 @@ struct ShapeMapOptions {
         "be 0 at the outer radius of the inner sphere around the object"};
   };
 
-  using common_options = tmpl::list<LMax, SizeInitialValues>;
+  using common_options = tmpl::list<LMax, InitialValues, SizeInitialValues>;
 
   using options = tmpl::conditional_t<
       IsCylindrical, common_options,
       tmpl::push_back<common_options, TransitionEndsAtCube>>;
 
   size_t l_max{};
+  std::optional<std::variant<sphere::KerrSchildFromBoyerLindquist>>
+      initial_values{};
   std::array<double, 3> initial_size_values{};
   bool transition_ends_at_cube{false};
 };
@@ -289,7 +301,9 @@ struct TimeDependentMapOptions {
    */
   std::unordered_map<std::string,
                      std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>
-  create_functions_of_time(const std::unordered_map<std::string, double>&
+  create_functions_of_time(const std::optional<double> inner_radius_A,
+                           const std::optional<double> inner_radius_B,
+                           const std::unordered_map<std::string, double>&
                                initial_expiration_times) const;
 
   /*!

src/Domain/Creators/CylindricalBinaryCompactObject.cpp

@@ -1096,7 +1096,7 @@ CylindricalBinaryCompactObject::functions_of_time(
     const {
   return time_dependent_options_.has_value()
              ? time_dependent_options_->create_functions_of_time(
-                   initial_expiration_times)
+                   radius_A_, radius_B_, initial_expiration_times)
              : std::unordered_map<
                    std::string,
                    std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>{};

support/Pipelines/Bbh/Inspiral.yaml

@@ -83,10 +83,12 @@ DomainCreator:
         InitialAngularVelocity: [0.0, 0.0, {{ InitialAngularVelocity }}]
       ShapeMapA:
         LMax: &LMax 10
+        InitialValues: Spherical
         SizeInitialValues: [0.0, 0.0, 0.0]
         TransitionEndsAtCube: true
       ShapeMapB:
         LMax: *LMax
+        InitialValues: Spherical
         SizeInitialValues: [0.0, 0.0, 0.0]
         TransitionEndsAtCube: true
 

tests/InputFiles/CurvedScalarWave/WorldtubeKerrSchild.yaml

@@ -75,10 +75,12 @@ DomainCreator:
         InitialAngularVelocity: [0., 0., 0.08944271909999159]
       ShapeMapA:
         LMax: 8
+        InitialValues: Spherical
         SizeInitialValues: [0., 0., 0.]
         TransitionEndsAtCube: false
       ShapeMapB:
         LMax: 8
+        InitialValues: Spherical
         SizeInitialValues: [0., 0., 0.]
         TransitionEndsAtCube: false
 

tests/InputFiles/ExportCoordinates/InputTimeDependent3D.yaml

@@ -72,10 +72,12 @@ DomainCreator:
         InitialAngularVelocity: [0.0, 0.0, 1.5264577062000000e-02]
       ShapeMapA:
         LMax: 8
+        InitialValues: Spherical
         SizeInitialValues: [0.0, 0.0, 0.0]
         TransitionEndsAtCube: true
       ShapeMapB:
         LMax: 8
+        InitialValues: Spherical
         SizeInitialValues: [0.0, 0.0, 0.0]
         TransitionEndsAtCube: true
 

tests/InputFiles/GeneralizedHarmonic/BinaryBlackHole.yaml

@@ -81,10 +81,12 @@ DomainCreator:
         InitialAngularVelocity: [0.0, 0.0, 0.0]
       ShapeMapA:
         LMax: &LMax 10
+        InitialValues: Spherical
         SizeInitialValues: [0.0, 0.0, 0.0]
         TransitionEndsAtCube: true
       ShapeMapB:
         LMax: *LMax
+        InitialValues: Spherical
         SizeInitialValues: [0.0, 0.0, 0.0]
         TransitionEndsAtCube: true
 

tests/InputFiles/GeneralizedHarmonic/CylindricalBinaryBlackHole.yaml

@@ -63,9 +63,11 @@ DomainCreator:
         InitialAngularVelocity: [0.0, 0.0, 1.5264577062000000e-02]
       ShapeMapA:
         LMax: &LMax 10
+        InitialValues: Spherical
         SizeInitialValues: [0.0, 0.0, 0.0]
       ShapeMapB:
         LMax: *LMax
+        InitialValues: Spherical
         SizeInitialValues: [0.0, 0.0, 0.0]
 
 EventsAndDenseTriggers:

tests/Unit/Domain/Creators/Test_BinaryCompactObject.cpp

@@ -368,11 +368,13 @@ std::string create_option_string(
             "      InitialAngularVelocity: [0.0, 0.0, -0.2]\n"s +
                 (excise_A ? "    ShapeMapA:\n"
                             "      LMax: 8\n"
+                            "      InitialValues: Spherical\n"
                             "      SizeInitialValues: [0.0, -0.1, 0.01]\n"
                             "      TransitionEndsAtCube: false\n"s
                           : "    ShapeMapA: None\n"s) +
                 (excise_B ? "    ShapeMapB:\n"
                             "      LMax: 8\n"
+                            "      InitialValues: Spherical\n"
                             "      SizeInitialValues: [0.0, -0.2, 0.02]\n"
                             "      TransitionEndsAtCube: true"s
                           : "    ShapeMapB: None"s)

tests/Unit/Domain/Creators/Test_BinaryCompactObjectHelpers.cpp

@@ -134,18 +134,20 @@ void test(const bool include_expansion, const bool include_rotation,
   const size_t l_max_A = 8;
   if (include_shape_a) {
     shape_map_a_options =
-        IsCylindrical
-            ? ShapeMapAOptions<IsCylindrical>{l_max_A, size_A_values}
-            : ShapeMapAOptions<IsCylindrical>{l_max_A, size_A_values, true};
+        IsCylindrical ? ShapeMapAOptions<IsCylindrical>{l_max_A, std::nullopt,
+                                                        size_A_values}
+                      : ShapeMapAOptions<IsCylindrical>{l_max_A, std::nullopt,
+                                                        size_A_values, true};
   }
 
   const std::array<double, 3> size_B_values{-0.001, -0.02, -0.3};
   const size_t l_max_B = 10;
   if (include_shape_b) {
     shape_map_b_options =
-        IsCylindrical
-            ? ShapeMapBOptions<IsCylindrical>{l_max_B, size_B_values}
-            : ShapeMapBOptions<IsCylindrical>{l_max_B, size_B_values, false};
+        IsCylindrical ? ShapeMapBOptions<IsCylindrical>{l_max_B, std::nullopt,
+                                                        size_B_values}
+                      : ShapeMapBOptions<IsCylindrical>{l_max_B, std::nullopt,
+                                                        size_B_values, false};
   }
 
   const double initial_time = 1.5;
@@ -238,7 +240,7 @@ void test(const bool include_expansion, const bool include_rotation,
           gsl::at(TimeDependentMapOptions<IsCylindrical>::shape_names, 1))};
 
   const auto functions_of_time =
-      time_dep_options.create_functions_of_time(expiration_times);
+      time_dep_options.create_functions_of_time(0.8, 0.5, expiration_times);
 
   if (include_expansion) {
     CHECK(functions_of_time.count(
@@ -413,22 +415,16 @@ void test(const bool include_expansion, const bool include_rotation,
 
     if ((not include_rotation) and (not include_expansion) and
         (not is_excised(excise_A))) {
-      CHECK_THROWS_WITH(
-          time_dep_options
-              .template grid_to_inertial_map<domain::ObjectLabel::A>(excise_A),
-          Catch::Matchers::ContainsSubstring(
-              "Requesting grid to inertial map without a distorted frame and "
-              "without a Rotation or Expansion map for object"));
+      CHECK(time_dep_options
+                .template grid_to_inertial_map<domain::ObjectLabel::A>(
+                    excise_A) == nullptr);
       continue;
     }
     if ((not include_rotation) and (not include_expansion) and
         (not is_excised(excise_B))) {
-      CHECK_THROWS_WITH(
-          time_dep_options
-              .template grid_to_inertial_map<domain::ObjectLabel::B>(excise_B),
-          Catch::Matchers::ContainsSubstring(
-              "Requesting grid to inertial map without a distorted frame and "
-              "without a Rotation or Expansion map for object"));
+      CHECK(time_dep_options
+                .template grid_to_inertial_map<domain::ObjectLabel::B>(
+                    excise_B) == nullptr);
       continue;
     }