Add YlmsFromSpEC to shape time dependent options

.clang-tidy

@@ -55,6 +55,8 @@ Checks: '*,
 -llvm-qualified-auto,
 # We are not developing LLVM libc,
 -llvmlibc-*,
+# Triggers when 'l' is used because it's "too similar" to 'I'
+-misc-confusable-identifiers,
 # thinks constexpr variables in header files cause ODR violations,
 -misc-definitions-in-headers,
 # false positives,

src/Domain/Creators/TimeDependentOptions/ShapeMap.cpp

@@ -4,17 +4,25 @@
 #include "Domain/Creators/TimeDependentOptions/ShapeMap.hpp"
 
 #include <array>
+#include <cmath>
+#include <fstream>
+#include <istream>
+#include <limits>
+#include <sstream>
 #include <string>
 #include <utility>
 #include <variant>
 
 #include "DataStructures/DataVector.hpp"
+#include "DataStructures/ModalVector.hpp"
 #include "DataStructures/Tensor/Tensor.hpp"
 #include "Domain/Structure/ObjectLabel.hpp"
+#include "FromVolumeFile.hpp"
 #include "NumericalAlgorithms/SphericalHarmonics/Spherepack.hpp"
 #include "NumericalAlgorithms/SphericalHarmonics/SpherepackIterator.hpp"
 #include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
 #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrHorizon.hpp"
+#include "Utilities/EqualWithinRoundoff.hpp"
 #include "Utilities/GenerateInstantiations.hpp"
 #include "Utilities/StdArrayHelpers.hpp"
 
