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]),

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"
@@ -145,34 +146,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 =
+            inner_radius -
+            get(gr::Solutions::kerr_schild_radius_from_boyer_lindquist(
+                inner_radius, ylm.theta_phi_points(), mass_and_spin.mass,
+                mass_and_spin.spin));
+        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_radii_[0].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_radii_[0],
+                             shape_names[0], size_names[0]);
+  }
+  if (shape_options_B_.has_value()) {
+    if (not inner_radii_[1].has_value()) {
+      ERROR(
+          "A shape map was specified for object B, but no inner radius is "
+          "available. The object must be enclosed by a sphere.");
+    }
+    build_shape_and_size_fot(shape_options_B_.value(), *inner_radii_[1],
+                             shape_names[1], size_names[1]);
   }
 
   return result;
@@ -206,6 +248,8 @@ void TimeDependentMapOptions<IsCylindrical>::build_maps(
             << ", but no time dependent map options were specified "
                "for that object.");
       }
+      // Store the inner radii for creating functions of time
+      gsl::at(inner_radii_, i) = radii[0];
 
       const size_t initial_l_max = i == 0 ? shape_options_A_.value().l_max
                                           : shape_options_B_.value().l_max;
@@ -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};
 };
@@ -394,6 +406,7 @@ struct TimeDependentMapOptions {
   std::optional<RotationMapOptions> rotation_options_{};
   std::optional<ShapeMapOptions<domain::ObjectLabel::A>> shape_options_A_{};
   std::optional<ShapeMapOptions<domain::ObjectLabel::B>> shape_options_B_{};
+  std::array<std::optional<double>, 2> inner_radii_{};
 
   // Maps
   std::optional<Expansion> expansion_map_{};

src/Domain/Creators/SphereTimeDependentMaps.cpp

@@ -78,10 +78,10 @@ TimeDependentMapOptions::create_functions_of_time(
       const auto& mass_and_spin = std::get<KerrSchildFromBoyerLindquist>(
           shape_map_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;
+          inner_radius -
+          get(gr::Solutions::kerr_schild_radius_from_boyer_lindquist(
+              inner_radius, ylm.theta_phi_points(), mass_and_spin.mass,
+              mass_and_spin.spin));
       shape_func = ylm.phys_to_spec(radial_distortion);
       // Transform from SPHEREPACK to actual Ylm for size func
       size_func[0] = shape_func[0] * sqrt(0.5 * M_PI);

src/Domain/Creators/TimeDependence/Shape.cpp

@@ -178,8 +178,8 @@ Shape<Label>::functions_of_time(const std::unordered_map<std::string, double>&
       std::array<DataVector, 4>{
           // Size holds the *actual* \lambda_00 spherical harmonic coefficient,
           // but shape holds Spherepack coefficients so we must convert between
-          // the two. Need to multiply lambda_00 by sqrt(2/pi)
-          {{M_SQRT1_2 * M_2_SQRTPI * radial_distortion_coefs[0]},
+          // the two. Need to multiply lambda_00 by sqrt(pi/2)
+          {{sqrt(0.5 * M_PI) * radial_distortion_coefs[0]},
            zeros_size,
            zeros_size,
            zeros_size}},

support/Pipelines/Bbh/InitialData.py

@@ -16,6 +16,11 @@
 ID_INPUT_FILE_TEMPLATE = Path(__file__).parent / "InitialData.yaml"
 
 
+def L1_distance(m1, m2, separation):
+    """Distance of the L1 Lagrangian point from m1, in Newtonian gravity"""
+    return separation * (0.5 - 0.227 * np.log10(m2 / m1))
+
+
 def id_parameters(
     mass_ratio: float,
     dimensionless_spin_a: Sequence[float],
@@ -49,6 +54,7 @@ def id_parameters(
     M_B = 1.0 / (1.0 + mass_ratio)
     x_A = separation / (1.0 + mass_ratio)
     x_B = x_A - separation
+    # Spins
     chi_A = np.asarray(dimensionless_spin_a)
     r_plus_A = M_A * (1.0 + np.sqrt(1 - np.dot(chi_A, chi_A)))
     Omega_A = -0.5 * chi_A / r_plus_A
@@ -57,13 +63,18 @@ def id_parameters(
     r_plus_B = M_B * (1.0 + np.sqrt(1 - np.dot(chi_B, chi_B)))
     Omega_B = -0.5 * chi_B / r_plus_B
     Omega_B[2] += orbital_angular_velocity
+    # Falloff widths of superposition
+    L1_dist_A = L1_distance(M_A, M_B, separation)
+    L1_dist_B = separation - L1_dist_A
+    falloff_width_A = 3.0 / 5.0 * L1_dist_A
+    falloff_width_B = 3.0 / 5.0 * L1_dist_B
     return {
         "MassRight": M_A,
         "MassLeft": M_B,
         "XRight": x_A,
         "XLeft": x_B,
-        "ExcisionRadiusRight": 0.89 * 2.0 * M_A,
-        "ExcisionRadiusLeft": 0.89 * 2.0 * M_B,
+        "ExcisionRadiusRight": 0.89 * r_plus_A,
+        "ExcisionRadiusLeft": 0.89 * r_plus_B,
         "OrbitalAngularVelocity": orbital_angular_velocity,
         "RadialExpansionVelocity": radial_expansion_velocity,
         "DimensionlessSpinRight_x": chi_A[0],
@@ -78,6 +89,8 @@ def id_parameters(
         "HorizonRotationLeft_x": Omega_B[0],
         "HorizonRotationLeft_y": Omega_B[1],
         "HorizonRotationLeft_z": Omega_B[2],
+        "FalloffWidthRight": falloff_width_A,
+        "FalloffWidthLeft": falloff_width_B,
         # Resolution
         "L": refinement_level,
         "P": polynomial_order,

support/Pipelines/Bbh/InitialData.yaml

@@ -31,24 +31,26 @@ Background: &background
     XCoords: [&x_left {{ XLeft }}, &x_right {{ XRight }}]
     ObjectLeft: &kerr_left
       KerrSchild:
-        Mass: {{ MassLeft }}
-        Spin:
+        Mass: &mass_left {{ MassLeft }}
+        Spin: &spin_left
           - {{ DimensionlessSpinLeft_x }}
           - {{ DimensionlessSpinLeft_y }}
           - {{ DimensionlessSpinLeft_z }}
         Center: [0., 0., 0.]
     ObjectRight: &kerr_right
       KerrSchild:
-        Mass: {{ MassRight }}
-        Spin:
+        Mass: &mass_right {{ MassRight }}
+        Spin: &spin_right
           - {{ DimensionlessSpinRight_x }}
           - {{ DimensionlessSpinRight_y }}
           - {{ DimensionlessSpinRight_z }}
         Center: [0., 0., 0.]
     AngularVelocity: {{ OrbitalAngularVelocity }}
     Expansion: {{ RadialExpansionVelocity }}
     LinearVelocity: [0., 0., 0.]
-    FalloffWidths: [4.8, 4.8]
+    FalloffWidths:
+      - {{ FalloffWidthLeft }}
+      - {{ FalloffWidthRight }}
 
 InitialGuess: *background
 
@@ -98,8 +100,8 @@ DomainCreator:
       ObjectBShell:     [{{ L }}, {{ L }}, {{ L }}]
       ObjectACube:      [{{ L }}, {{ L }}, {{ L }}]
       ObjectBCube:      [{{ L }}, {{ L }}, {{ L }}]
-      Envelope:         [{{ L }}, {{ L }}, {{ L }}]
-      OuterShell:       [{{ L }}, {{ L }}, {{ L + 2}}]
+      Envelope:         [{{ L }}, {{ L }}, {{ L + 1 }}]
+      OuterShell:       [{{ L }}, {{ L }}, {{ L + 2 }}]
     # This p-refinement represents a crude manual optimization of the domain. We
     # will need AMR to optimize the domain further.
     InitialGridPoints:
@@ -109,7 +111,24 @@ DomainCreator:
       ObjectBCube:    [{{ P + 1}}, {{ P + 1}}, {{ P + 2}}]
       Envelope:       [{{ P + 1}}, {{ P + 1}}, {{ P + 1}}]
       OuterShell:     [{{ P + 1}}, {{ P + 1}}, {{ P + 1}}]
-    TimeDependentMaps: None
+    TimeDependentMaps:
+      InitialTime: 0.
+      ExpansionMap: None
+      RotationMap: None
+      ShapeMapA:
+        LMax: 20
+        InitialValues:
+          Mass: *mass_right
+          Spin: *spin_right
+        SizeInitialValues: [0, 0, 0]
+        TransitionEndsAtCube: False
+      ShapeMapB:
+        LMax: 20
+        InitialValues:
+          Mass: *mass_left
+          Spin: *spin_left
+        SizeInitialValues: [0, 0, 0]
+        TransitionEndsAtCube: False
 
 Amr:
   Verbosity: Verbose
@@ -150,7 +169,7 @@ LinearSolver:
     ConvergenceCriteria:
       MaxIterations: 100
       RelativeResidual: 1.e-3
-      AbsoluteResidual: 1.e-10
+      AbsoluteResidual: 1.e-9
     Verbosity: Quiet
 
   Multigrid:
@@ -188,7 +207,7 @@ RadiallyCompressedCoordinates:
   Compression: *outer_shell_distribution
 
 EventsAndTriggers:
-  - Trigger: HasConverged
+  - Trigger: Always
     Events:
       - ObserveFields:
           SubfileName: VolumeData

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