From 1274f505b47a13cf2cfbd5d0fae548c262ede9be Mon Sep 17 00:00:00 2001
From: Nils Vu <owls@nilsvu.de>
Date: Wed, 4 Dec 2024 13:56:51 -0800
Subject: [PATCH] Fix shape map loaded from file

The size of the domain didn't match the Strahlkorper in the
file because the loaded l=0, m=0 coefficient was interpreted
as a deformation/scale. Fixed by subtracting the size of the
original sphere.
---
 .../Creators/TimeDependentOptions/ShapeMap.cpp   | 14 +++++++++-----
 .../Creators/TimeDependentOptions/ShapeMap.hpp   |  6 ++++--
 .../TimeDependentOptions/Test_ShapeMap.cpp       | 16 ++++++++++------
 3 files changed, 23 insertions(+), 13 deletions(-)

diff --git a/src/Domain/Creators/TimeDependentOptions/ShapeMap.cpp b/src/Domain/Creators/TimeDependentOptions/ShapeMap.cpp
index a86cc1530783..d0a7209cc8e6 100644
--- a/src/Domain/Creators/TimeDependentOptions/ShapeMap.cpp
+++ b/src/Domain/Creators/TimeDependentOptions/ShapeMap.cpp
@@ -59,7 +59,7 @@ template <bool IncludeTransitionEndsAtCube, domain::ObjectLabel Object>
 std::pair<std::array<DataVector, 3>, std::array<DataVector, 4>>
 initial_shape_and_size_funcs(
     const ShapeMapOptions<IncludeTransitionEndsAtCube, Object>& shape_options,
-    const double inner_radius) {
+    const double deformed_radius) {
   const DataVector shape_zeros{
       ylm::Spherepack::spectral_size(shape_options.l_max, shape_options.l_max),
       0.0};
@@ -76,9 +76,9 @@ initial_shape_and_size_funcs(
       const auto& mass_and_spin = std::get<KerrSchildFromBoyerLindquist>(
           shape_options.initial_values.value());
       const DataVector radial_distortion =
-          inner_radius -
+          deformed_radius -
           get(gr::Solutions::kerr_schild_radius_from_boyer_lindquist(
-              inner_radius, ylm.theta_phi_points(), mass_and_spin.mass,
+              deformed_radius, ylm.theta_phi_points(), mass_and_spin.mass,
               mass_and_spin.spin));
       shape_funcs[0] = ylm.phys_to_spec(radial_distortion);
       // Transform from SPHEREPACK to actual Ylm for size func
@@ -115,7 +115,11 @@ initial_shape_and_size_funcs(
                        this_strahlkorper.ylm_spherepack());
         // Transform from SPHEREPACK to actual Ylm for size func
         gsl::at(size_funcs, i)[0] =
-            gsl::at(shape_funcs, i)[0] * sqrt(0.5 * M_PI);
+            gsl::at(shape_funcs, i)[0] * sqrt(0.5 * M_PI) +
+            // Account for the size of the original sphere, since the shape/size
+            // coefficients are deformations from the original sphere.
+            // The factor 2 sqrt(pi) is 1/Y_00.
+            deformed_radius * 2.0 * sqrt(M_PI);
         // Set l=0 for shape map to 0 because size control will adjust l=0
         gsl::at(shape_funcs, i)[0] = 0.0;
         if (set_l1_coefs_to_zero) {
@@ -259,7 +263,7 @@ initial_shape_and_size_funcs(
   initial_shape_and_size_funcs<INCLUDETRANSITION(data), OBJECT(data)>(     \
       const ShapeMapOptions<INCLUDETRANSITION(data), OBJECT(data)>&        \
           shape_options,                                                   \
-      double inner_radius);
+      double deformed_radius);
 
 GENERATE_INSTANTIATIONS(INSTANTIATE, (true, false),
                         (domain::ObjectLabel::A, domain::ObjectLabel::B,
diff --git a/src/Domain/Creators/TimeDependentOptions/ShapeMap.hpp b/src/Domain/Creators/TimeDependentOptions/ShapeMap.hpp
index 69f1d50d5020..0d181fb2ae59 100644
--- a/src/Domain/Creators/TimeDependentOptions/ShapeMap.hpp
+++ b/src/Domain/Creators/TimeDependentOptions/ShapeMap.hpp
@@ -217,7 +217,9 @@ struct ShapeMapOptions {
     static constexpr Options::String help = {
         "Initial value and two derivatives of the 00 coefficient. Specify "
         "'Auto' to use the 00 coefficient specified in the 'InitialValues' "
-        "option."};
+        "option. If you specify 'Auto', the deformed sphere will match the "
+        "'InitialValues' surface exactly, and the original radius will only "
+        "set the radius of the sphere in the grid frame (before deformation)."};
   };
 
   struct TransitionEndsAtCube {
@@ -262,5 +264,5 @@ template <bool IncludeTransitionEndsAtCube, domain::ObjectLabel Object>
 std::pair<std::array<DataVector, 3>, std::array<DataVector, 4>>
 initial_shape_and_size_funcs(
     const ShapeMapOptions<IncludeTransitionEndsAtCube, Object>& shape_options,
-    double inner_radius);
+    double deformed_radius);
 }  // namespace domain::creators::time_dependent_options
diff --git a/tests/Unit/Domain/Creators/TimeDependentOptions/Test_ShapeMap.cpp b/tests/Unit/Domain/Creators/TimeDependentOptions/Test_ShapeMap.cpp
index 9e49dd61d3ab..5c4ec4c44868 100644
--- a/tests/Unit/Domain/Creators/TimeDependentOptions/Test_ShapeMap.cpp
+++ b/tests/Unit/Domain/Creators/TimeDependentOptions/Test_ShapeMap.cpp
@@ -377,12 +377,16 @@ void test_funcs(const gsl::not_null<Generator*> generator) {
       for (size_t i = 0; i < shape_funcs.size(); i++) {
         CHECK(gsl::at(shape_funcs, i)[0] == 0.0);
       }
-      CHECK(size_funcs ==
-            std::array{DataVector{-1.0 * strahlkorpers[0].coefficients()[0] *
-                                  sqrt(0.5 * M_PI)},
-                       DataVector{-1.0 * strahlkorpers[1].coefficients()[0] *
-                                  sqrt(0.5 * M_PI)},
-                       DataVector{0.0}, DataVector{0.0}});
+      CHECK_ITERABLE_APPROX(
+          size_funcs,
+          (std::array{
+              DataVector{(inner_radius - ylm::Spherepack::average(
+                                             strahlkorpers[0].coefficients())) *
+                         2.0 * sqrt(M_PI)},
+              DataVector{(inner_radius - ylm::Spherepack::average(
+                                             strahlkorpers[1].coefficients())) *
+                         2.0 * sqrt(M_PI)},
+              DataVector{0.0}, DataVector{0.0}}));
     }
 
     if (file_system::check_if_file_exists(test_filename)) {