@@ -37,6 +45,16 @@ YlmsFromFile::YlmsFromFile(std::string h5_filename_in,
       set_l1_coefs_to_zero(set_l1_coefs_to_zero_in),
       check_frame(check_frame_in) {}
 
+YlmsFromSpEC::YlmsFromSpEC() = default;
+YlmsFromSpEC::YlmsFromSpEC(std::string dat_filename_in,
+                           const double match_time_in,
+                           const std::optional<double> match_time_epsilon_in,
+                           bool set_l1_coefs_to_zero_in)
+    : dat_filename(std::move(dat_filename_in)),
+      match_time(match_time_in),
+      match_time_epsilon(match_time_epsilon_in),
+      set_l1_coefs_to_zero(set_l1_coefs_to_zero_in) {}
+
 template <bool IncludeTransitionEndsAtCube, domain::ObjectLabel Object>
 std::pair<std::array<DataVector, 3>, std::array<DataVector, 4>>
 initial_shape_and_size_funcs(
@@ -65,7 +83,7 @@ initial_shape_and_size_funcs(
       shape_funcs[0] = ylm.phys_to_spec(radial_distortion);
       // Transform from SPHEREPACK to actual Ylm for size func
       size_funcs[0][0] = shape_funcs[0][0] * sqrt(0.5 * M_PI);
-      // Set l=0 for shape map to 0 because size is going to be used
+      // Set l=0 for shape map to 0 because size control will adjust l=0
       shape_funcs[0][0] = 0.0;
     } else if (std::holds_alternative<YlmsFromFile>(
                    shape_options.initial_values.value())) {
@@ -98,14 +116,125 @@ initial_shape_and_size_funcs(
         // 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);
-        // Set l=0 for shape map to 0 because size is going to be used
+        // 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) {
           for (int m = -1; m <= 1; m++) {
             gsl::at(shape_funcs, i)[iter.set(1_st, m)()] = 0.0;
           }
         }
       }
+    } else if (std::holds_alternative<YlmsFromSpEC>(
+                   shape_options.initial_values.value())) {
+      const auto& spec_option =
+          std::get<YlmsFromSpEC>(shape_options.initial_values.value());
+      const std::string& dat_filename = spec_option.dat_filename;
+      const double match_time = spec_option.match_time;
+      const double match_time_epsilon =
+          spec_option.match_time_epsilon.value_or(1e-12);
+      const bool set_l1_coefs_to_zero = spec_option.set_l1_coefs_to_zero;
+
+      std::ifstream dat_file(dat_filename);
+      if (not dat_file.is_open()) {
+        ERROR("Unable to open SpEC dat file " << dat_filename);
+      }
+      std::string line{};
+      size_t total_col = 0;
+      std::optional<size_t> l_max{};
+      std::array<double, 3> center{};
+      ModalVector coefficients{};
+      // This will be actually set below
+      ylm::SpherepackIterator file_iter{2, 2};
+
+      // We have to parse the dat file manually
+      while (std::getline(dat_file, line)) {
+        // Avoid comment lines. The SpEC file puts the legend in comments at the
+        // top of the file, so we count how many columns the dat file has based
+        // on the number of comment lines that are the legend (ends in ')')
+        if (line.starts_with("#")) {
+          if (line.starts_with("# [") and line.ends_with(")")) {
+            ++total_col;
+          }
+          continue;
+        }
+
+        std::stringstream ss(line);
+
+        double time = 0.0;
+        ss >> time;
+
+        // Set scale to current time plus 1 just in case time = 0
+        if (not equal_within_roundoff(time, match_time, match_time_epsilon,
+                                      time + 1.0)) {
+          continue;
+        }
+
+        if (l_max.has_value()) {
+          ERROR("Found more than one time in the SpEC dat file "
+                << dat_filename << " that is within a relative epsilon of "
+                << match_time_epsilon << " of the time requested " << time);
+        }
+
+        // Casting to an integer floors a double, so we add 0.5 before we take
+        // the sqrt to avoid any rounding issues
+        const auto l_max_plus_one =
+            static_cast<size_t>(sqrt(static_cast<double>(total_col) + 0.5));
+        if (l_max_plus_one == 0) {
+          ERROR(
+              "Invalid l_max from SpEC dat file. l_max + 1 was computed to be "
+              "0");
+        }
+        l_max = l_max_plus_one - 1;
+
+        ss >> center[0];
+        ss >> center[1];
+        ss >> center[2];
+
+        coefficients.destructive_resize(
+            ylm::Spherepack::spectral_size(l_max.value(), l_max.value()));
+
+        file_iter = ylm::SpherepackIterator{l_max.value(), l_max.value()};
+
+        for (int l = 0; l <= static_cast<int>(l_max.value()); l++) {
+          for (int m = -l; m <= l; m++) {
+            ss >> coefficients[file_iter.set(static_cast<size_t>(l), m)()];
+          }
+        }
+      }
+
+      if (not l_max.has_value()) {
+        ERROR_NO_TRACE("Unable to find requested time "
+                       << time << " within an epsilon of " << match_time_epsilon
+                       << " in SpEC dat file " << dat_filename);
+      }
+
+      const ylm::Strahlkorper<Frame::Inertial> file_strahlkorper{
+          l_max.value(), l_max.value(), coefficients, center};
+      const ylm::Strahlkorper<Frame::Inertial> this_strahlkorper{
+          shape_options.l_max, 1.0, std::array{0.0, 0.0, 0.0}};
+      ylm::SpherepackIterator iter{shape_options.l_max, shape_options.l_max};
+
+      shape_funcs[0] =
+          -1.0 * file_strahlkorper.ylm_spherepack().prolong_or_restrict(
+                     file_strahlkorper.coefficients(),
+                     this_strahlkorper.ylm_spherepack());
+      // Transform from SPHEREPACK to actual Ylm for size func
+      size_funcs[0][0] = shape_funcs[0][0] * sqrt(0.5 * M_PI);
+      // Set l=0 for shape map to 0 because size control will adjust l=0
+      shape_funcs[0][0] = 0.0;
+      if (set_l1_coefs_to_zero) {
+        for (int m = -1; m <= 1; m++) {
+          shape_funcs[0][iter.set(1_st, m)()] = 0.0;
+        }
+      }
+    } else if (std::holds_alternative<FromVolumeFile<names::ShapeSize<Object>>>(
+                   shape_options.initial_values.value())) {
+      const auto& volume_file_options =
+          std::get<FromVolumeFile<names::ShapeSize<Object>>>(
+              shape_options.initial_values.value());
+
+      shape_funcs = volume_file_options.shape_values;
+      size_funcs = volume_file_options.size_values;
     }
   }
 

src/Domain/Creators/TimeDependentOptions/ShapeMap.hpp

