Allow prolongation/restriction for shape FoT from vol file

src/Domain/Creators/TimeDependentOptions/ShapeMap.cpp

@@ -9,6 +9,7 @@
 #include <istream>
 #include <limits>
 #include <memory>
+#include <optional>
 #include <sstream>
 #include <string>
 #include <utility>
@@ -27,6 +28,7 @@
 #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrHorizon.hpp"
 #include "Utilities/EqualWithinRoundoff.hpp"
 #include "Utilities/GenerateInstantiations.hpp"
+#include "Utilities/Gsl.hpp"
 #include "Utilities/StdArrayHelpers.hpp"
 
 namespace domain::creators::time_dependent_options {
@@ -60,21 +62,26 @@ YlmsFromSpEC::YlmsFromSpEC(std::string dat_filename_in,
 
 template <ObjectLabel Object>
 FromVolumeFileShapeSize<Object>::FromVolumeFileShapeSize(
+    const std::optional<size_t>& l_max_in,
     const bool transition_ends_at_cube_in, std::string h5_filename,
     std::string subfile_name)
     : FromVolumeFile(std::move(h5_filename), std::move(subfile_name)),
       transition_ends_at_cube(transition_ends_at_cube_in) {
-  const auto shape_fot_map = retrieve_function_of_time(
-      {std::string{"Shape" + name(Object)}}, std::nullopt);
-  const auto& shape_fot = shape_fot_map.at("Shape" + name(Object));
-
-  const double initial_time = shape_fot->time_bounds()[0];
-  const auto function = shape_fot->func(initial_time);
-
-  // num_components = 2 * (l_max + 1)**2 if l_max == m_max which it is for the
-  // shape map. This is why we can divide by 2 and take the sqrt without
-  // worrying about odd numbers or non-perfect squares
-  l_max = -1 + sqrt(function[0].size() / 2);
+  if (l_max_in.has_value()) {
+    l_max = l_max_in.value();
+  } else {
+    const auto shape_fot_map = retrieve_function_of_time(
+        {std::string{"Shape" + name(Object)}}, std::nullopt);
+    const auto& shape_fot = shape_fot_map.at("Shape" + name(Object));
+
+    const double initial_time = shape_fot->time_bounds()[0];
+    const auto function = shape_fot->func(initial_time);
+
+    // num_components = 2 * (l_max + 1)**2 if l_max == m_max which it is for the
+    // shape map. This is why we can divide by 2 and take the sqrt without
+    // worrying about odd numbers or non-perfect squares
+    l_max = -1 + sqrt(function[0].size() / 2);
+  }
 }
 
 template <bool IncludeTransitionEndsAtCube, domain::ObjectLabel Object>
@@ -130,9 +137,43 @@ FunctionsOfTimeMap get_shape_and_size(
     check_fot.template operator()<2>(shape_name);
     check_fot.template operator()<3>(size_name);
 
-    result[shape_name] =
+    // Not const in case we move it below
+    auto temporary_shape_fot =
         volume_fots.at(shape_name)
             ->create_at_time(initial_time, shape_expiration_time);
+    const auto shape_funcs =
+        temporary_shape_fot->func_and_2_derivs(initial_time);
+
+    const size_t file_l_max = -1 + sqrt(shape_funcs[0].size() / 2);
+
+    // Prolong or restrict if necessary, otherwise just use the exact function
+    // of time from the volume file
+    if (file_l_max != l_max) {
+      ModalVector shape_coefs{};
+      std::array<DataVector, 3> new_shape_funcs{};
+      for (size_t i = 0; i < shape_funcs.size(); i++) {
+        shape_coefs.set_data_ref(
+            const_cast<DataVector&>(gsl::at(shape_funcs, i)).data(),  // NOLINT
+            gsl::at(shape_funcs, i).size());
+        // The centers and radii of the following don't matter. Only their l_max
+        // and their coefficients
+        const ylm::Strahlkorper<Frame::Distorted> file_strahlkorper{
+            file_l_max, file_l_max, shape_coefs, std::array{0.0, 0.0, 0.0}};
+        const ylm::Strahlkorper<Frame::Distorted> new_strahlkorper{
+            l_max, 1.0, std::array{0.0, 0.0, 0.0}};
+
+        gsl::at(new_shape_funcs, i) =
+            file_strahlkorper.ylm_spherepack().prolong_or_restrict(
+                file_strahlkorper.coefficients(),
+                new_strahlkorper.ylm_spherepack());
+      }
+
+      result[shape_name] =
+          std::make_unique<domain::FunctionsOfTime::PiecewisePolynomial<2>>(
+              initial_time, std::move(new_shape_funcs), shape_expiration_time);
+    } else {
+      result[shape_name] = std::move(temporary_shape_fot);
+    }
     result[size_name] = volume_fots.at(size_name)->create_at_time(
         initial_time, size_expiration_time);
 

src/Domain/Creators/TimeDependentOptions/ShapeMap.hpp

@@ -193,15 +193,24 @@ struct TransitionEndsAtCube {
  * \brief Specialized version of `FromVolumeFile` for the shape map
  *
  * \details This is needed because the regular `FromVolumeFile` doesn't have
- * options for domain settings like `TransitionEndsAtCube`.
+ * options for domain settings like `TransitionEndsAtCube` or `LMax`.
  */
 template <ObjectLabel Object>
 struct FromVolumeFileShapeSize : public FromVolumeFile {
-  using options =
-      tmpl::push_front<FromVolumeFile::options, detail::TransitionEndsAtCube>;
+ public:
+  struct LMax {
+    using type = Options::Auto<size_t>;
+    static constexpr Options::String help = {
+        "LMax used for the number of spherical harmonic coefficients of the "
+        "distortion map. If set to 'Auto', will use the LMax from the shape "
+        "function of time in the volume file."};
+  };
+  using options = tmpl::push_front<FromVolumeFile::options, LMax,
+                                   detail::TransitionEndsAtCube>;
 
   FromVolumeFileShapeSize() = default;
-  FromVolumeFileShapeSize(bool transition_ends_at_cube_in,
+  FromVolumeFileShapeSize(const std::optional<size_t>& l_max_in,
+                          bool transition_ends_at_cube_in,
                           std::string h5_filename, std::string subfile_name);
 
   size_t l_max{};

support/Pipelines/Bbh/Inspiral.yaml

@@ -139,10 +139,12 @@ DomainCreator:
       ShapeMapA:
         H5Filename: "{{ FotFilename }}"
         SubfileName: "VolumeData"
+        LMax: Auto
         TransitionEndsAtCube: True
       ShapeMapB:
         H5Filename: "{{ FotFilename }}"
         SubfileName: "VolumeData"
+        LMax: Auto
         TransitionEndsAtCube: True
 {% else %}
     TimeDependentMaps:

tests/Unit/Domain/Creators/TimeDependentOptions/Test_ShapeMap.cpp

@@ -299,6 +299,7 @@ void test_funcs(const gsl::not_null<Generator*> generator) {
       const auto shape_map_options = TestHelpers::test_option_tag<
           domain::creators::time_dependent_options::ShapeMapOptions<
               false, domain::ObjectLabel::B>>(
+          "LMax: 6\n"
           "TransitionEndsAtCube: false\n"
           "H5Filename: TotalEclipseOfTheHeart.h5\n"
           "SubfileName: VolumeData\n");
@@ -308,6 +309,10 @@ void test_funcs(const gsl::not_null<Generator*> generator) {
               FromVolumeFileShapeSize<domain::ObjectLabel::B>>(
           shape_map_options.value()));
 
+      CHECK(std::get<FromVolumeFileShapeSize<domain::ObjectLabel::B>>(
+                shape_map_options.value())
+                .l_max == 6);
+
       const FunctionsOfTimeMap shape_and_size = get_shape_and_size(
           shape_map_options.value(), 0.2, 1.1, 1.2, inner_radius);
 
@@ -316,8 +321,32 @@ void test_funcs(const gsl::not_null<Generator*> generator) {
       CHECK(shape_and_size.at(size_name)->time_bounds() ==
             std::array{0.2, 1.2});
 
-      CHECK(shape_and_size.at(shape_name)->func_and_2_derivs(0.2) ==
-            functions_of_time.at(shape_name)->func_and_2_derivs(0.2));
+      const auto shape = shape_and_size.at(shape_name)->func_and_2_derivs(0.2);
+      const auto expected_shape = functions_of_time.at(shape_name)
+                                      ->create_at_time(0.2, 1.1)
+                                      ->func_and_2_derivs(0.2);
+
+      // We restricted from LMax = 8 to LMax = 6
+      const size_t expected_l_max = l_max - 2;
+      const size_t expected_size =
+          ylm::Spherepack::spectral_size(expected_l_max, expected_l_max);
+      ylm::SpherepackIterator file_iter{l_max, l_max};
+      ylm::SpherepackIterator iter{expected_l_max, expected_l_max};
+
+      for (size_t i = 0; i < 3; i++) {
+        CAPTURE(i);
+        CHECK(gsl::at(shape, i).size() == expected_size);
+
+        // Loop pointwise so we only check the coefficients that matter
+        for (size_t l = 0; l <= expected_l_max; l++) {
+          CAPTURE(l);
+          for (int m = -static_cast<int>(l); m <= static_cast<int>(l); m++) {
+            CAPTURE(m);
+            CHECK(gsl::at(shape, i)[iter.set(l, m)()] ==
+                  gsl::at(expected_shape, i)[file_iter.set(l, m)()]);
+          }
+        }
+      }
       CHECK(shape_and_size.at(size_name)->func_and_2_derivs(0.2) ==
             functions_of_time.at(size_name)->func_and_2_derivs(0.2));
     }