@@ -11,6 +11,7 @@
 #include <variant>
 
 #include "DataStructures/DataVector.hpp"
+#include "Domain/Creators/TimeDependentOptions/FromVolumeFile.hpp"
 #include "Domain/Structure/ObjectLabel.hpp"
 #include "NumericalAlgorithms/SphericalHarmonics/IO/ReadSurfaceYlm.hpp"
 #include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
@@ -129,6 +130,54 @@ struct YlmsFromFile {
   bool check_frame{true};
 };
 
+struct YlmsFromSpEC {
+  struct DatFilename {
+    using type = std::string;
+    static constexpr Options::String help =
+        "Name of the SpEC dat file holding the coefficients. Note that this "
+        "isn't a Dat file within an H5 file. This must be an actual '.dat' "
+        "file on disk.";
+  };
+
+  struct MatchTime {
+    using type = double;
+    static constexpr Options::String help =
+        "Time in the H5File to get the coefficients at. Will likely be the "
+        "same as the initial time";
+  };
+
+  struct MatchTimeEpsilon {
+    using type = Options::Auto<double>;
+    static constexpr Options::String help =
+        "Look for times in the H5File within this epsilon of the match time. "
+        "This is to avoid having to know the exact time to all digits. Default "
+        "is 1e-12.";
+  };
+
+  struct SetL1CoefsToZero {
+    using type = bool;
+    static constexpr Options::String help =
+        "Whether to set the L=1 coefs to zero or not. This may be desirable "
+        "because L=1 is degenerate with a translation of the BH.";
+  };
+
+  using options =
+      tmpl::list<DatFilename, MatchTime, MatchTimeEpsilon, SetL1CoefsToZero>;
+
+  static constexpr Options::String help = {
+      "Read the Y_lm coefficients of a Strahlkorper from file and use those to "
+      "initialize the coefficients of a shape map."};
+  YlmsFromSpEC();
+  YlmsFromSpEC(std::string dat_filename_in, double match_time_in,
+               std::optional<double> match_time_epsilon_in,
+               bool set_l1_coefs_to_zero_in);
+
+  std::string dat_filename{};
+  double match_time{};
+  std::optional<double> match_time_epsilon{};
+  bool set_l1_coefs_to_zero{};
+};
+
 /*!
  * \brief Class to be used as an option for initializing shape map coefficients.
  *
@@ -154,9 +203,10 @@ struct ShapeMapOptions {
   };
 
   struct InitialValues {
-    using type =
-        Options::Auto<std::variant<KerrSchildFromBoyerLindquist, YlmsFromFile>,
-                      Spherical>;
+    using type = Options::Auto<
+        std::variant<KerrSchildFromBoyerLindquist, YlmsFromFile, YlmsFromSpEC,
+                     FromVolumeFile<names::ShapeSize<Object>>>,
+        Spherical>;
     static constexpr Options::String help = {
         "Initial Ylm coefficients for the shape map. Specify 'Spherical' for "
         "all coefficients to be initialized to zero."};
@@ -187,7 +237,9 @@ struct ShapeMapOptions {
   ShapeMapOptions() = default;
   ShapeMapOptions(
       size_t l_max_in,
-      std::optional<std::variant<KerrSchildFromBoyerLindquist, YlmsFromFile>>
+      std::optional<
+          std::variant<KerrSchildFromBoyerLindquist, YlmsFromFile, YlmsFromSpEC,
+                       FromVolumeFile<names::ShapeSize<Object>>>>
           initial_values_in,
       std::optional<std::array<double, 3>> initial_size_values_in =
           std::nullopt,
@@ -198,7 +250,9 @@ struct ShapeMapOptions {
         transition_ends_at_cube(transition_ends_at_cube_in) {}
 
   size_t l_max{};
-  std::optional<std::variant<KerrSchildFromBoyerLindquist, YlmsFromFile>>
+  std::optional<
+      std::variant<KerrSchildFromBoyerLindquist, YlmsFromFile, YlmsFromSpEC,
+                   FromVolumeFile<names::ShapeSize<Object>>>>
       initial_values{};
   std::optional<std::array<double, 3>> initial_size_values{};
   bool transition_ends_at_cube{false};